[gmx-users] Is anyone also using lammps?s
Ran Friedman
r.friedman at bioc.uzh.ch
Thu Oct 29 17:21:38 CET 2009
Hi Peng,
Note that you're not using any bond constraints in Gromacs and a
timestep of 2fs may be too long.
Also, tau_t=0.02 seems too short for me.
With 1fs timescale the agreement seem good enough, but you didn't
include estimated errors so it's hard to tell. Also, I assume you run
GMX in single and LAMMPS in double precision. Did you check for convergence?
Ran
Peng Yi wrote:
>
> 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
More information about the gromacs.org_gmx-users
mailing list