[MMTK] Difference in number of Torsion parameters between MMTK and Leap

Mark Williamson Mark.Williamson at imperial.ac.uk
Wed May 2 14:15:26 UTC 2007


Dear All,

I've recently been comparing the energy outputs from MMTK and sander for
identical systems and noticed a difference with the dihedrals terms. Let
me initially explain the protocol of the test that I'm carrying out.

I'm using AMBER-9 with all patches up to level 17 applied and MMTK
2.5.19. The test system is the PDB bala1 distributed with MMTK.


On the AMBER side of things, I used the following leap script to
construct a topology and initial crd:

[mjw99 at mjw leap]$ more leap.bat
verbosity 10

# Load the structure
bala = loadPdb bALA1.pdb

#Save it all
saveAmberParm bala bala.prmtop bala.crd

quit


This is then executed thus:

$AMBERHOME/exe/tleap -f leaprc.ff94 -f leap.bat

Next, I then run some minimization on this output using these directives:

[mjw99 at mjw md]$ more stage1.in
Stage 1 minimisation
 &cntrl
  imin   = 1,
  maxcyc = 10000,
  ncyc   = 500,
  ntb    = 0,
  ntr    = 0,
  cut    = 50
 &end

and execute using sander:

        $AMBERHOME/exe/sander -O \
                -i stage1.in \
                -o stage1.log \
                -p bala.prmtop \
                -c bala.crd \
                -r stage1.crd


Finally, I convert the final structure to a PDB using ptraj:

$AMBERHOME/exe/ptraj bala.prmtop << EO
        trajin stage1.crd
        trajout stage1.pdb pdb
        go
EO
mv stage1.pdb.1 stage1.pdb


The final energy output from sander log file is as follows:

                    FINAL RESULTS

NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
9499      -2.8276E+01     9.1502E-05     2.7444E-04     CH3        19

 BOND    =   0.5081  ANGLE   =     1.0879  DIHED      =        4.5292
 VDWAALS =  -0.5892  EEL     =   -82.5769  HBOND      =        0.0000
 1-4 VDW =   2.7354  1-4 EEL =    46.0298  RESTRAINT  =        0.0000

and it should correspond to that of the configuration in the PDB.





On the MMTK side of things, I used script to load the MINIMISED pdb that
I've just generated above to evaluate its energy and display how many
different terms are present in the force field:



[mjw99 at mjw mmtk]$ more go_small.py
#!/usr/bin/python

from MMTK import *
from MMTK.Proteins import Protein
from MMTK.ForceFields import Amber94ForceField
from MMTK.PDB import PDBConfiguration
import sys

universe = InfiniteUniverse(forcefield=Amber94ForceField())
conf = PDBConfiguration('stage1.pdb')
chain = conf.peptide_chains[0].createPeptideChain(c_terminus=1)
universe.addObject(Protein([chain]))

amber_energy_unit = Units.kcal/Units.mol

