[gmx-users] Is anyone also using lammps?s
Peng Yi
pengyi at MIT.EDU
Thu Oct 29 17:09:56 CET 2009
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
More information about the gromacs.org_gmx-users
mailing list