[gmx-users] Is anyone also using lammps?s
Peng Yi
pengyi at MIT.EDU
Tue Nov 3 16:12:26 CET 2009
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
-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
>
More information about the gromacs.org_gmx-users
mailing list