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 120Elj+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, andthe 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 80Elj+corr: -10711 -11007 100 P(atm): 3200 2730 700Would 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 berather 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. MarkHi, 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 pleasesee below: It is an NPT simulation of a melt of 240 n-octane molecules usingunited-atom model, i.e., CHx group is considered as one atom. There arebond, angle, torsion and LJ interactions. T=300K and P=1atm. Lammps uses nose-hoover thermostat and barostat, and Gromacs usesnose-hoover thermostat and Parranello-Rahman barostat. Time constantsfor 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.5where 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 GromacsEbond(kJ/mol): 2133 2700 Eangle: 1799 1640Etors: 2552 2200Elj+corr: -9292 -9886 Volume: 66.7 64.0The 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-offrvdw-switch = 0 ; default rvdw = 1.0025 ; default 1 nmDispCorr = 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 10.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-usersPlease 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-usersPlease 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-usersPlease 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-usersPlease 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
