[gmx-users] Is anyone also using lammps?s
Ran Friedman
r.friedman at bioc.uzh.ch
Thu Oct 29 18:33:38 CET 2009
Hi Peng,
The time scale should be much shorter than the fastest vibration. A rule
of thumb from the reference below is a factor of ten, but it would
depend on the precision. Running with double precision is shorter but I
didn't make benchmarks (perhaps other users have).
Appropriate values of tau_t and tau_p have also been discussed in this
list (search for references by Berk). I tend to use something like
tau_t=0.2 and tau_p=1.0.
I advise you to follow David's suggestion as well and run with NVT.
Ran.
Reference:
@book{Becker2001,
Author = {Becker, O. M. and MacKerell, A. D. Jr. and Roux, B. and
Watanabe, M.},
Title = {Computational biochemistry and biophysics},
Publisher = {Dekker, M.},
Address = {New York},
Year = {2001} }
Peng Yi wrote:
>
> 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
>>
> _______________________________________________
> 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