<div>Hi there,</div><div><br></div><div>I have a very simple case of a tri alanine (AAA). You can use tleap from ambertools to create such peptide and save the pdb and amber's parameters.</div><div><br></div><div>Then I do a single point calculation to get the potential energy. I am using mdin (amber input file) like:</div>
<div><br></div><div>cat << EOF >| mdin</div><div>Single point</div><div>&cntrl</div><div>imin=0, maxcyc=0,</div><div>ntmin=2,</div><div>ntb=0,</div><div>igb=0,</div><div>cut=999,</div><div>/</div><div>EOF</div>
<div><br></div><div>Then I test the same input parameters with Namd2, using namd.conf (Amber is not freely available, but AmberTools and Namd are):</div><div><br></div><div>cat << EOF >| AAAamb_namd.conf</div><div>
outputEnergies 1 # Energy output frequency</div><div>DCDfreq 1 # Trajectory file frequency</div><div>timestep 2 # in unit of fs</div><div>temperature 300 # Initial temp for velocity assignment</div><div>
cutoff 999</div><div>switching off # Turn off the switching functions</div><div>PME off # Use PME for electrostatic calculation</div><div>amber on # Specify this is AMBER force field</div>
<div>parmfile AAAamb.prmtop # Input PARM file</div><div>ambercoor AAAamb.inpcrd # Input coordinate file</div><div>outputname AAAamb # Prefix of output files</div><div>exclude scaled1-4</div><div>
1-4scaling 0.833333 # =1/1.2, default is 1.0</div>
<div>minimize 0</div><div>EOF</div><div><br></div><div>I have only a 0.001 kcal/mol absolute diff (or relative 0.0037%) for the Elec term, everything else is absolutely the same.</div><div><br></div><div>Then comes gromacs. I convert my prmtop and inpcrd to gromacs top and gro with either acpype or amb2gmx and then doing a single point with file:</div>
<div><br></div><div>cat << EOF >| SPE.mdp</div><div>integrator = md</div><div>nsteps = 0</div><div>dt = 0.001</div><div>constraints = none</div>
<div>
emtol = 10.0</div><div>emstep = 0.01</div><div>nstcomm = 1</div><div>ns_type = simple</div><div>nstlist = 0</div><div>rlist = 0</div>
<div>rcoulomb = 0</div><div>rvdw = 0</div><div>Tcoupl = no</div><div>Pcoupl = no</div><div>gen_vel = no</div><div>nstxout = 1</div>
<div>pbc = no</div><div>nstlog = 1</div><div>nstenergy = 1</div><div>nstvout = 1</div><div>nstfout = 1</div><div>nstxtcout = 1</div><div>comm_mode = ANGULAR</div><div>EOF</div><div><br></div><div>I got something like:</div>
<div><div> Energies (kJ/mol)</div><div> Bond Angle Ryckaert-Bell. LJ-14 Coulomb-14</div><div> 1.77662e+02 3.19281e+01 4.55707e+01 2.48926e+01 9.40060e+02</div><div> LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy</div>
<div> -9.21735e+00 -1.05005e+03 1.60848e+02 8.84934e+00 1.69697e+02</div><div> Temperature Pressure (bar)</div><div> 2.28887e+01 0.00000e+00</div></div><div><br></div><div>First question, why I have temperature and kinetic energy? Remember, it's a single point energy calculation, without any sort of coupling, 0 velocities and and gen_vel = no.</div>
<div><br></div><div>Secondly, converting gromacs energies (kJ/mol) by dividing the terms by 4.184 and they don't match amber/namd results. See relative diffs (abs(g-a)/max(abs(a),abs(g))):</div><div><br></div><div>Bond: 5.87%</div>
<div>Angle: 3.49%</div><div>Dihe: 0.14%</div><div>vdW: 0.52%</div><div>Elec: 0.17%</div><div>Pot.Tot: 5.90%</div><div><br></div><div>I am very surprised to find this difference, in special for Bond since it's for this term essentially the same equation and parameters. Any ideas of what could be missing here?</div>
<div><br></div><div>Many thanks,</div><div><br></div><div>Alan</div><div><br></div><br>-- <br>Alan Wilter S. da Silva, D.Sc. - CCPN Research Associate<br>Department of Biochemistry, University of Cambridge. <br>80 Tennis Court Road, Cambridge CB2 1GA, UK.<br>
>><a href="http://www.bio.cam.ac.uk/~awd28">http://www.bio.cam.ac.uk/~awd28</a><<<br>