[gmx-users] Inconsistence of SPC/E and TIP4P water in Gromacs

Thu May 11 17:30:38 CEST 2017

On 5/11/17 1:53 AM, Mark Abraham wrote:
> Hi,
>
> On Thu, May 11, 2017 at 3:31 AM ZUO Taisen <zuots at ihep.ac.cn> wrote:
>
>>
>> Hi guys:
>>
>>      I have compared the SPC/E and TIP4P water of the opls force field in
>> Gromacs.But there is something inconsistent of the potential energy in the
>> system
>>
>>
>> The potential energy of SPC/E water is -4.71400e+04/1000=-47.14kJ/mol
>> which is far from the results from literature (Peter G. Kusalik,
>> Science,vol,265,26,p1219-1221,1994)
>> The potential energy TIP4P water is -4.17875e+04/1000=41.787kJ/mol which
>> is very close to the literature -41.8kJ/mol(Peter G. Kusalik,
>> Science,vol,265,26,p1219-1221,1994)
>>
>>
>> I have attached all the files simulating the SPC/E and TIP4P water(.trr
>> .xtc files are delected) and also the related paper in my last email but
>> was suspended because of the files were too big.So no files was attached
>> this time but the .mdp file was posted as following.
>>
>
> We shouldn't have a mailing list where thousands of recipients can receive
> large emails :-)
>
>
>> I'm eager to know why! Thank you very much!
>>
>
> Force fields, including water models, are parameterized to reproduce
> certain theoretical or experimental data. To the extent that those targets
> resemble other observables, you might expect agreement. You should expect
> the level of such agreement to differ between models.
>

Especially when the methods applied are different.  The original Berendsen 1987 
SPC/E paper used very different methods relative to what the OP shows in the 
.mdp file.  In the original paper, the cutoffs were 0.9 nm, no PME, Berendsen's 
weak coupling method for temperature and pressure, temperature at 300 K, SHAKE 
instead of SETTLE (doubt that will make much of a difference, honestly), and the 
simulations were only a whopping 27.5 ps.  So there's no doubt in my mind that 
one will get different results.  Also, every paper that involves SPC/E in some 
form of comparison reports a different diffusion coefficient (though this is not 
unique to SPC/E), and most are not corrected for system size effects. 
Reproducibility is a considerable challenge...

-Justin

> Mark
>
>
>> TIP4P simulation results:
>> Statistics over 30000001 steps using 300001 frames
>>
>>    Energies (kJ/mol)
>>         LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.      Potential
>>           7.71821e+03   -2.05973e+02   -4.93956e+04    9.58136e+01
>>  -4.17875e+04
>>     Kinetic En.   Total Energy    Temperature Pres. DC (bar) Pressure (bar)
>>     7.30494e+03   -3.44826e+04    2.93006e+02   -1.14009e+02    1.41117e+00
>>
>>
>> SPC/E simulation results:
>>
>> Statistics over 30000001 steps using 300001 frames
>>
>>
>>    Energies (kJ/mol)
>>         LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.      Potential
>>     9.23352e+03   -2.12206e+02   -5.62576e+04    9.63633e+01   -4.71400e+04
>>     Kinetic En.   Total Energy    Temperature Pres. DC (bar) Pressure (bar)
>>     7.30494e+03   -3.98350e+04    2.93005e+02   -1.18040e+02    1.64654e+00
>>
>>
>>
>>
>> title= OPLS MD waterethanol siimulation
>> ; Run parameters
>> integrator= md; leap-frog integrator
>> nsteps= 30000000; 1 * 20000000 = 20000 ps (20 ns)
>> dt    = 0.001; 1 fs
>> ; Output control
>> nstxout        = 100000; save coordinates every 10.0 ps
>> nstvout        = 100000; save velocities every 10.0 ps
>> nstenergy        = 100000; save energies every 10.0 ps
>> nstlog        = 100000; update log file every 10.0 ps
>> nstxout-compressed  = 100000      ; save compressed coordinates every 10.0
>> ps
>> define                  =-DEFLEXIBLE                         ; flexible
>> water
>> compressed-x-grps   = System    ; replaces xtc-grps
>> ; Bond parameters
>>  continuation        = yes; Restarting after NPT
>>  constraint_algorithm    = lincs    ; holonomic constraints
>>  constraints            = none; H bonds (even heavy atom-H bonds)
>> constrained
>>  lincs_iter            = 1    ; accuracy of LINCS
>>  lincs_order            = 4    ; also related to accuracy
>> ; Neighborsearching
>> cutoff-scheme       = Verlet
>> ns_type    = grid; search neighboring grid cells
>> nstlist    = 10    ; 20 fs, largely irrelevant with Verlet scheme
>> rcoulomb    = 1.2; short-range electrostatic cutoff (in nm)
>> rvdw    = 1.2; short-range van der Waals cutoff (in nm)
>> ; Electrostatics
>> coulombtype    = PME; Particle Mesh Ewald for long-range electrostatics
>> pme_order    = 4    ; cubic interpolation
>> fourierspacing= 0.12; grid spacing for FFT
>> ; Temperature coupling is on
>> tcoupl= V-rescale            ; modified Berendsen thermostat
>> tc-grps= SOL          ; two coupling groups - more accurate
>> tau_t= 0.1        ; time constant, in ps
>> ref_t= 293          ; reference temperature, one for each group, in K
>> ; Pressure coupling is on
>> pcoupl        = Parrinello-Rahman    ; Pressure coupling on in NPT
>> pcoupltype        = isotropic            ; uniform scaling of box vectors
>> tau_p        = 2.0            ; time constant, in ps
>> ref_p        = 1.0            ; reference pressure, in bar
>> compressibility     = 4.5e-5            ; isothermal compressibility of
>> water, bar^-1
>> ; Periodic boundary conditions
>> pbc= xyz; 3-D PBC
>> ; Dispersion correction
>> DispCorr= EnerPres; account for cut-off vdW scheme
>> ; Velocity generation
>> gen_vel= no; Velocity generation is off
>>
>>
>>
>> --
>>
>>
>> Taisen Zuo
>>
>> China Spallation Neutron Source，Institute of High Energy Physics， Chinese
>> Academy of Science
>> A1-510, Zhongziyuan road NO.1, Dongguan, Guangdong, PR China. 523770
>> Tel: 86-0769-89156495
>> Cell: 13650469795
>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-request at gromacs.org.

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

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

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