[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