Dear Peng, Did you also try to run GMX in double precision at some point?
Ran Peng Yi wrote: > > I turned off the torsion interaction. The difference between Lammps > and Gromacs at integration time step 2fs was reduced. Details below: > > A melt of 240 n-octane (united-atom model), NVT, T=300K, V=55.46nm^3. > Both Lammps and Gromacs use berendsen thermostat with tau_t=1ps. > > Integration time step 1fs: > Lammps Gromacs Std. Err. (for both) > Ebond(kJ/mol): 2092 2109 100 > Eangle: 1778 1754 80 > Elj+corr: -10501 -10553 100 > T(K): 300 299 5 > P(atm): 3188 3016 700 > > integration time step 2fs: > Lammps Gromacs Std.Err. (for both) > Ebond: 2133 2232 100 > Eangle: 1803 1737 80 > Elj+corr: -10501 -10623 100 > T: 300 298 5 > P: 3133 2955 600 > > Lammps results remain almost unchanged when dt increases from 1fs to 2fs, > and Ebond : Eangle = 7 : 6, which is the ratio between # of bonds and > # of angles. > > Gromacs results change more significantly when dt goes from 1fs to 2fs. > and the trend of Ebond and Eangle are opposite. It is more significant > when torsion interaction is present. > > > On Tue, 3 Nov 2009, David van der Spoel wrote: > >> Peng Yi wrote: >>> >>> Hi, David, >>> >>> I used Berendsen thermostat and a bigger tau_t=1ps to redo the >>> simulations. >>> The general conclusion is the same. The std err is the same in both >>> packages. And during each simulation, the integration time does not >>> change. Details below: >>> >>> Integration time step 1fs: >>> Lammps Gromacs Std.Err. (for both) >>> Ebond(kJ/mol): 2112 2170 100 >>> Eangle: 1799 1770 100 >>> Etors: 2552 2490 100 >>> Elj+corr: -10711 -10777 100 >>> P(atm): 3250 3216 500 >>> >>> Integration time step 2fs: >>> Lammps Gromacs Std.Err. (for both) >>> Ebond: 2154 2654 120 >>> Eangle: 1840 1645 120 >>> Etors: 2573 2236 120 >>> Elj+corr: -10711 -11019 100 P(atm): >>> 3250 2590 600 >>> >> >> How about the temperature in both systems? Was Lammps also run with >> Berendsen? It could also still be a topology error. Maybe you can >> turn off the torsion potential to test this. >> >>> -Peng >>> >>> On Sun, 1 Nov 2009, David van der Spoel wrote: >>> >>>> Peng Yi wrote: >>>>> >>>>> Thank for your reply! I have done some NVT runs per your >>>>> suggestion, and >>>>> the results are similar to NPT runs, i.e., Gromacs results is more >>>>> affected by changing integration timestep than Lammps. Details >>>>> below: >>>>> >>>>> A melt of 240 octane chains by united-atom model. T=300K, V=55.46 >>>>> nm^3. >>>>> Both Gromacs and Lammps use Nose-Hoover thermostat with tau_t=0.2 ps. >>>> >>>> As some people have note the tau_t is short for Nose Hoover. Are >>>> you sure this means the same in Lammps and in Gromacs? For one >>>> thing, there is no tau_t in the NH algorithm as far as I know, and >>>> Gromacs converts it to an appropriate weight or whatever that is >>>> called. What does Lammps do with this tau_t. To be a the safe side >>>> you could run both with Berendsen as well. Is the std err identical >>>> in both packages? And in the 2 fs run, are both simulation >>>> equlibrated with this time step as well? >>>> >>>>> All other parameters in .top and .mdp files are the same as >>>>> previously >>>>> attached.. >>>>> >>>>> If I use integration time 1fs, Lammps and Gromacs produce >>>>> consistent results: >>>>> Lammps Gromacs std. err. >>>>> Ebond(kJ/mol): 2133 2160 100 >>>>> Eangle: 1757 1780 80 >>>>> Etors: 2531 2510 80 >>>>> Elj+corr: -10711 -10767 90 >>>>> P(atm): 3500 3250 500 >>>>> >>>>> if I use integration time 2fs, Lammps results remain unchanged, but >>>>> Gromacs results change significantly, particularly bonded energy: >>>>> >>>>> Lammps Gromacs std. err. >>>>> Ebond(kJ/mol): 2175 2710 100 >>>>> Eangle: 1799 1640 70 >>>>> Etors: 2573 2230 80 >>>>> Elj+corr: -10711 -11007 100 P(atm): >>>>> 3200 2730 700 >>>>> >>>>> Would that be a result of using different integrator between >>>>> Lammps and Gromacs? Lammps uses Velocity-Verlet, and Gromacs uses >>>>> Leap-frog. >>>>> Thanks, >>>>> -Peng >>>>> >>>>> On Thu, 29 Oct 2009, David van der Spoel wrote: >>>>> >>>>>> 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@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 >>>>>> >>>>>> >>>>>> -- >>>>>> 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 >>>>>> sp...@xray.bmc.uu.se sp...@gromacs.org http://folding.bmc.uu.se >>>>>> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ >>>>>> _______________________________________________ >>>>>> 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 >>>> >>>> _______________________________________________ >>>> 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 >> >> >> -- >> David van der Spoel, Ph.D., Professor of Biology >> Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University. >> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. Fax: >> +4618511755. >> sp...@xray.bmc.uu.se sp...@gromacs.org http://folding.bmc.uu.se >> _______________________________________________ >> 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 > _______________________________________________ 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