[gmx-users] Is anyone also using lammps?s
Peng Yi
pengyi at MIT.EDU
Thu Oct 29 18:24:54 CET 2009
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.
>
I tried NPT simulation with only LJ interaction. With tau_t=0.2ps for
both Gromacs and Lammps, all the other parameters same as I mentioned
before. At integration time 2fs, Gromacs and Lammps produce same results
and agree with published data. (N=4000, T=2.0, P=0.5 -> rou=0.3). That's
why I focused my attention to the difference on bonded energy.
I think I used a simple cutoff plus tail correction for energy and
pressure for both Gromacs and Lammps.
>>
>> 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