[gmx-users] Is anyone also using lammps?s
David van der Spoel
spoel at xray.bmc.uu.se
Sun Nov 1 19:19:45 CET 2009
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
More information about the gromacs.org_gmx-users
mailing list