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).
>
> 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.
> 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.
>
> 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.
Anyway, that's my take on this.
HTH,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu May 22 2014 - 15:00:03 PDT
