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@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-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
_______________________________________________
gmx-users mailing list gmx-users@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-requ...@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
sp...@xray.bmc.uu.se sp...@gromacs.org http://folding.bmc.uu.se
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
_______________________________________________
gmx-users mailing list gmx-users@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-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
_______________________________________________
gmx-users mailing list gmx-users@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-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
