[gmx-users] Gromacs and Lammps: energy conservation very different
James
jamesresearching at gmail.com
Tue May 6 13:20:12 CEST 2014
Dear Michael and Mark,
As an update, it seems there was a units error in my Lammps input file. Now
both Lammps and Gromacs have catastrophic energy drift.
Apologies for taking your time with this unconventional "validation" of
Gromacs against Lammps. Thank you again for the discussion.
With best regards,
James
On 6 May 2014 18:42, James <jamesresearching at gmail.com> wrote:
> Dear Michael and Mark,
>
> Thanks a lot for your replies.
> The final 100ps of the 500ps equilibration for Gromacs and Lammps yield
> the following average energies (all kJ/mol):
>
> Total energy per molecule: Gromacs = -37.1, Lammps = -36.8
> PE per molecule: Gromacs = -44.5, Lammps = -44.3
> KE per molecule: Gromacs = 7.5, Lammps = 7.5
>
> So they match remarkably well (within 1 standard deviation of the
> fluctuations).
>
> The drift rate for Lammps is 20 kJ/mol/ns/molecule. That for Gromacs is
> exponential, taking about 50ps to raise the initial 20 kJ/mol/molecule!
>
> I have tried with the default LINCS and also shake (the latter with the
> same tolerance as Lammps), but this made no difference, unfortunately.
>
> I will ask the Lammps mailing list to see if anything goes on behind the
> scenes, although my initial intuition is that 20 kJ/mol/ns/molecule is
> fairly reasonable (ie, it's poor, but it's not exponential). I will let you
> know the response.
>
> I'm interested in this, not because I want to simulate water with a simple
> cut-off per se, but because I'm doing other calculations utilising
> short-term auto-correlations which seem to work with Lammps but give an
> over-estimate in Gromacs, and I'm wondering if it's connected to this. At
> least if I can rule out a trivial error in my input files, there can be
> confidence to look deeper.
>
> Thanks again for your help.
> With best regards,
>
> James
>
>
>
>
> On 6 May 2014 15:09, Mark Abraham <mark.j.abraham at gmail.com> wrote:
>
>> On May 6, 2014 4:16 AM, "Michael Shirts" <mrshirts at gmail.com> wrote:
>> >
>> > What are the energies with the same initial configuration? Perhaps
>> LAMMPS
>> > is doing something under the hood like shifting the cutoff which Gromacs
>> > won't do unless you ask. Knowing the step=0 energy will help identify if
>> it
>> > is a difference in dynamics or in the energy function.
>>
>> Using the potential-shift modifier is indeed a good idea, but it should
>> only affect the noisiness, not the drift.
>>
>> > Simulating water with a cutoff is generally a bad idea (LJ = maybe OK,
>> > water = bad). It is interesting that lammps fails less badly, but this
>> is
>> > guaranteed to fail in either program. Lammps working almost suggests
>> that
>> > it's doing something behind the scenes.
>>
>> Also, using the LINCS settings intended for conservation might help there.
>>
>> The OP reporting the magnitude of the drift would be nice ;-) Any drift
>> looks horrific on the wrong scale!
>>
>> Mark
>>
>> > On Mon, May 5, 2014 at 8:21 PM, James <jamesresearching at gmail.com>
>> wrote:
>> >
>> > > Dear all,
>> > >
>> > > (resending due to hitting 50kb limit on mailing messages)
>> > >
>> > > Considering a simple SPC/E water simulation, I am getting catastrohpic
>> > > energy drift in Gromacs but only minor drift in Lammps for the same
>> system
>> > > (see attached graph).
>> > >
>> > > I guess I'm making an embarrassingly trivial mistake somehow with
>> Gromacs,
>> > > but if anyone might be able to point out where I'm going wrong it
>> would
>> be
>> > > very much appreciated!
>> > >
>> > > Details are below, and the md.log output can be found here:
>> > > https://gist.github.com/jamesresearching/908e4b8804c96a49c61c
>> > >
>> > > The systems start from the same configuration, are equilibrated for
>> 500ps
>> > > with NVT, and then run under NVE. Both use a simple cut-off. I
>> understand
>> > > this will not conserve energy perfectly, but it shouldn't be as bad as
>> > > Gromacs suggests.
>> > >
>> > > Both systems
>> > > - Update neighbour list at every step
>> > > - Apply COM velocity correction every 10 steps
>> > > - Same temperature
>> > > - Same LJ combining rule
>> > > - Same cut-off
>> > > - Same time-step
>> > >
>> > > Running Gromacs in double-precision doesn't help.
>> > > Switching from LINCS to Shake (constraint-algorithm = SHAKE,
>> shake-tol =
>> > > 0.000001) in Gromacs doesn't help either (Lammps uses shake with the
>> same
>> > > tolerance). Example files below are with LINCS (default).
>> > >
>> > > The equilibration files are as follows:
>> > >
>> > > ___top.top___
>> > > #include "FF.itp"
>> > > #include "mol.itp"
>> > >
>> > > [ system ]
>> > > ; Name
>> > > SPCEWater
>> > >
>> > > [ molecules ]
>> > > ; Compound #mols
>> > > SPCE 901
>> > > ------------------------------
>> > >
>> > >
>> > > ___FF.itp___
>> > > [ defaults ]
>> > > ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
>> > > 1 2 no 0.5 0.5
>> > >
>> > > [ atomtypes ]
>> > > ; name bond_type mass charge ptype sigma
>> epsilon
>> > > OW OW 15.9994 -0.8476 A 0.3166 0.65
>> > > HW HW 1.008 0.4238 A 0.0 0.0
>> > > ------------------------------
>> > >
>> > >
>> > > ___mol.itp___
>> > > ; molname nrexcl
>> > > SPCE 2
>> > >
>> > > [ atoms ]
>> > > ; id at type res nr residu name at name cg nr charge mass
>> > > 1 OW 1 SOL OW 1 -0.8476 15.9994
>> > > 2 HW 1 SOL HW 2 0.4238 1.008
>> > > 3 HW 1 SOL HW 3 0.4238 1.008
>> > >
>> > >
>> > > [ settles ]
>> > > ; OW funct doh dhh
>> > > 1 1 0.1000 0.1632991
>> > >
>> > > [ exclusions ]
>> > > 1 2 3
>> > > 2 1 3
>> > > 3 1 2
>> > > --------------------------------
>> > >
>> > >
>> > > ___parameters file (Equilibration)___
>> > > cutoff-scheme = group
>> > > constraints = all-angles
>> > > integrator = md-vv
>> > > dt = 0.002
>> > > nsteps = 250000
>> > > nstlist = 1 ; neighbour-list update period
>> > > ns_type = grid ; restrict building of new neighbour lists
>> to
>> > > those atoms in neighbouring cells
>> > > coulombtype = cut-off
>> > >
>> > > rlist = 1.2
>> > > rcoulomb = 1.2 ; Maximum distance for short-range coulomb
>> > > calculation
>> > > rvdw = 1.2 ; Maximum distance for short-range vdw
>> > > calculation
>> > >
>> > >
>> > > nstxout = 0 ; period with which to write out co-ordinates
>> to a
>> > > trajectory file
>> > > nstvout = 0 ; " velocities
>> > > nstfout = 0 ; " forces
>> > > nstlog = 500 ; period with which to write out energies to
>> log
>> > > file
>> > > nstenergy = 500 ; period with which to write out energies to
>> > > energy file
>> > > energygrps = System ; the groups to write to the energy file
>> > > nstxtcout = 500 ; period with which to write out positions
>> to
>> an
>> > > xtc file
>> > > xtc-precision = 10000 ; xtc precision
>> > >
>> > > nstcomm = 10 ; number of steps between COM removal
>> > > comm_mode = Linear ; remove translational COM (other options
>> > > include angular)
>> > > comm_grps = System ; remove COM for the whole system
>> > >
>> > > Tcoupl = nose-hoover ; Thermostatting
>> > > tc-grps = System ; can also couple groups separately to
>> bath
>> > > tau_t = 0.5 ; nose-hoover time-constant (ps)
>> > > ref_t = 300.0 ; Temperature
>> > >
>> > > DispCorr = EnerPres
>> > >
>> > > Pcoupl = no
>> > >
>> > > gen_vel = yes
>> > > gen_temp = 300.0
>> > > gen_seed = 173529
>> > > -------------------------------------
>> > >
>> > >
>> > > The production calculation is then performed as:
>> > >
>> > > ___parameters file (production)___
>> > > cutoff-scheme = group
>> > > constraints = all-angles
>> > > integrator = md-vv
>> > > dt = 0.002
>> > > nsteps = 250000
>> > > nstlist = 1 ; neighbour-list update period
>> > > ns_type = grid ; restrict building of new neighbour lists
>> to
>> > > those atoms in neighbouring cells
>> > > coulombtype = cut-off
>> > >
>> > > rlist = 1.2
>> > > rcoulomb = 1.2 ; Maximum distance for short-range coulomb
>> > > calculation
>> > > rvdw = 1.2 ; Maximum distance for short-range vdw
>> > > calculation
>> > >
>> > >
>> > > nstxout = 0 ; period with which to write out co-ordinates
>> to a
>> > > trajectory file
>> > > nstvout = 0; " velocities
>> > > nstfout = 0 ; " forces
>> > > nstlog = 100 ; period with which to write out energies to
>> log
>> > > file
>> > > nstenergy = 100 ; period with which to write out energies to
>> > > energy file
>> > > energygrps = System ; the groups to write to the energy file
>> > > nstxtcout = 0 ; period with which to write out positions to
>> an
>> > > xtc file
>> > > xtc-precision = 10000 ; xtc precision
>> > >
>> > > nstcomm = 10 ; number of steps between COM removal
>> > > comm_mode = Linear ; remove translational COM (other options
>> > > include angular)
>> > > comm_grps = System ; remove COM for the whole system
>> > >
>> > > Tcoupl = no
>> > > DispCorr = EnerPres
>> > > Pcoupl = no
>> > > -------------------------------
>> > >
>> > > Production is after the equilibration simulation and run with grompp
>> before
>> > > execution:
>> > > grompp -f parameters.mdp -c ../Equil/confout.gro -p ../Equil/top.top
>> -t
>> > > ../Equil/state.cpt
>> > >
>> > > For reference, the lammps simulation runs with the following input
>> file:
>> > >
>> > > ___Lammps input file___
>> > > units real
>> > > atom_style full
>> > >
>> > > read_data lammps.data # Same position configuration as Gromacs
>> > >
>> > > pair_style lj/cut/coul/cut 12
>> > >
>> > > pair_coeff 1 1 0.1554 3.166 # kcal/mol = 0.65 kJ/mol
>> > > pair_coeff * 2 0.0000 0.0000
>> > >
>> > > bond_style harmonic
>> > > angle_style harmonic
>> > >
>> > > bond_coeff 1 1000.00 1.0
>> > > angle_coeff 1 100.0 109.47
>> > >
>> > > neighbor 0.0 bin
>> > > neigh_modify every 1
>> > >
>> > > timestep 0.002
>> > >
>> > > fix 1 all shake 0.000001 20 0 b 1 a 1
>> > > fix NVT all nvt temp 300.0 300.0 0.5
>> > >
>> > > velocity all create 300 432567 dist uniform
>> > > fix linmom all momentum 10 linear 1 1 1
>> > >
>> > > thermo_style custom step temp etotal ke pe
>> > > variable mystep equal step*0.002
>> > > variable myetotal equal etotal
>> > > variable mype equal pe
>> > > variable myke equal ke
>> > > variable myevdwl equal evdwl
>> > > variable myecoul equal ecoul
>> > > variable mytemp equal temp
>> > >
>> > > fix printenergy all print 100 "${mystep} ${myetotal} ${mype} ${myke}
>> > > ${myevdwl} ${myecoul}" file sim.energy.equil
>> > > fix printtemp all print 100 "${mystep} ${mytemp}" file sim.temp.equil
>> > >
>> > > # NVT equilibration
>> > > run 250000
>> > >
>> > > # NVE production
>> > > unfix NVT
>> > > fix NVE all nve
>> > > reset_timestep 0
>> > >
>> > > fix printenergy all print 100 "${mystep} ${myetotal} ${mype} ${myke}
>> > > ${myevdwl} ${myecoul}" file sim.energy.prod
>> > > fix printtemp all print 100 "${mystep} ${mytemp}" file sim.temp.prod
>> > >
>> > > run 250000
>> > > -------------------------------
>> > >
>> > > Best regards,
>> > >
>> > > James
>> > >
>> > > --
>> > > 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.
>> > >
>> > >
>> > --
>> > 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.
>> --
>> 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.
>>
>
>
More information about the gromacs.org_gmx-users
mailing list