[gmx-users] Parrinello-Rahman and Simulated Annealing

[Justin Lemkul]

On 5/21/13 7:55 AM, Baptiste Demoulin wrote:
> Dear Gromacs users,
>
> I am trying to simulate a large systems embeded in a membrane (roughly
> 360,000 atoms). I would like to equilibrate it at 300K and 1 atm. To this
> mean, I use simulated annealing to slowly warm up the system from 50 to
> 300K, under NPT conditions.
>
> For the simulation, I use position restraints on the heavy atoms on the
> protein, and stronger constraints on the metallic core, and LINCS to
> constrain all the bonds, and 2 fs timestep.
>
> My problem is that when I use Parrinello-Rahman barostat for the
> simulation, I got LINCS warning concerning the restrained metallic core.
> However, using Berendsen barostat for the same system, same input, seems to
> work fine, (at least for the 200,900 steps I did until now). I thus tried
> various modifications, such as changing the coupling tau_p, running first a
> pure NVT at 50K for 2 ps, and then warming, but that did not work... Is
> there an incompatibility between Parrinello-Rahman and simulated annealing
> ?
>

Parrinello-Rahman allows for much wider fluctuations than the Berendsen 
algorithm, so I can see how warming the system may cause problems unless it is 
very well equilibrated beforehand.  Assuming that the annealing step is a 
preparatory step, there is no compelling need to use P-R, is there?  If there 
is, you may just have to do the annealing over a much longer time period such 
that the change in temperature is (very) long with respect to tau_p, to allow 
for the barostat to adjust to changing conditions.

-Justin

> Here is my input:
>
> define               = -DPOSRES -DPOSRES_CORE
> integrator           = md
> continuation         = no           ; first dynamics run
> nsteps               = 1500000
> dt                   = 0.002
>
> nstcomm              = 1
> comm-mode            = linear
> comm-grps            = Protein_POPC SOL_NA_CL
>
> pbc                  = xyz
> nstlist              = 10
> ns-type              = grid
> rlist                = 1.0
> rlistlong            = 1.6 ; the force on a particle which falls within
> rlist and rlistlong is calculated during the pair list update and retained
> during nslist steps
> coulombtype          = pme
> rcoulomb             = 1.0
> vdwtype              = cut-off
> rvdw                 = 1.4
> DispCorr             = EnerPres
>
> nstxtcout            = 10
> nstxout              = 0
> nstvout              = 0
> nstfout              = 0
> nstlog               = 10
>
> constraints          = all-bonds
> constraint-algorithm = lincs
>
> tcoupl               = V-rescale
> tc-grps              = Protein POPC SOL_NA_CL
> tau_t                = 0.1  0.1  0.1
> ref_t                = 300  300  300
>
> pcoupl               = Parrinello-Rahman
> pcoupltype           = semiisotropic
> compressibility      = 4.5e-5 4.5e-5
> tau_p                = 10.0
> ref_p                = 1.0 1.0
>
> refcoord_scaling     = com
>
> gen_vel              = no
> gen_temp             = 50
> gen_seed             = -1
>
> annealing            = single single single
> annealing_npoints    = 2 2 2
> annealing_time       = 0 2000 0 2000 0 2000
> annealing_temp       = 50 300 50 300 50 300
>
> Thank you very much for your help !
>
> Baptiste
>

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================