[gmx-users] [Fwd: Energy conservation]

[Janne Hirvi]

> > I would be really pleased if you could help me with energy conservation
> > problem I have struggle for so long time. First of all I would like to
> > thank you for sending me your parameter file related to your interesting
> > article "The Origin of Layer Structure Artifacts in Simulations of Liquid
> > Water". Those parameters where pretty equal to mine and I couldn't achieve
> > good energy conservation and would like to know if I am doing something
> > really stupid?
> > 
> > -------------------------------------------
> > ;    File 'mdout.mdp' was generated
> > ;    By user: spoel (291)
> > ;    On host: matisse
> > ;    At date: Sun Aug 11 22:20:40 2002
> > 
> > dt                       = 0.002
> > nstlist                  = 5
> > ns-type                  = Grid
> > pbc                      = xyz
> > rlist                    = 0.9
> > coulombtype              = PME
> > rcoulomb                 = 0.9
> > vdw-type                 = Cut-off
> > rvdw                     = 0.9
> > fourierspacing           = 0.12
> > pme_order                = 4
> > ewald_rtol               = 1e-05
> > tcoupl                   = Berendsen
> > tau-t                    = 0.1
> > ref-t                    = 298.15
> > Pcoupl                   = Berendsen
> > Pcoupltype               = Isotropic
> > tau-p                    = 1.0
> > compressibility          = 5e-5
> > ref-p                    = 1
> > -------------------------------------------
> > 
> > If I use previous parameters (without pressure and temperature coupling) to
> > simulate bulk water with 1372 SPC/E water molecules and with double
> > precision (first equilibriated) I observe drift of approximately +800kJ/mol
> > in the total energy during 200ps. If I try to remove the effect of van der
> > Waals cut-off using switch function energy drift increases from before to
> > +1900kJ/mol/200ps, which was quite a surprise. Still enhancing the
> > parameters by decreasing the value of ewald_rtol to 1E-6 finally decreases
> > energy drift for +200kJ/mol/200ps. Good energy conservation seems to need
> > even much smaller value of ewald_rtol and thus larger value of beta.
> > 
> > ------------------------------------------------------
> > Cut-off(vdW)(0.9)           1E-5   +793kJ/mol/200ps
> > Switch(vdW) (0.8/0.85/0.9)  1E-5  +1878kJ/mol/200ps
> > Switch(vdW) (0.8/0.85/0.9)  1E-6   +236kJ/mol/200ps
> > ------------------------------------------------------
> 
> Did you use switch or shift? I would definitely recommend the shifted 
> potential (for Van der Waals) over a slightly larger range, at least 2 
> angstrom
>
I used switch function and tested also larger range (0.7/0.8/0.9), which
actually increased the energy drift to +2100kJ/mol/200ps. I also just tested
that the same result was achieved with shift function (0.7/0.8/0.9). However
this interesting behaviour isnt that actual concern because I managed to
achieve pretty good energy conservation with much smaller value of ewald_rtol
atleast when using parameters (rlist=rcoulomb=rvdw=1.2, nstlist=5). So I can
now achieve good energy conservation (if necessary) but I am worried about my
earlier simulations with non energy conservative parameters.
>
> Did you use SETTLE for the constraints?
> 
Yes, I used SETTLE as 'a standard option' from the SPC/E model file.
>
> > 
> > Previous values are almost equal or maybe little bit larger than energy 
> > drifts when using parameters I have employed in my studies 
> > (rlist=rcoulomb=rvdw=1.2, nstlist=10). Use of van der Waals cut-off always 
> > decreases the positive energy drift or changes its sign to negative.
> > 
> > Now I think that I understand much more about energy conservation when 
> > using PME, but I am worried about the my earlier simulations with non
> > energy conservative parameters and have some questions:
> > 
> > 1)First of all I would like to know if you agree with this large energy 
> > drift with these parameters or have I done something stupid?
> First off all you did not divide by the number of molecules, so your 
> drift really is 0.1 to  0.5 kJ/mol
>
Thats true, I understand that, but however drift of 0.2-1.3kJ/mol (total
~45kJ/mol) seemed really large to me atleast when I want to simulate longer
times than 200ps.
> 
> > 
> > 2)If you agree with these values do you see any problems simulating with 
> > non energy conservative parameters when temperature coupling is used?
> Barring other problems, I think that is fine.
> 
I am relieved if you agree that temperature coupling can be used to mask not so
perfect energy conservation (even as large drifts as mentioned earlier).
> 
> > 
> > 3)Do you think/agree that when van der Waals cut-off is employed there 
> > is no need to be so demanding about the PME parameters?
> It depends on the application. For normal simulation with T,P coupling 
> the default parameters are fine. For e.g. free energy calculations more 
> accuracy may be needed.
> 
I have made bulk simulations and wetting simulations (water droplet spreading on
a surface) and I have to admit that I dont know what makes free energy
calculations so special? I meant that even with perfect PME parameters we can't
achieve energy conservation because of van der Waals cut-off so why to use the
computationally heaviest parameters for PME? 
> 
> > 
> > I will be very delighted if you have time to comment these issues because I
> > think you are one of the real experts on this area and I would really need
> > experts clarification to my own thoughts.
> > 
> > Thanks,
> > 
> > Janne
> > 
> >
>
------------------------------------------------------------------------------
> > Janne Hirvi, MSc(Physical Chemistry), Researcher
> > University of Joensuu, Department of Chemistry, P.O.Box 111 80101 
> > Joensuu, FI
> > Tel: +358 13 2514544 & +358 50 3474223
> > E-mail: Janne.Hirvi at joensuu.fi & hirvi at cc.joensuu.fi
>
------------------------------------------------------------------------------
> 
> David.
> ________________________________________________________________________
> David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
> Dept. of Cell and Molecular Biology, Uppsala University.
> Husargatan 3, Box 596,  	75124 Uppsala, Sweden
> phone:	46 18 471 4205		fax: 46 18 511 755
> spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++