[gmx-users] Is anyone also using lammps?s

Peng Yi pengyi at MIT.EDU
Sun Nov 1 17:53:18 CET 2009


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.
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
>



More information about the gromacs.org_gmx-users mailing list