[gmx-users] Is anyone also using lammps?s
Ran Friedman
r.friedman at bioc.uzh.ch
Wed Nov 4 10:13:08 CET 2009
Dear Peng,
Did you also try to run GMX in double precision at some point?
Ran
Peng Yi wrote:
>
> I turned off the torsion interaction. The difference between Lammps
> and Gromacs at integration time step 2fs was reduced. Details below:
>
> A melt of 240 n-octane (united-atom model), NVT, T=300K, V=55.46nm^3.
> Both Lammps and Gromacs use berendsen thermostat with tau_t=1ps.
>
> Integration time step 1fs:
> Lammps Gromacs Std. Err. (for both)
> Ebond(kJ/mol): 2092 2109 100
> Eangle: 1778 1754 80
> Elj+corr: -10501 -10553 100
> T(K): 300 299 5
> P(atm): 3188 3016 700
>
> integration time step 2fs:
> Lammps Gromacs Std.Err. (for both)
> Ebond: 2133 2232 100
> Eangle: 1803 1737 80
> Elj+corr: -10501 -10623 100
> T: 300 298 5
> P: 3133 2955 600
>
> Lammps results remain almost unchanged when dt increases from 1fs to 2fs,
> and Ebond : Eangle = 7 : 6, which is the ratio between # of bonds and
> # of angles.
>
> Gromacs results change more significantly when dt goes from 1fs to 2fs.
> and the trend of Ebond and Eangle are opposite. It is more significant
> when torsion interaction is present.
>
>
> On Tue, 3 Nov 2009, David van der Spoel wrote:
>
>> Peng Yi wrote:
>>>
>>> Hi, David,
>>>
>>> I used Berendsen thermostat and a bigger tau_t=1ps to redo the
>>> simulations.
>>> The general conclusion is the same. The std err is the same in both
>>> packages. And during each simulation, the integration time does not
>>> change. Details below:
>>>
>>> Integration time step 1fs:
>>> Lammps Gromacs Std.Err. (for both)
>>> Ebond(kJ/mol): 2112 2170 100
>>> Eangle: 1799 1770 100
>>> Etors: 2552 2490 100
>>> Elj+corr: -10711 -10777 100
>>> P(atm): 3250 3216 500
>>>
>>> Integration time step 2fs:
>>> Lammps Gromacs Std.Err. (for both)
>>> Ebond: 2154 2654 120
>>> Eangle: 1840 1645 120
>>> Etors: 2573 2236 120
>>> Elj+corr: -10711 -11019 100 P(atm):
>>> 3250 2590 600
>>>
>>
>> How about the temperature in both systems? Was Lammps also run with
>> Berendsen? It could also still be a topology error. Maybe you can
>> turn off the torsion potential to test this.
>>
>>> -Peng
>>>
>>> On Sun, 1 Nov 2009, David van der Spoel wrote:
>>>
>>>> Peng Yi wrote:
>>>>>
>>>>> Thank for your reply! I have done some NVT runs per your
>>>>> suggestion, and
>>>>> the results are similar to NPT runs, i.e., Gromacs results is more
>>>>> affected by changing integration timestep than Lammps. Details
>>>>> below:
>>>>>
>>>>> A melt of 240 octane chains by united-atom model. T=300K, V=55.46
>>>>> nm^3.
>>>>> Both Gromacs and Lammps use Nose-Hoover thermostat with tau_t=0.2 ps.
>>>>
>>>> As some people have note the tau_t is short for Nose Hoover. Are
>>>> you sure this means the same in Lammps and in Gromacs? For one
>>>> thing, there is no tau_t in the NH algorithm as far as I know, and
>>>> Gromacs converts it to an appropriate weight or whatever that is
>>>> called. What does Lammps do with this tau_t. To be a the safe side
>>>> you could run both with Berendsen as well. Is the std err identical
>>>> in both packages? And in the 2 fs run, are both simulation
>>>> equlibrated with this time step as well?
>>>>
>>>>> All other parameters in .top and .mdp files are the same as
>>>>> previously
>>>>> attached..
>>>>>
>>>>> If I use integration time 1fs, Lammps and Gromacs produce
>>>>> consistent results:
>>>>> Lammps Gromacs std. err.
>>>>> Ebond(kJ/mol): 2133 2160 100
>>>>> Eangle: 1757 1780 80
>>>>> Etors: 2531 2510 80
>>>>> Elj+corr: -10711 -10767 90
>>>>> P(atm): 3500 3250 500
>>>>>
>>>>> if I use integration time 2fs, Lammps results remain unchanged, but
>>>>> Gromacs results change significantly, particularly bonded energy:
>>>>>
>>>>> Lammps Gromacs std. err.
>>>>> Ebond(kJ/mol): 2175 2710 100
>>>>> Eangle: 1799 1640 70
>>>>> Etors: 2573 2230 80
>>>>> Elj+corr: -10711 -11007 100 P(atm):
>>>>> 3200 2730 700
>>>>>
>>>>> Would that be a result of using different integrator between
>>>>> Lammps and Gromacs? Lammps uses Velocity-Verlet, and Gromacs uses
>>>>> Leap-frog.
>>>>> Thanks,
>>>>> -Peng
>>>>>
>>>>> On Thu, 29 Oct 2009, David van der Spoel wrote:
>>>>>
>>>>>> aherz wrote:
>>>>>>> Hey,
>>>>>>>
>>>>>>> are you running single or double precision gromacs?
>>>>>>> Afaik, depending on the circumstances the energy drift in
>>>>>>> gromacs can be
>>>>>>> rather bad for single precision.
>>>>>>
>>>>>> Please refer to the gromacs 4.0 paper for a discussion of the drift.
>>>>>> If you want to compare energies you need the same density, which
>>>>>> you do not have, you may need to run NVT for that.
>>>>>>
>>>>>> Note that your integration time step is quite large, and the
>>>>>> temperature coupling constant is very small.
>>>>>>
>>>>>> You could try a shifted LJ + dispersion correction, it is not
>>>>>> clear to me how LAMMPS treats cutoffs, couldn't find it in the
>>>>>> manual.
>>>>>>
>>>>>>>
>>>>>>> Alex
>>>>>>>
>>>>>>>
>>>>>>> Peng Yi schrieb:
>>>>>>>> On Wed, 28 Oct 2009, Mark Abraham wrote:
>>>>>>>>
>>>>>>>>> Peng Yi wrote:
>>>>>>>>>> I am trying to simulate alkane melt and found out that
>>>>>>>>>> gromacs and
>>>>>>>>>> lammps gave different results, particularly the bonded
>>>>>>>>>> interaction
>>>>>>>>>> energy.
>>>>>>>>>> I wonder if anyone has such experience. Thanks,
>>>>>>>>> Even two installations of the same version of GROMACS can give
>>>>>>>>> different results. The question is whether when using comparable
>>>>>>>>> model physics you observe the same ensemble averages.
>>>>>>>>>
>>>>>>>>> Mark
>>>>>>>> Hi, Mark,
>>>>>>>>
>>>>>>>> Thanks for reply! The difference is statistically
>>>>>>>> significant. And I am
>>>>>>>> wondering if it is caused by the integrator: Leap-frog for
>>>>>>>> Gromacs and
>>>>>>>> Velocity-verlet for Lammps. Detail description of the
>>>>>>>> comparison please
>>>>>>>> see below:
>>>>>>>>
>>>>>>>> It is an NPT simulation of a melt of 240 n-octane molecules using
>>>>>>>> united-atom model, i.e., CHx group is considered as one atom.
>>>>>>>> There are
>>>>>>>> bond, angle, torsion and LJ interactions. T=300K and P=1atm.
>>>>>>>>
>>>>>>>> Lammps uses nose-hoover thermostat and barostat, and Gromacs uses
>>>>>>>> nose-hoover thermostat and Parranello-Rahman barostat. Time
>>>>>>>> constants
>>>>>>>> for
>>>>>>>> thermostat and barostat are 0.02ps and 2.0ps, respectively.
>>>>>>>>
>>>>>>>> If I use integration time 1fs, Lammps and Gromacs gave consistent
>>>>>>>> results:
>>>>>>>> Lammps Gromacs
>>>>>>>> Ebond(kJ/mol): 2092 2146
>>>>>>>> Eangle: 1757 1760
>>>>>>>> Etors: 2510 2500
>>>>>>>> Elj+corr: -9238 -9350
>>>>>>>> Volume(nm^3): 66.7 66.5
>>>>>>>>
>>>>>>>> where energy fluctuation is 100 kJ/mol and volume fluctuation
>>>>>>>> is 1 nm^3,
>>>>>>>> Elj+corr is the total LJ energy including tail correction.
>>>>>>>>
>>>>>>>> However, if I use integration time 2fs, Lammps results do not
>>>>>>>> change
>>>>>>>> much, but Gromacs results changed a lot:
>>>>>>>>
>>>>>>>> Lammps Gromacs
>>>>>>>> Ebond(kJ/mol): 2133 2700
>>>>>>>> Eangle: 1799 1640
>>>>>>>> Etors: 2552 2200
>>>>>>>> Elj+corr: -9292 -9886
>>>>>>>> Volume: 66.7 64.0
>>>>>>>>
>>>>>>>> The results given by Lammps is more reasonable because the
>>>>>>>> Ebond should
>>>>>>>> be equal to the total # of bonds times 1/2k_BT and Eangle
>>>>>>>> should be equal
>>>>>>>> to the total # of angles times 1/2k_BT. At T=300K,
>>>>>>>> 1/2k_BT=1.25 kJ/mol.
>>>>>>>> 240 n-octanes have total 1680 bonds and 1440 angles.
>>>>>>>>
>>>>>>>> The bond and angle interactions are both harmonic functions. Bond
>>>>>>>> interaction constant kl=292880 kJ/mol/nm^2, corresponding to a
>>>>>>>> bond
>>>>>>>> ossilation period 16 fs.
>>>>>>>>
>>>>>>>> Is there something related to the integrator?
>>>>>>>>
>>>>>>>> Here I attached my grompp.mdp and topol.top files.
>>>>>>>>
>>>>>>>> ##########
>>>>>>>> grompp.mdp
>>>>>>>> ##########
>>>>>>>>
>>>>>>>> ; VARIOUS PREPROCESSING OPTIONS
>>>>>>>> title = Yo
>>>>>>>> cpp = /usr/bin/cpp
>>>>>>>> include = define =
>>>>>>>>
>>>>>>>> ; RUN CONTROL PARAMETERS
>>>>>>>> integrator = md
>>>>>>>> tinit = 0
>>>>>>>> dt = 0.001
>>>>>>>> nsteps = 2000000
>>>>>>>> init_step = 0
>>>>>>>> comm-mode = Linear
>>>>>>>> nstcomm = 1
>>>>>>>> comm-grps =
>>>>>>>>
>>>>>>>> ; OUTPUT CONTROL OPTIONS
>>>>>>>> nstxout = 5000
>>>>>>>> nstvout = 5000
>>>>>>>> nstfout = 5000
>>>>>>>> nstcheckpoint = 10000
>>>>>>>> nstlog = 1000
>>>>>>>> nstenergy = 1000
>>>>>>>> nstxtcout = 5000
>>>>>>>> xtc-precision = 1000
>>>>>>>> xtc-grps = energygrps =
>>>>>>>>
>>>>>>>> ; NEIGHBORSEARCHING PARAMETERS
>>>>>>>> nstlist = 10
>>>>>>>> ns_type = grid
>>>>>>>> pbc = xyz
>>>>>>>> rlist = 1.0025
>>>>>>>> domain-decomposition = no
>>>>>>>>
>>>>>>>> ; OPTIONS FOR ELECTROSTATICS AND VDW
>>>>>>>> coulombtype = Cut-off
>>>>>>>> rcoulomb-switch = 0
>>>>>>>> rcoulomb = 1.0025
>>>>>>>> epsilon-r = 1
>>>>>>>> vdw-type = Cut-off
>>>>>>>> rvdw-switch = 0 ; default rvdw = 1.0025
>>>>>>>> ; default 1 nm
>>>>>>>> DispCorr = EnerPres
>>>>>>>> ;table-extension = 1.5
>>>>>>>> fourierspacing = 0.12
>>>>>>>> fourier_nx = 0
>>>>>>>> fourier_ny = 0
>>>>>>>> fourier_nz = 0
>>>>>>>> pme_order = 4
>>>>>>>> ewald_rtol = 1e-05
>>>>>>>> ewald_geometry = 3d
>>>>>>>> epsilon_surface = 0
>>>>>>>> optimize_fft = no
>>>>>>>>
>>>>>>>>
>>>>>>>> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
>>>>>>>> Tcoupl = nose-hoover
>>>>>>>> tc-grps = System
>>>>>>>> tau_t = 0.02
>>>>>>>> ref_t = 300.0
>>>>>>>> Pcoupl = Parrinello-Rahman
>>>>>>>> Pcoupltype = isotropic
>>>>>>>> tau_p = 2.0
>>>>>>>> compressibility = 4.5e-5
>>>>>>>> ref_p = 1.0
>>>>>>>> andersen_seed = 815131
>>>>>>>>
>>>>>>>> ; GENERATE VELOCITIES FOR STARTUP RUN
>>>>>>>> gen_vel = yes
>>>>>>>> gen_temp = 300
>>>>>>>> gen_seed = 2009
>>>>>>>>
>>>>>>>> ; OPTIONS FOR BONDS constraints = none
>>>>>>>> constraint-algorithm = Lincs
>>>>>>>> unconstrained-start = no
>>>>>>>> Shake-SOR = no
>>>>>>>> shake-tol = 1e-04
>>>>>>>> lincs-order = 4
>>>>>>>> lincs-iter = 1
>>>>>>>> lincs-warnangle = 30
>>>>>>>> morse = no
>>>>>>>>
>>>>>>>> ; ENERGY GROUP EXCLUSIONS
>>>>>>>> ; Pairs of energy groups for which all non-bonded interactions are
>>>>>>>> excluded
>>>>>>>> energygrp_excl =
>>>>>>>>
>>>>>>>> ; NMR refinement stuff disre = No
>>>>>>>> disre-weighting = Conservative
>>>>>>>> disre-mixed = no
>>>>>>>> disre-fc = 1000
>>>>>>>> disre-tau = 0
>>>>>>>> nstdisreout = 100
>>>>>>>> orire = no
>>>>>>>> orire-fc = 0
>>>>>>>> orire-tau = 0
>>>>>>>> orire-fitgrp = nstorireout = 100
>>>>>>>> dihre = No
>>>>>>>> dihre-fc = 1000
>>>>>>>> dihre-tau = 0
>>>>>>>> nstdihreout = 100
>>>>>>>>
>>>>>>>> #########
>>>>>>>> topol.top
>>>>>>>> #########
>>>>>>>>
>>>>>>>> #include "ffG53a6.itp"
>>>>>>>>
>>>>>>>> [atom-types]
>>>>>>>> ;name mass charge ptype V/c6 W/c12
>>>>>>>> CH2 14.0 0.00 A 0.0 0.0
>>>>>>>> CH3 15.0 0.00 A 0.0 0.0
>>>>>>>>
>>>>>>>> [nonbond-params]
>>>>>>>> ; i j func V/c6 W/c12
>>>>>>>> CH2 CH2 1 0.0078 3.24e-5
>>>>>>>> CH2 CH3 1 0.0078 3.24e-5
>>>>>>>> CH3 CH3 1 0.0078 3.24e-5
>>>>>>>>
>>>>>>>> [ moleculetype ]
>>>>>>>> ; name nrexcl
>>>>>>>> Octane1 3
>>>>>>>>
>>>>>>>> [ atoms ]
>>>>>>>> ; nr type resnr residu atom cgnr charge
>>>>>>>> 1 CH3 1 C8 CH3 1 0.0
>>>>>>>> 2 CH2 1 C8 CH2 2 0.0
>>>>>>>> 3 CH2 1 C8 CH2 3 0.0
>>>>>>>> 4 CH2 1 C8 CH2 4 0.0
>>>>>>>> 5 CH2 1 C8 CH2 5 0.0
>>>>>>>> 6 CH2 1 C8 CH2 6 0.0
>>>>>>>> 7 CH2 1 C8 CH2 7 0.0
>>>>>>>> 8 CH3 1 C8 CH3 8 0.0
>>>>>>>>
>>>>>>>> [ bonds ]
>>>>>>>> ; ai aj funct c0(nm) c1(kJ/mol/nm^2)
>>>>>>>> 1 2 1 0.153 292880.0
>>>>>>>> 2 3 1 0.153 292880.0
>>>>>>>> 3 4 1 0.153 292880.0
>>>>>>>> 4 5 1 0.153 292880.0
>>>>>>>> 5 6 1 0.153 292880.0
>>>>>>>> 6 7 1 0.153 292880.0
>>>>>>>> 7 8 1 0.153 292880.0
>>>>>>>>
>>>>>>>> [ pairs ]
>>>>>>>> ; ai aj funct c0 c1
>>>>>>>> ; 1 4 1 0.000000e+00 0.000000e+00 ; 2 5 1
>>>>>>>> 0.000000e+00 0.000000e+00 ; 3 6 1 0.000000e+00
>>>>>>>> 0.000000e+00
>>>>>>>> ; 4 7 1 0.000000e+00 0.000000e+00 ; 5 8 1
>>>>>>>> 0.000000e+00 0.000000e+00
>>>>>>>>
>>>>>>>> [ angles ]
>>>>>>>> ; ai aj ak funct c0(degree) c1(kJ/mol/rad^-2)
>>>>>>>> 1 2 3 1 109.5 502.08
>>>>>>>> 2 3 4 1 109.5 502.08
>>>>>>>> 3 4 5 1 109.5 502.08
>>>>>>>> 4 5 6 1 109.5 502.08
>>>>>>>> 5 6 7 1 109.5 502.08
>>>>>>>> 6 7 8 1 109.5 502.08
>>>>>>>>
>>>>>>>> [ dihedrals ]
>>>>>>>> ; ai aj ak al funct c0 c1 c2
>>>>>>>> c3 c4 c5
>>>>>>>> 1 2 3 4 3 6.4977 16.9868 3.6275
>>>>>>>> -27.112 0 0
>>>>>>>> 2 3 4 5 3 6.4977 16.9868 3.6275
>>>>>>>> -27.112 0 0
>>>>>>>> 3 4 5 6 3 6.4977 16.9868 3.6275
>>>>>>>> -27.112 0 0
>>>>>>>> 4 5 6 7 3 6.4977 16.9868 3.6275
>>>>>>>> -27.112 0 0
>>>>>>>> 5 6 7 8 3 6.4977 16.9868 3.6275
>>>>>>>> -27.112 0 0
>>>>>>>>
>>>>>>>> [ system ]
>>>>>>>> octane melt
>>>>>>>>
>>>>>>>> [ molecules ]
>>>>>>>> Octane1 240
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> gmx-users mailing list gmx-users at gromacs.org
>>>>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>>>>>> Please search the archive at http://www.gromacs.org/search before
>>>>>>>> posting!
>>>>>>>> Please don't post (un)subscribe requests to the list. Use the www
>>>>>>>> interface or send it to gmx-users-request at gromacs.org.
>>>>>>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> gmx-users mailing list gmx-users at gromacs.org
>>>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>>>>> Please search the archive at http://www.gromacs.org/search
>>>>>>> before posting!
>>>>>>> Please don't post (un)subscribe requests to the list. Use the
>>>>>>> www interface or send it to gmx-users-request at gromacs.org.
>>>>>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>>>>>
>>>>>>
>>>>>> --
>>>>>> David.
>>>>>> ________________________________________________________________________
>>>>>> David van der Spoel, PhD, Professor of Biology
>>>>>> 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
>>>>>> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
>>>>>> _______________________________________________
>>>>>> gmx-users mailing list gmx-users at gromacs.org
>>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>>>> Please search the archive at http://www.gromacs.org/search before
>>>>>> posting!
>>>>>> Please don't post (un)subscribe requests to the list. Use the www
>>>>>> interface or send it to gmx-users-request at gromacs.org.
>>>>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>>>>>
>>>>> _______________________________________________
>>>>> gmx-users mailing list gmx-users at gromacs.org
>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>>> Please search the archive at http://www.gromacs.org/search before
>>>>> posting!
>>>>> Please don't post (un)subscribe requests to the list. Use the www
>>>>> interface or send it to gmx-users-request at gromacs.org.
>>>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>>>
>>>> _______________________________________________
>>>> gmx-users mailing list gmx-users at gromacs.org
>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>> Please search the archive at http://www.gromacs.org/search before
>>>> posting!
>>>> Please don't post (un)subscribe requests to the list. Use the www
>>>> interface or send it to gmx-users-request at gromacs.org.
>>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>>>
>>> _______________________________________________
>>> gmx-users mailing list gmx-users at gromacs.org
>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> Please search the archive at http://www.gromacs.org/search before
>>> posting!
>>> Please don't post (un)subscribe requests to the list. Use the www
>>> interface or send it to gmx-users-request at gromacs.org.
>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>
>>
>> --
>> David van der Spoel, Ph.D., Professor of Biology
>> Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University.
>> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. Fax:
>> +4618511755.
>> spoel at xray.bmc.uu.se spoel at gromacs.org http://folding.bmc.uu.se
>> _______________________________________________
>> gmx-users mailing list gmx-users at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> Please search the archive at http://www.gromacs.org/search before
>> posting!
>> Please don't post (un)subscribe requests to the list. Use the www
>> interface or send it to gmx-users-request at gromacs.org.
>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>
> _______________________________________________
> gmx-users mailing list gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before
> posting!
> Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
More information about the gromacs.org_gmx-users
mailing list