Hi All,
This message concerns some LEaP behaviour that I have found when following
up the message by Angela Liu that was posted to the amber mailing list
this
morning. This message was entitled: "AMBER: amber 8 leap topology file
dihedral term problem"
I have confirmed the issue in this message and can confirm that Leap gives
different numbers of dihedral terms for a copy of a unit. My initial
conclusions were that the copy command was broken but having dug further I
have found things get considerably more complicated.
Firstly I have a question concerning how torsion terms involving wild
types
should be interpreted. It would be nice to get a definitive answer on
this.
I.e. In parm99.dat there are the following two terms:
X -CT-CT-X 9 1.400 0.0 3.
OS-CT-CT-H1 1 0.250 0.0 1.
This combination appears for example in a DNA base. Now, in this
situation,
for any dihedrals involving the types OS-CT-CT-H1 what should be written
to
the prmtop file? Should it contain just the 0.250 0.0 1. term or should
it
contain both the 1 fold and the 3 fold term? (0.156 0.0 3.)?
In other words if a torsion angle is explicitly specified in the parm file
does it override ALL wild card torsion terms that also match or does it
only
override ones that have the same periodicity?
I have always believed it was the former while Dave Case thinks it should
be
the latter? Can one of the force field developers please let me know which
is correct.
For the record Leap typically (BUT NOT ALWAYS) seems to employ the former
case, I.e. the 3 fold term is NOT used for the OS-CT-CT-H1 system.
Although
this seems strange since I would think that one wants to preserve the 3
fold
rotational barrier for this system... Comments?
Secondly I just want to ask what the true benefit of wildcards is? Do we
really need them? Although it would result in a larger parmxx.dat file
surely it would make things a LOT less ambiguous if we specified every
dihedral explicitly. The reason I ask this is that leap does not seems to
implement the wild type torsion parameters in a consistent fashion. They
also cause problems with the copy command in Leap.
E.g. for a single base of DNA (nuc.pdb attached) if you load that in Leap
and do the following:
source leaprc.ff99
dna = loadpdb nuc.pdb
dna2 = copy dna
saveamberparm dna dna.prmtop dna.inpcrd
saveamberparm dna2 dna2.prmtop dna2.inpcrd
You will find that dna2.prmtop has 3 extra dihedral terms. These
correspond
to 3 fold X-CT-CT-X being added for O4'-C4'-C5'-H5'1, O4'-C4'-C5'-H5'2 and
O4'-C4'-C3'-H3' in addition to the explicitly specified 1 fold OS-CT-CT-H1
term. However, other examples where this X-CT-CT-X term should be added if
things were consistent would be H4'-C4'-C3'-O3' (there are many other
examples as well). In both the original AND the copy only the 1 fold term
(0.250 0.00 1.0) is present. Strangely if you do:
dna3 = copy dna2
saveamberparm dna3 dna3.prmtop dna3.inpcrd
you get back to dna.prmtop.
So in this case, the file dna.prmtop is correct assuming that ANY
explicitly
specified torsion parameter should overide the wild card parameters even
if
the phase does not match. While dna2.prmtop is definately wrong as it only
has SOME of the wild card matches added.
However, the story does not end here. If you load nma.pdb (attached) and
do:
nma = loadpdb nma.pdb
nma2 = copy nma
saveamberparm nma nma.prmtop nma.inpcrd
saveamberparm nma2 nma2.prmtop nma2.inpcrd
You will find that in this case the nma.prmtop has MORE dihedral terms
than
the nma2.prmtop. Although this time the two prmtop files give the same
energies since the extra dihedral terms have zero barrier height. However,
this is fortuitous in my opinion and indicative of a deeper problem. If
you
print out the dihedrals with rdparm you'll see that the first prmtop has
the
following two extra terms:
10: 0.000 0.00 2.00 :1.HH32 :1.CH3 :1.C :1.O (3,2,5,6)
15: 0.000 0.00 2.00 :1.HH31 :1.CH3 :1.C :1.O (1,2,5,6)
Two things are very strange here. Firstly there is are two explicit
HC-CT-C-O dihedrals:
HC-CT-C-O 1 0.80 0.0 -1.
HC-CT-C-O 1 0.08 180.0 3.
Both of which are included. Now, if leap was being consistent in it's
treatment of wild type parameters as in the DNA case above then that
should
be all there is in the prmtop file. In the one created from the COPIED
unit
this is true but in the uncopied one, the one created by loading the pdb
file leap has included the X-C -CT- X 2 fold term. BUT it has only
included
it for 2 out of the 3 hydrogens!!! Why??? Fortunately here it doesn't make
any difference since the barrier height is zero but if you were to make
the
X-C -CT- X term non-zero you would have a serious problem.
These are just two simple examples, I'm sure there are many more and more
serious examples.
I think we should check very carefully that the parameters to which the
force fields are being fit are the same as those generated for the prmtop
file. Should we perhaps do away with the wild cards all together to avoid
these ambiguities???
Your comments....
All the best
Ross
/\
\/
|\oss Walker
| Department of Molecular Biology TPC15 |
| The Scripps Research Institute |
| Tel:- +1 858 784 8889 | EMail:- ross.rosswalker.co.uk |
|
http://www.rosswalker.co.uk/ | PGP Key available on request |
Note: Electronic Mail is not secure, has no guarantee of delivery, may not
be read every day, and should not be used for urgent or sensitive issues.
- application/octet-stream attachment: nuc.pdb
- application/octet-stream attachment: nma.pdb
Received on Wed Apr 05 2006 - 23:49:56 PDT
