[gmx-users] Recommended protocol for cooling down temperature in 50 K steps

[Justin A. Lemkul]

Andrew DeYoung wrote:
> Hi, 
> 
> I have a system that I would like to equilibrate in the NPT ensemble at high
> temperature, and then, in subsequent equilibration steps, cool down in steps
> of 50 K.  Thus, in the first step, I will use "gen_vel = yes", while in all
> subsequent steps, I will use "gen_vel = no".
> 
> Suppose that I want to run 200 ps at 700 K, 200 ps at 650 K, and 200 ps at
> 600 K.  I have two possible protocols, Protocol 1 and Protocol 2.  If you
> have time, can you please help me see what is the difference between these
> two protocols?  For my application (cooling down a system in steps of 50 K),
> which protocol do you think is more appropriate?  Or do you have any other
> suggestions?  
> 
> In Protocol 1, both the output configuration (.gro) and the checkpoint
> (.cpt) resulting from mdrun of the previous step are fed to grompp of the
> next step (however, I am not appending files, because I want separate .trr
> files for each step):  
> 
> ***BEGIN PROTOCOL 1***
> STEP 1 (700 K, 200 ps):
> grompp -f a1.mdp -c inputconf.gro -p topology.top -o a1.tpr -po a1out.mdp
> ;where a1.mdp contains "gen_vel = yes"
> mdrun -deffnm a1 -cpo a1.cpt -noappend
> 
> STEP 2 (650 K, 200 ps):
> grompp -f a2.mdp -c a1.gro -p topology.top -t a1.cpt -o a2.tpr -po a2out.mdp
> ;where a2.mdp contains "gen_vel = no"
> mdrun -deffnm a2 -cpo a2.cpt -noappend
> 
> STEP 3 (600 K, 200 ps):
> grompp -f a3.mdp -c a2.gro -p topology.top -t a2.cpt -o a3.tpr -po a3out.mdp
> ;where a3.mdp contains "gen_vel = no"
> mdrun -deffnm a3 -cpo a3.cpt -noappend
> ***END PROTOCOL 1***
> 
> In Protocol 2, only the output configuration (.gro) resulting from mdrun of
> the previous step are fed to grompp of the next step:
> 
> ***BEGIN PROTOCOL 2***
> STEP 1 (700 K, 200 ps):
> grompp -f a1.mdp -c inputconf.gro -p topology.top -o a1.tpr -po a1out.mdp
> ;where a1.mdp contains "gen_vel = yes"
> mdrun -deffnm a1
> 
> STEP 2 (650 K, 200 ps):
> grompp -f a2.mdp -c a1.gro -p topology.top -o a2.tpr ;where a2.mdp contains
> "gen_vel = no"
> mdrun -deffnm a2
> 
> STEP 3 (600 K, 200 ps):
> grompp -f a3.mdp -c a2.gro -p topology.top -o a3.tpr ;where a3.mdp contains
> "gen_vel = no"
> mdrun -deffnm a3
> ***END PROTOCOL 2***
> 
> Which is preferable?  In both protocols, the end velocities of the previous
> step are used to begin the next step, since the .gro files contain the
> velocities.  However, only in Protocol 1 do I use the checkpoint files to
> continue the run.  Does doing so just cause the energy statistics to be more
> consistent between steps?
> 

I doubt either protocol really does the best job.  Your procedure does not 
produce a true "cooling" but instead does three separate simulations, and in the 
final two steps you're initializing the simulation with the wrong ensemble and 
trying to let the system recover.

The best approach is simulated annealing.  You can either do a linear decrease 
in temperature over time, or decrease for a period, then hold the temperature 
constant, make another change, and repeat as necessary.  The details are 
described in the manual.

> Also, if you have time, should I use "continuation = yes" in a2.mdp and
> a3.mdp?  My system does not have any constraints per se; atoms are either
> frozen (using a freeze group to simulate an immobilized surface) or their
> bonds are flexible.  Thus, I am not really using any constraints, unless
> freeze groups count as constraints.  However, I am not sure what this means
> as far as what to set "continuation" to in a2.mdp and a3.mdp.
> 

The "continuation" keyword applies to bond constraints and should not affect 
frozen groups.

-Justin

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

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

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