print 'Units in Kcal / Mol'
print 'Lennard Jones = ',universe.energyTerms()['Lennard-Jones']  /
amber_energy_unit
print 'Electrostatic = ',universe.energyTerms()['electrostatic'] /
amber_energy_unit
print 'Dihedral Angles = ',universe.energyTerms()['cosine dihedral
angle'] / amber_energy_unit
print 'Bond Length = ',universe.energyTerms()['harmonic bond']  /
amber_energy_unit
print 'Bond Angle = ',universe.energyTerms()['harmonic bond angle'] /
amber_energy_unit
print 'Total Energy = ',universe.energy() / amber_energy_unit

parameters = universe.energyEvaluatorParameters()

for p_type in parameters:
 print p_type, len(parameters[p_type])



and it's output:

[mjw99 at mjw mmtk]$ ./go_small.py
Units in Kcal / Mol
Lennard Jones =  2.14783482849
Electrostatic =  -36.5517186704
Dihedral Angles =  4.51897906175
Bond Length =  0.506530763856
Bond Angle =  1.09535608693
Total Energy =  -28.2830179294

nonbonded 3
lennard_jones 4
electrostatic 4
harmonic_distance_term 21
harmonic_angle_term 36
cosine_dihedral_term 32




Ok, the energy results are in good agreement with AMBER FF94 including
the dihedral term, but if I repeat the above this on a bigger system
such as crambin, all energy terms EXCEPT the dihedral term are in
general agreement between the two MD codes.
Keeping with this bala1 test case here, I used the rdparm util to look
at the topology present in the prmtop file:

/opt/amber-9.17/exe/rdparm bala.prmtop

Using the commands angle and bond I get the same number of terms as in
the second part of the MMTK output (21 and 36), but the dihedrals output
gives a different number 52 vs 32. This is the output from rdparm:
(justified a little for readability)


RDPARM MENU:  dihedrals
Mask [] represents 22 atoms

Dihedral    pk     phase pn                atoms
  1:   2.000  0.00  1.0    :1 at O    :1 at C    :2 at N    :2 at H    (6,5,7,8)
E 2:   2.500  3.14  2.0    :1 at O    :1 at C    :2 at N    :2 at H    (6,5,7,8)
  3:   0.000  0.00  2.0    :1 at C    :2 at N    :2 at CA   :2 at HA   (5,7,9,10)
  4:   0.000  0.00  2.0    :1 at HH33 :1 at CH3  :1 at C    :1 at O    (4,2,5,6)
  5:   0.000  0.00  2.0    :1 at HH33 :1 at CH3  :1 at C    :2 at N    (4,2,5,7)
  6:   0.000  0.00  2.0    :1 at HH32 :1 at CH3  :1 at C    :1 at O    (3,2,5,6)
  7:   0.000  0.00  2.0    :1 at HH32 :1 at CH3  :1 at C    :2 at N    (3,2,5,7)
  8:   2.500  3.14  2.0    :1 at CH3  :1 at C    :2 at N    :2 at H    (2,5,7,8)
  9:   0.000  0.00  2.0    :1 at HH31 :1 at CH3  :1 at C    :1 at O    (1,2,5,6)
 10:   0.000  0.00  2.0    :1 at HH31 :1 at CH3  :1 at C    :2 at N    (1,2,5,7)
 11:   2.000  0.00  1.0    :2 at O    :2 at C    :3 at N    :3 at H    (16,15,17,18)
E12:   2.500  3.14  2.0    :2 at O    :2 at C    :3 at N    :3 at H    (16,15,17,18)
 13:   0.000  0.00  2.0    :2 at C    :3 at N    :3 at CH3  :3 at HH31 (15,17,19,20)
 14:   0.000  0.00  2.0    :2 at C    :3 at N    :3 at CH3  :3 at HH32 (15,17,19,21)
 15:   0.000  0.00  2.0    :2 at C    :3 at N    :3 at CH3  :3 at HH33 (15,17,19,22)
 16:   0.156  0.00  3.0    :2 at HB3  :2 at CB   :2 at CA   :2 at C    (14,11,9,15)
 17:   0.156  0.00  3.0    :2 at HB2  :2 at CB   :2 at CA   :2 at C    (13,11,9,15)
 18:   0.156  0.00  3.0    :2 at HB1  :2 at CB   :2 at CA   :2 at C    (12,11,9,15)
 19:   0.156  0.00  3.0    :2 at HA   :2 at CA   :2 at CB   :2 at HB1  (10,9,11,12)
 20:   0.156  0.00  3.0    :2 at HA   :2 at CA   :2 at CB   :2 at HB2  (10,9,11,13)
 21:   0.156  0.00  3.0    :2 at HA   :2 at CA   :2 at CB   :2 at HB3  (10,9,11,14)
 22:   0.000  0.00  2.0    :2 at HA   :2 at CA   :2 at C    :2 at O    (10,9,15,16)
 23:   0.000  0.00  2.0    :2 at HA   :2 at CA   :2 at C    :3 at N    (10,9,15,17)
 24:   2.500  3.14  2.0    :2 at CA   :2 at C    :3 at N    :3 at H    (9,15,17,18)
 25:   0.000  0.00  2.0    :2 at H    :2 at N    :2 at CA   :2 at HA   (8,7,9,10)
 26:   0.000  0.00  2.0    :2 at H    :2 at N    :2 at CA   :2 at CB   (8,7,9,11)
 27:   0.000  0.00  2.0    :2 at H    :2 at N    :2 at CA   :2 at C    (8,7,9,15)
 28:   0.156  0.00  3.0    :2 at N    :2 at CA   :2 at CB   :2 at HB1  (7,9,11,12)
 29:   0.156  0.00  3.0    :2 at N    :2 at CA   :2 at CB   :2 at HB2  (7,9,11,13)
 30:   0.156  0.00  3.0    :2 at N    :2 at CA   :2 at CB   :2 at HB3  (7,9,11,14)
 31:   0.000  0.00  2.0    :3 at H    :3 at N    :3 at CH3  :3 at HH31 (18,17,19,20)
 32:   0.000  0.00  2.0    :3 at H    :3 at N    :3 at CH3  :3 at HH32 (18,17,19,21)
 33:   0.000  0.00  2.0    :3 at H    :3 at N    :3 at CH3  :3 at HH33 (18,17,19,22)
B34:   1.000  3.14  2.0    :1 at C    :2 at CA   :2 at N    :2 at H    (5,9,7,8)
B35:   1.000  3.14  2.0    :2 at C    :3 at CH3  :3 at N    :3 at H    (15,19,17,18)

 36:   2.500  3.14  2.0    :1 at O    :1 at C    :2 at N    :2 at CA   (6,5,7,9)
 37:   0.530  0.00  1.0    :1 at C    :2 at N    :2 at CA   :2 at CB   (5,7,9,11)
E38:   0.150  3.14  3.0    :1 at C    :2 at N    :2 at CA   :2 at CB   (5,7,9,11)
E39:   0.500  3.14  4.0    :1 at C    :2 at N    :2 at CA   :2 at CB   (5,7,9,11)
 40:   0.200  3.14  2.0    :1 at C    :2 at N    :2 at CA   :2 at C    (5,7,9,15)
 41:   2.500  3.14  2.0    :1 at CH3  :1 at C    :2 at N    :2 at CA   (2,5,7,9)
 42:   2.500  3.14  2.0    :2 at O    :2 at C    :3 at N    :3 at CH3  (16,15,17,19)
 43:   0.000  0.00  2.0    :2 at CB   :2 at CA   :2 at C    :2 at O    (11,9,15,16)
 44:   0.070  0.00  2.0    :2 at CB   :2 at CA   :2 at C    :3 at N    (11,9,15,17)
E45:   0.100  0.00  4.0    :2 at CB   :2 at CA   :2 at C    :3 at N    (11,9,15,17)
 46:   2.500  3.14  2.0    :2 at CA   :2 at C    :3 at N    :3 at CH3  (9,15,17,19)
 47:   0.000  0.00  2.0    :2 at N    :2 at CA   :2 at C    :2 at O    (7,9,15,16)
 48:   0.750  3.14  1.0    :2 at N    :2 at CA   :2 at C    :3 at N    (7,9,15,17)
E49:   1.350  3.14  2.0    :2 at N    :2 at CA   :2 at C    :3 at N    (7,9,15,17)
E50:   0.400  3.14  4.0    :2 at N    :2 at CA   :2 at C    :3 at N    (7,9,15,17)
B51:  10.500  3.14  2.0    :1 at CH3  :2 at N    :1 at C    :1 at O    (2,7,5,6)
B52:  10.500  3.14  2.0    :2 at CA   :3 at N    :2 at C    :2 at O    (9,17,15,16)

One can hopefully see why I decided to use this smaller test case over
crambin :)
I think (it does not say this in the AMBER manual) that the dihedrals
prefixed with E are additional terms in Fourier series to the dihedral
on the line before. I also think that the terms with B are impropers.

Therefore, there are 52 - 7(E) = 45 dihedrals here, 4 of which are
improper. This still does not match MMTK's number of dihedrals? This
difference in the number of dihedrals is also seen with crambin.

Can this be explained? Is this real or a problem in my comparison protocol?

I am hoping to use MMTK as a tool for some AMBER parameter development
and I need to be sure that I am comparing like with like before I
progress onto the next stage.

regards,

Mark





More information about the mmtk mailing list