[gmx-users] Re: Improving Energy Conservation in NVE Simulation of Water

[Justin A. Lemkul]

Andrew DeYoung wrote:
> Hi, 
> 
> I am sorry to bring up a several-days-old post, but if you have time, I am
> wondering if you can please help me with two concepts related to the pdf
> file at http://www.andrew.cmu.edu/user/adeyoung/july27/withnvt.pdf and
> Justin's comments on it.  
> 
> 1.  For temperature coupling in NVT simulations, what is a reasonable value
> for the time constant for coupling (tau_t)?  In my previous simulation in
> which I used tcoupl = v-rescale, I chose tau_t = 2.0 ps.  Justin kindly
> suggested that this may be a bit long for V-rescale.  Would a value of, for
> example, tau_t = 1.0 ps be more reasonable, or should I use something even
> shorter?
> 

The minimum for tau_t is a magnitude 10x your timestep.  Values of tau_t are 
often as low as 0.1 - 0.5 ps.  It is easy to test if this improves your results. 
  Pressure coupling, however, is different and generally these values would be 
too small.  For temperature coupling with weak coupling methods, though, there 
shouldn't be a problem.  Coupling too tightly can result in fixed kinetic energy 
which is unreasonable, so don't set tau_t = dt :)

> 2.  How do you determine, as a rule of thumb, what to use for vdw cutoffs?
> For other simulations of water (i.e., another configuration), I had been
> using cutoffs such as rlist = rvdw = 1.0 nm, and I had been getting good
> results.  However, a few weeks ago, I switched to a new, more interesting
> configuration, and I have been having some trouble with it, as described in
> my previous posts in this thread.  Someone (not on this mailing list)
> suggested that I increase the vdw cutoffs in case the "long-range vdw
> interactions" would make a significant difference.  So, in my last
> simulation, for example, I used rlist = rvdw = 2.3 nm.  I chose this because
> my simulation box is a cube of edgelength 4.71 nm, and the cutoffs must be
> less than half the largest box length in order to obey the minimum image
> convention.  Is it always (or usually) OK to use as large of cutoffs as
> possible without violating the minimum image convention?  Of course, by
> making cutoffs longer than necessary, I am taking a huge hit in terms of
> performance.  But performance aside, is there ever any danger (in terms of
> accuracy) of making the cutoffs too long?  I wouldn't think so, but I am not
> sure.  
> 

Force field models are parameterized to work under a given set of circumstances 
that include cutoffs and other parameters.  Arbitrarily increasing cutoffs can 
absolutely induce artifacts, and I would say it is far more likely that such 
changes cause problems rather than preventing them.  A few relevant posts to 
consider:

http://lists.gromacs.org/pipermail/gmx-users/2009-December/047586.html
http://lists.gromacs.org/pipermail/gmx-users/2009-December/047593.html

Refer to the primary literature for whatever parameters you're using.  Rule of 
thumb: unless you can demonstrate that your conditions better represent reality 
(or are at least similarly accurate), don't mess with them.

> Or, does making the cutoffs too long lead to the "flying ice cube" problem?
> Looking at the article here: http://en.wikipedia.org/wiki/Flying_ice_cube,
> it seems that the flying ice cube problem is more likely related to the
> choice of temperature coupling (it seems my choice of V-rescale was poor; in
> my next simulation perhaps I should try Nose-Hoover) and possibly
> intergration algorithm, rather than related to cutoff lengths.  Do you think
> that this is a correct interpretation?  
> 

My mistake.  In my head I equated the flying ice cube with the artificial 
ordering and density increases that can happen from arbitrary increases to 
cutoff values.  There is a relationship between cutoffs, thermostats, and 
equipartition of energy, but it's not immediately relevant here.

-Justin

