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