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