[MMTK] Does MMTK support negative AMBER PHASE values?
Mark.Williamson at imperial.ac.uk
Thu May 31 13:18:42 UTC 2007
Recently I've been using MMTK's AMBER functionality to model an unusual
organic molecule. The system has a cationic sp2 carbon bound to three
groups, hence some unusual dihedrals. I've managed to get GAFF to work
with MMTK with some minor modifications, but in this process I think I
may have noticed a possible issue with MMTK's handling of negative
dihedral PHASE values.
Therefore, I set up a test to demonstrate this. Here it is:
A model of ethane created and parameterised using GAFF:
All hc: -0.01911
All c3: 0.00637
A minimisation was a carried out using SANDER 9 and the final structure
was converted to pdb using ptraj. A MMTK database file was created for
ethane based the about and the PDB configurations was pointed at the
ptraj generated pdb file. This was then read into MMTK (2.5.19) and an
energy evaluation was carried out. This was compared back with the final
AMBER energy value from the SANDER output. In theory, these two values
should match and indeed they did for TEST #1.
This involved the "default" parameters from gaff.dat. The hc-c3-c3-hc
parameter is highlighted here since it will be varied in test #2:
IPT , JPT , KPT , LPT , IDIVF , PK , PHASE , PN
hc-c3-c3-hc 1 0.15 0.0 3.
Even though this is the default value from gaff.dat, it is loaded in via
a frcmod file to ensure consistency between the two tests.
ALL UNITS IN Kcal / Mol
Lennard Jones = 0.0709 0.0711
Electrostatic = 0.0377 0.0377
Dihedral = 0.0000 0.0000
Bond Length = 0.0013 0.0009
Bond Angle = 0.0129 0.0125
Total Energy = 0.1227 0.1222
Ok, so that's all good and well, there is good agreement.
If I modify the torsion PHASE parameter of hc-c3-c3-hc, so that
it is negative, (I've selected a random non physical value here):
hc-c3-c3-hc 1 0.15 -90.0 3.
regenerate the AMBER topology file and then repeat the above process, I
get the following:
Lennard Jones = 0.1481 0.1470
Electrostatic = 0.0378 0.0378
Dihedral = 2.6970 0.0029 <----
Bond Length = 0.0032 0.0020
Bond Angle = 0.0132 0.0149
Total Energy = 2.8993 0.2045
MMTK is not giving the same energy as SANDER for the dihedral term. More
experimenting seems to indicate a deviation in behaviour away from
SANDER when a dihedral PHASE value is less than 0 degrees or greater
than 180 degrees. Do you know why this is?
I think it has something to do with the following:
seems to have some form of checking for this situation:
if dihedral > Numeric.pi:
dihedral = dihedral - 2.*Numeric.pi
if phase < 0.:
phase = phase + 2.*Numeric.pi
but this function:
does not. Using some crude print()s to debug this, I see that the b)
function is being called. This same b) function is also called when
using a Amber94ForceField forcefield. Should Amber94ForceField be
function a) or should b) have a similar form of checking as a) ?
More information about the mmtk