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

[Mark Abraham]

On 9/03/2012 1:08 PM, Justin A. Lemkul wrote:
>
>
> 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.

Agreed. Also, even though .gro files can save velocities, they are 
insufficiently accurate for preserving a simulation ensemble. See 
http://www.gromacs.org/Documentation/How-tos/Doing_Restarts. It's the 
least significant issue with your workflow, however.

Mark