Hi Jason,
Very useful appreciated hints!
Below my few comments.
On 22 May 2014 23:42, Jason Swails <jason.swails.gmail.com> wrote:
>
> On May 22, 2014, at 3:52 PM, Massimiliano Porrini <
> m.porrini.iecb.u-bordeaux.fr> wrote:
>
> > Hi all,
> >
> > I am attempting to run in vacuo simulations and I would like to know
> which
> > would be
> > the best setting with sander module.
> >
> > As stated at p.60 of Amber14 user guide, igb=6 would correspond to an in
> > vacuo
> > simulation and I quote here what it is written "This option is logically
> > equivalent to
> > setting igb=0 and eedmeth=4, although the implementation (and
> computational
> > efficiency)
> > is not the same".
> >
> > It was not clear to me in which of the 2 cases (igb=6 or ntb=0, igb=0,
> > eedmeth=4)
> > the computational efficiency should improve, therefore I made some tests
> > with the
> > 3 different temperature couplings available in Amber12, here are the
> > results:
> >
> > *igb=6*
> > Berendsen: 113.7 ns/day
> > Andersen(vrand=1000): 70.12 ns/day
> > Langevin(gamma_ln=2): 61.79 ns/day
> >
> > *ntb=0, igb=0, eedmeth=4*
> > Berendsen: 66.48 ns/day
> > Andersen(vrand=1000): 68.94 ns/day
> > Langevin(gamma_ln=2): 72.70 ns/day
> >
> > I ran a simulation of a 758 atoms molecule, on a blade with 8 CPUs dual
> core
> > (hence 16 cores) and hyper-threading (hence 32 threads), however I
> > utilised only 16 cores of such a blade and with the following
> > sander input file (example from the simulation with Berendsen and igb=6):
>
> First a few words of caution about running vacuum dynamics -- the
> parameters (including the standard charge sets and am1-bcc charge
> derivation schemes) are designed to reproduce solvent-induced polarization
> (particularly polar solvents). So you should be careful about how you
> approach vacuum dynamics in general (for instance, dominant protonation
> states are completely different -- zwitterions are unstable in vacuum but
> very stable in solution).
>
Yes, I am aware of this risk, as removing the solvent (the very polar
water) the potential energy surface is totally
different and in principle one should recalculate the whole set of
parameters for the solute.
However, I just want the solute to explore the accessible microstates in
absence of water, and compare
the associated collision cross sections with ion mobility mass spectrometry
experimental values.
> >
> > dyn in vacuo (igb=6)
> > &cntrl
> > imin = 0, irest = 1, ntx = 5,
> > nstlim = 100000, dt = 0.001,
> > ntc = 2, ntf = 2,
> > cut = 1000.0,
> > igb = 6,
> > ntpr = 5000, ntwr = 10000, ntwx = 2000,
> > ntt = 1,
> > temp0 = 300.0,
> > ioutfm = 1
> > /
> >
> > Therefore I would obviously be tempted to utilise Berendsen thermostat
> with
> > igb=6,
> > as it is the fastest one (113.7 ns/day)
>
> I question the use of any thermostat at all. The idea behind a thermostat
> is that there is some "unlimited heat bath at a constant temperature"
> surrounding the system of interest that is able to provide or absorb heat
> energy as needed. This makes far more sense in a condensed phase
> simulation than it does in vacuum or gas. If I was running a gas phase
> simulation, for whatever reason, I would probably use a short time step,
> stringent SHAKE tolerance (if I was using SHAKE at all), and if I didn't
> observe acceptable energy conservation, a Berendsen thermostat with a very
> weak coupling constant (on the order of 10 to 20 /ps) to smooth over
> integration errors.
>
I totally forgot to mention 2 things:
1) the simulation in microcanonical ensemble was in my "to-do" list, I have
just sent this email prematurely and I agree that in vacuo
(no external bath) one should implement pure Newton dynamics, however (and
here it comes the 2nd thing)
2) the simulation in not actually in the "high" vacuum, but there is a
buffer gas at (relatively) low pressure ranges
(with respect to the atmospheric one) and temperature T (which, in my case,
is set equal to 300K), hence I would
pursue the Berendsen coupling choice, with the settings suggested by you
(short time-step, very low SHAKE tolerance
and weak coupling).
> > But, provided that Berendsen coupling is kind of deprecated for several
> > reasons (Amber14
> > user guide, p.290), especially for simulations in implicit solvent, in
> > which a uniform
> > temperature is not guaranteed because collisions with water molecules are
> > missing, do you think that it would be a sensible choice to use Berendsen
> > coupling
> > for in vacuo simulations?
> >
> > My answer would be yes for two reasons:
> > 1) performance efficiency
> > 2) with in vacuo simulations one does not actually want
> > any "random kick" (like it occurs with Andersen or Langevin
> > coupling) that would emulate the collisions with the solvent.
>
> The weak coupling algorithm (for either temperature or pressure) has
> well-documented shortcomings and is fundamentally _known_ not to sample
> from the desired distribution. For one thing, fluctuations are too small,
> and the thermostat doesn't even keep correct partitioning between
> translational, rotational, and vibrational motion (flying block of ice
> effect in explicit solvent). I would personally not use the Berendsen
> thermostat here except as I described above.
>
I was aware of these Berendsen thermostat drawbacks too, this is one main
reason why I wanted any suggestion about
running dynamics in gas phase (not actually vacuum).
>
> >
> > Am I right/wrong?
> > Any comments/hints on this issue would be very appreciated.
> >
> > Many thanks!
> >
> > PS Just to be meticulous (running 1 step of minimisation) I have also
> > checked if the potential functions utilised were
> > actually the same between igb=6 and igb=0, ntb=0, eedmeth=4 settings and
> > they
> > perfectly match :) as you can see below (therefore igb=6 with Berendsen
> is
> > actually
> > the most efficient implementation):
>
> The difference is, as the manual states, in the implementation. I would
> use igb=6 when doing gas-phase energy calculations since the electrostatic
> code used to compute the energy in that case respects the exclusion list
> stored in the topology file whereas the PME code (eedmeth=4, igb=0)
> computes the exclusions from the bond list. This only really matters if
> you add custom exclusions with ParmEd.
>
I am slightly puzzled here, in that I can not use PME where there are no
periodic boundaries
(like for molecules in gas phase).
In fact, if igb>0 then ntb=0 by default (thus sander module does not
implement PME algorithm).
Therefore, if one wants to utilise exactly the same setting, one needs to
impose ntb=0 when [eedmeth=4, igb=0],
otherwise ntb would be equal to 1 by default (igb=0).
As a matter of fact, my simulation with [eedmeth=4, igb=0] terminated
because pmemd complained about the missing
box.
> Anyway, that's my take on this.
>
> HTH,
> Jason
>
Thanks a lot again, TDH (this definitely helps)!
>
> --
> Jason M. Swails
> BioMaPS,
> Rutgers University
> Postdoctoral Researcher
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
--
Dr Massimiliano Porrini
Valérie Gabelica Team
U869 ARNA - Inserm / Bordeaux University
Institut Européen de Chimie et Biologie (IECB)
2, rue Robert Escarpit
33607 Pessac Cedex
FRANCE
Tel : 33 (0)5 40 00 63 31
http://www.iecb.u-bordeaux.fr/teams/GABELICA
Emails: m.porrini.iecb.u-bordeaux.fr
mozz76.gmail.com
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri May 23 2014 - 03:30:02 PDT
