<html>
<head>
<style><!--
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Tahoma
}
--></style>
</head>
<body class='hmmessage'>
Hi,<br><br>The kinetic energy is non-zero because Gromacs uses a leap frog integrator and the kinetic energy<br>at t=0 is the average of the kinetic energies at -dt/2 and +dt/2 (as described in the manual).<br><br>Eric Sorin did comparisons between Amber and Gromacs (uses his Amber force field port in Gromacs)<br>and found differences only of the size of roughly the machine precision. So I don't think there is an issue<br>in Gromacs. But I have no experience with the route you took.<br><br>Berk<br><br>> Date: Tue, 27 Jul 2010 10:43:30 -0400<br>> From: chris.neale@utoronto.ca<br>> To: gmx-users@gromacs.org<br>> Subject: [gmx-users] gromacs energies very different from other MD engine for the very same system and conditions<br>> <br>> Dear Alan:<br>> <br>> I don't have a solution for you, but I do have a simpler test case: a <br>> single spc water molecule. I have not addressed your noted energy <br>> difference between programs because I am not familiar with them <br>> (although you should try gromacs compiled in double precision).<br>> <br>> I can reproduce the non-zero kinetic energy with a single spc water,<br>> as long as I add define=-DFLEXIBLE to the .mdp file. The output .gro<br>> contains all zeroes for velocities (as long as the input .gro does not <br>> contain any velocities!):<br>> <br>> title<br>> 3<br>> 1SOL OW 1 0.230 0.628 0.113 0.0000 0.0000 0.0000<br>> 1SOL HW1 2 0.137 0.626 0.150 0.0000 0.0000 0.0000<br>> 1SOL HW2 3 0.231 0.589 0.021 0.0000 0.0000 0.0000<br>> 10.00000 10.00000 10.00000<br>> <br>> but the run indicates a non-zero kinetic energy (log and .edr file, gromacs<br>> 4.0.5 and 3.3.3).<br>> <br>> The reported kinetic energy is small though, 1.16470e-03, with a<br>> reported temperature of 9.33871e-02.<br>> <br>> I added unconstrainted_start=yes to your .mdp file but that did not help.<br>> <br>> I tried in double precision, and I get the same reported kinetic <br>> energy (1.16473e-03)<br>> <br>> Chris.<br>> <br>> -- original message --<br>> <br>> Hi there,<br>> <br>> I have a very simple case of a tri alanine (AAA). You can use tleap from<br>> ambertools to create such peptide and save the pdb and amber's parameters.<br>> <br>> Then I do a single point calculation to get the potential energy. I am using<br>> mdin (amber input file) like:<br>> <br>> cat << EOF >| mdin<br>> Single point<br>> &cntrl<br>> imin=0, maxcyc=0,<br>> ntmin=2,<br>> ntb=0,<br>> igb=0,<br>> cut=999,<br>> /<br>> EOF<br>> <br>> Then I test the same input parameters with Namd2, using namd.conf (Amber is<br>> not freely available, but AmberTools and Namd are):<br>> <br>> cat << EOF >| AAAamb_namd.conf<br>> outputEnergies 1 # Energy output frequency<br>> DCDfreq 1 # Trajectory file frequency<br>> timestep 2 # in unit of fs<br>> temperature 300 # Initial temp for velocity assignment<br>> cutoff 999<br>> switching off # Turn off the switching functions<br>> PME off # Use PME for electrostatic calculation<br>> amber on # Specify this is AMBER force field<br>> parmfile AAAamb.prmtop # Input PARM file<br>> ambercoor AAAamb.inpcrd # Input coordinate file<br>> outputname AAAamb # Prefix of output files<br>> exclude scaled1-4<br>> 1-4scaling 0.833333 # =1/1.2, default is 1.0<br>> minimize 0<br>> EOF<br>> <br>> I have only a 0.001 kcal/mol absolute diff (or relative 0.0037%) for the<br>> Elec term, everything else is absolutely the same.<br>> <br>> Then comes gromacs. I convert my prmtop and inpcrd to gromacs top and gro<br>> with either acpype or amb2gmx and then doing a single point with file:<br>> <br>> cat << EOF >| SPE.mdp<br>> integrator = md<br>> nsteps = 0<br>> dt = 0.001<br>> constraints = none<br>> emtol = 10.0<br>> emstep = 0.01<br>> nstcomm = 1<br>> ns_type = simple<br>> nstlist = 0<br>> rlist = 0<br>> rcoulomb = 0<br>> rvdw = 0<br>> Tcoupl = no<br>> Pcoupl = no<br>> gen_vel = no<br>> nstxout = 1<br>> pbc = no<br>> nstlog = 1<br>> nstenergy = 1<br>> nstvout = 1<br>> nstfout = 1<br>> nstxtcout = 1<br>> comm_mode = ANGULAR<br>> EOF<br>> <br>> I got something like:<br>> Energies (kJ/mol)<br>> Bond Angle Ryckaert-Bell. LJ-14 Coulomb-14<br>> 1.77662e+02 3.19281e+01 4.55707e+01 2.48926e+01 9.40060e+02<br>> LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy<br>> -9.21735e+00 -1.05005e+03 1.60848e+02 8.84934e+00 1.69697e+02<br>> Temperature Pressure (bar)<br>> 2.28887e+01 0.00000e+00<br>> <br>> First question, why I have temperature and kinetic energy? Remember, it's a<br>> single point energy calculation, without any sort of coupling, 0 velocities<br>> and and gen_vel = no.<br>> <br>> Secondly, converting gromacs energies (kJ/mol) by dividing the terms by<br>> 4.184 and they don't match amber/namd results. See relative diffs<br>> (abs(g-a)/max(abs(a),abs(g))):<br>> <br>> Bond: 5.87%<br>> Angle: 3.49%<br>> Dihe: 0.14%<br>> vdW: 0.52%<br>> Elec: 0.17%<br>> Pot.Tot: 5.90%<br>> <br>> I am very surprised to find this difference, in special for Bond since it's<br>> for this term essentially the same equation and parameters. Any ideas of<br>> what could be missing here?<br>> <br>> Many thanks,<br>> <br>> Alan<br>> <br>> <br>> <br>> -- <br>> gmx-users mailing list gmx-users@gromacs.org<br>> http://lists.gromacs.org/mailman/listinfo/gmx-users<br>> Please search the archive at http://www.gromacs.org/search before posting!<br>> Please don't post (un)subscribe requests to the list. Use the <br>> www interface or send it to gmx-users-request@gromacs.org.<br>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php<br> </body>
</html>