[gmx-users] Is anyone also using lammps?s
Peng Yi
pengyi at MIT.EDU
Thu Oct 29 18:10:20 CET 2009
Hi, Ran,
I didn't use bond restraints. I checked that the bond length had a
Gaussian-like distributes, and the length range looked normal.
I estimated the fastest timescale in the system, which is the bond
ossilation period, around 16fs. Would that require as integration
timestep much smaller than 1fs?
With the parameters I have, could you recommend a set of tau_t and tau_p?
I did mention the fluctuation, 100 kJ/mol for energy and 1 nm^3 for
volume. And I ran GMX in single. Not sure about Lammps, should be
double. All measured physical quantities converged well. Would you
expect differece if I compile GMX in double? Would that be much slower?
-Peng
On Thu, 29 Oct 2009, Ran Friedman wrote:
> 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
>
> _______________________________________________
> 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