> Thank you very much.
> 
> Andrew DeYoung
> Carnegie Mellon University
> 
> Andrew DeYoung wrote:
>> Hi, 
>>
>> Thank you.  I tried again, this time doing 100 ps of NVT equlibration,
>> followed by 500 ps of NVE.  I have posted my results here:
>> http://www.andrew.cmu.edu/user/adeyoung/july27/withnvt.pdf
>>
>> But, they are not much better.  For the NVT equilibration, I used
> v-rescale
>> temperature coupling with tau_t = 2 ps.  What kind of temperature coupling
>> do you typically use or recommend?
>>
> 
> That sounds reasonable, although a 2-ps coupling time is a bit long for 
> V-rescale.  Your results indicate that the system still is not stable, even
> in 
> NVT.  Temperature should not fluctuate very much after the first few ps.
> Try a 
> longer NVT.
> 
> What motivated your choice of cutoffs?  They seem arbitrary and too long.
> You 
> may be seeing effects of the classic "flying ice cube" problem where
> isolated 
> elements of the system start to (essentially) freeze.
> 
> -Justin
> 
>> Thanks,
>>
>> Andrew
>>
>> Andrew DeYoung wrote:
>>> Hi, 
>>>
>>> I have been simulating 1000 SPC/E water molecules in the NVE ensemble (by
>>> using tcoupl = no and pcoupl = no).  In these simulations, I am
>> considering
>>> only Lennard-Jones interactions; I have "turned charges off" by setting
>> the
>>> partial charges on oxygen and hydrogen to zero in a local copy of
> spce.itp
>>> that I include in my topology file.  
>>>
>>> If you have time, I have posted a pdf file containing a summary of what I
>>> have tried so far:
>> http://www.andrew.cmu.edu/user/adeyoung/july27/july27.pdf
>>> At first, I tried using vdwtype = Cut-off.  However, Gromacs gave me a
>> note
>>> message: 
>>> ---
>>> NOTE:
>>> You are using a cut-off for VdW interactions with NVE, for good energy
>>> conservation use vdwtype = Shift (possibly with DispCorr)
>>> ---
>>>
>>> I ignored this note message and ran the simulation anyway.  Indeed, as
>>> Gromacs warned me, the (total) energy conservation is very bad; the total
>>> energy decreases seemingly significantly in what appears to be an almost
>>> linear fashion.  I have done energy minimization, equilibration, and
>>> production dynamics, but for simplicity, I have drawn the equilibration
>> and
>>> production dynamics results in the same color in the pdf file that I have
>>> posted.  
>>>
>>> (As far as energy minimization, I did a series of two energy minimization
>>> steps.  First, I used integrator = l-bfgs with dt = 0.001 ps, nsteps =
>>> 10000, emstep = 0.01 nm, and emtol = 10.0 kJ mol^-1 nm^-1.  I got the
>>> following results:
>>>
>>> Low-Memory BFGS Minimizer converged to Fmax < 10 in 0 steps
>>> Potential Energy  = -9.3022546e+02
>>> Maximum force     =  5.5933684e-01 on atom 1123
>>> Norm of force     =  1.9135195e-01
>>>
>>> So, it seems that the energy of the configuration is quite low--the water
>>> molecules are relatively far apart from one another (the simulation box
> is
>>> quite large).  
>>>
>>> Next, I used integrator = steep (passing the ending configuration from
> the
>>> l-bfgs minimization to grompp using -t followed by the trajectory file)
>> with
>>> dt = 0.001 ps, nsteps = 10000, emstep = 0.01 nm, and emtol = 10.0 kJ
>> mol^-1
>>> nm^-1.  I got the following results:
>>>
>>> Steepest Descents converged to Fmax < 10 in 1 steps
>>> Potential Energy  = -1.2968502e+03
>>> Maximum force     =  7.0817396e-02 on atom 820
>>> Norm of force     =  2.6016856e-02
>>>
>>> This configuration seems probably acceptable.  So much for energy
>>> minimization.)
>>>
>>> Secondly, I used vdwtype = Shift as the note message recommended.  I
> tried
>>> two different sets of values of rlist, rvdw_switch, rvdw, and rcoulomb.
>>> But, in both cases, energy conservation is bad--almost as bad as when I
>> used
>>> vdwtype = Cut-off.  
>>>
>>> If you have time, can you please look at the pdf file that I have posted
>>> (http://www.andrew.cmu.edu/user/adeyoung/july27/july27.pdf)?  Do you have
>>> any suggestions of other vdw types or vdw parameters that I should try to
>>> hopefully improve energy conservation?  On page 1, I have displayed the
>>> total energy results.  On pages 2 and 3, I have pasted a typical .mdp
> file
>>> that I used and also the results of other quantities from the
> simulations.
>>> One thing that I notice is, on the bottom of page 2, it appears that the
>>> _potential_ energy conservation is good, whereas the kinetic energy
>>> conservation is quite bad.  Do you know, what sort of clue should this
>> fact
>>> give me as to what is going on with my system?  
>>>
>> What seems most obvious to me is that your system is cooling down because
>> your 
>> molecules are slowing down.  Your gen_temp is set to 298 K, but your
> initial
>> velocity distribution produces a temperature over 305 K, so you're not
>> exploring 
>> the energy surface you think you are.  Certainly in NVE temperature is not
> 
>> necessarily a conserved quantity (the total energy is), but you're not
>> starting 
>> under the conditions you want.  The temperature in all of these cases is
>> only 
>> approaching the desired 298 K towards the end of the 500 ps timeframe
> you've
>> shown.  Try an NVT equilibration at 298 K first, then switch to NVE.
>>
>> The other thing to consider is whether or not this should even work.
> You've
>> taken a condensed-phase water model that is largely held together by
>> hydrogen 
>> bonding and transformed it into a semi-gas state and removed the most
> potent
>> interactions that keep the molecules together.  Force fields must balance
> LJ
>> and 
>> electrostatic terms to reproduce some desired behavior.  In this case,
>> you're 
>> eliminating half of the equation (or more, depending on how the LJ and 
>> electrostatic terms are balanced to contribute to stability).
>>
>> I would suggest that if you have to hack this model under these
> conditions, 
>> start with NVT equilibration before moving to NVE.  Even better would be
> to
>> use 
>> a simpler system, like an ideal gas that only has LJ terms, and determine
>> the 
>> proper protocol for obtaining a stable NVE ensemble.  That way, you're not
> 
>> stabbing in the dark as to what the appropriate settings should be in the
>> .mdp file.
>>
>> -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

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