[gmx-users] Is anyone also using lammps?s
David van der Spoel
spoel at xray.bmc.uu.se
Thu Oct 29 17:21:01 CET 2009
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
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
More information about the gromacs.org_gmx-users
mailing list