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.
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