Hi!
I'm building saccharide monomer 0BC from Glycam 06h and found
that tleap does not put atoms of improper dihedral in correct
order. Moreover it isn't even consistent within one molecule
where atoms of the same improper type get ordered differently.
According to AmberTools and PREP file documentation atoms of
impropers are ordered as follows:
- Central atom is always in 3rd position.
- The three atoms bonded to the central atom are ordered by
their atom type.
- In case of identical atom types, they are sorted by atom
index.
But tleap produces this:
B 65: 1.500 3.14 2.0 :2.C2N :2.H2N :2.N2 :2.C2 (30,29,28,5)
B 66: 1.500 3.14 2.0 :2.C4N :2.C4 :2.N4 :2.H4N (20,9,18,19)
With Glycam 06g it put both of them in the same order:
B 65: 1.500 3.14 2.0 :2.C2 :2.C2N :2.N2 :2.H2N (26,30,28,29)
B 66: 1.500 3.14 2.0 :2.C4 :2.C4N :2.N4 :2.H4N (12,16,14,15)
The difference is that Glycam 06g had this improper specified
using wildcard atoms (X X N H) while in version 06h it is given
explicitly (C H Ng Cg).
So tleap actually has another rule for atom ordering in impropers:
- Atoms matched by wildcards are moved at the beginning of the
list and sorted by index.
Is it intentional and as intended? If so documentation should be
amended.
I think that atoms should be ordered by their type and index as
tie-breaker regardless if it happened to be matched by improper
with wildcards in it or not. That would be consistent with
documentation, but tleap should be corrected.
In any case, there is an inconsistency somewhere which results
in different order of the same improper type in one molecule.
The correct order would be this:
B 65: 1.500 3.14 2.0 :2.C2N :2.C2 :2.N2 :2.H2N (30,5,28,29)
B 66: 1.500 3.14 2.0 :2.C4N :2.C4 :2.N4 :2.H4N (20,9,18,19)
Or if ordered as in parameter file:
B 65: 1.500 3.14 2.0 :2.C2N :2.H2N :2.N2 :2.C2 (30,29,28,5)
B 66: 1.500 3.14 2.0 :2.C4N :2.H4N :2.N4 :2.C4 (20,19,18,9)
The cause of different order within the same molecule in case of
Glycam 06h probably is a problem in
leap/parmSet.c:ParmSetImproperOrderAtoms() or in the loop in
leap/unitio.c:zbUnitIOIndexTorsionParameters() from which it is
called.
Here is the explanation as I understand it:
- Improper parameters in leap are stored in TORSIONPARMt struct.
- This struct has 'sOrder' field.
- sOrder is basically a map between atom types and atom indexes
in the improper or in case of the parameters it maps sorted
atom types to the order they were in the original parameter
file.
- When improper parameters are read, their atom types are sorted
according to the rules above and sOrder is changed to reflect
how the order changed.
- When improper parameters are added to specific torsion,
atom types of the relevant atoms are sorted according to the
rules above and sOrder is changed. Since only the order of
atom types was changed, but not the order of the atoms, sOrder
now tells which atom types correspond to which atoms. sOrder
is stored in the torsion with the new atom type order and old
atom order.
- After that ParmSetImproperOrderAtoms() is called and should
change atom order to correspond to the order of atom types as
shown by sOrder.
- In this function the atoms which are given explicitly (without
wildcards) are ordered first without looking at sOrder at all.
- Then both, atoms and their types, are reordered according to
sOrder.
- After that atoms matching wildcards are moved to the beginning
and sorted by their index.
- And finally atoms with identical atom types are sorted by
their index.
- After this in unitio.c:zbUnitIOIndexTorsionParameters() the
function iParmSetAddImproperTerm() is called which resets
sOrder to '0123' since now it receives already sorted atom
types. In the next loop atom order will not be changed to the
one in original parameter file.
So the issue is that sOrder is changed in between the two calls
to ParmSetImproperOrderAtoms() for the same improper. For some
reason iParmSetAddImproperTerm() is called again (the first time
was when loading parameter file) and now it finds sOrder='0123'
when in the first case it was '0321'. So atoms of the second
improper are not reordered leading to different atom orders for
the two impropers.
I'm not familiar with inner workings of leap so I don't know
what would be the correct fix for this.
I have not checked if there is similar issue to proper dihedrals
and other terms.
Any comments and suggestions would be appreciated. Also,
pointing out if I got something wrong would be helpful.
Best Regards,
Reinis Danne
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue May 07 2013 - 06:30:02 PDT
