Improper ordering in MDAnalysis

[lilyminium]

Expected behavior

Impropers, where the central atom is typically located in the third index located in varying non-symmetric places (see below), are faithfully represented in MDAnalysis.

Actual behavior

MDAnalysis shows a representation where Improper.values[0] < Improper.values[-1].

>>> top.atoms[[31,34,43,35]].impropers.indices
array([[31, 34, 43, 35],  # <--- this one
       [34, 37, 35, 36],
       [34, 41, 43, 44],
       [35, 39, 37, 38],
       [39, 43, 41, 42]], dtype=int32)

Looking at the topology impropers:

>>> top._topology.impropers.values
[(4, 12, 10, 11), (14, 27, 25, 26), (19, 23, 22, 24), (29, 47, 45, 46), (35, 43, 34, 31),
...]

(35, 43, 34, 31) is now (31, 34, 43, 35).

Code to reproduce the behavior

Show us how to reproduce the failiure. If you can, use trajectory files from the test data.

import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PRM

u = mda.Universe(PRM)
print(u.atoms[[31,34,43,35]].impropers.indices)
print(u._topology.impropers.values)

Currently version of MDAnalysis

• Which version are you using? (run python -c "import MDAnalysis as mda; print(mda.__version__)") 0.20.1
• Which version of Python (python -V)? 3.7.3
• Which operating system? MacOS

[lilyminium]

There's more to this...

As has been discussed before:

CHARMM (PSF) and GROMACS (ITP/TPR) have the central atom in an improper first, while AMBER (PRMTOP) has it third. None of the related parsers rearrange Improper atom ordering. I think this will complicate outputting to a standard topology structure like Parmed in #2404.

[richardjgowers]

So MDA currently uses this definition of improper: http://cbio.bmt.tue.nl/pumma/uploads/Theory/improper.png

This looks like what gromacs also uses:

http://manual.gromacs.org/documentation/2020-beta2/reference-manual/functions/bonded-interactions.html?highlight=improper#improper-dihedrals

So the angle between the plane of IJK and plane of JKL. With this definition the order can be reversed.

Amber looks like it didn't have a standard: ?
openmm/openmm#220

So I think reversing it seems harmless, but I'm not sure why it was reversed...
https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/topology/TOPParser.py#L583

[IAlibay]

Sorry, trying to catch up here, but I'll happily take the blame for any non-standard improper ordering on reading parm7 files in mda.

If this helps in any way, I believe the attempt in the reader was just to read the atom indices for each dihedral (be it improper or not) purely based on the order that they are written in the parm7 file.

[richardjgowers]

Ah I just reread all this, so it looks like if Amber has the central atom as k, then we need to rearrange to get the same angle.

[IAlibay]

So on the note of reversal of the i,j,k,l order, this discussion might be worth noting: https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2013-April/080824.html

Appending the central atom as the first entry isn't too much of a change, but going by this comment openmm/openmm#220 (comment) consistently ordering the remaining atoms is not immediately obvious to me (I am probably just missing something obvious).

I was thinking just following whatever logic parmed takes when converting between parm7 to top files might be the best idea here, but I'm not actually sure that any re-ordering actually takes place.

Creating an { ACE PHE NME } peptide in tleap, the following improper dihedral is written in the parm7 file:

39      60     -66     -69

which relates to the following parmed-generated human readable entry (just divide by 3 and add one to the numbers to get the corresponding atom number):

I      14   CG (  CA)       21  CE2 (  CA)       23  CD2 (  CA)       24  HD2 (  HA)

now converting this via parmed (3.2.0) to a .top file gives you the following dihedral entry:

[ dihedrals ]
;    ai     aj     ak     al funct         c0         c1         c2         c3         c4         c5
     14     21     23     24     4  180.0000771  4.6024000  2

Again, I should point out that is has been a long day, and I am most likely misunderstanding something.

phedi.zip

[lilyminium]

Looking through Parmed a bit, it puts central atoms first for harmonic improper dihedrals, but third for periodic impropers (which all of AMBER's are).

GROMACS periodic impropers (dihedrals type 4) are identical to periodic propers:

Proper dihedral angles are defined according to the IUPAC/IUB convention, where φ is the angle between the ijk and the jkl planes, with zero corresponding to the cis configuration (i and l on the same side).

Interestingly, according to the DIHEDRALS_INC_HYDROGEN section of the prmtop format:

NOTE: The first atom has an index of zero. Since 0 cannot be negative and the 3rd and 4th atom indexes are tested for their sign to determine if 1-4 terms are calculated, the first atom in the topology file must be listed as either the first or second atom in whatever torsions it is defined in. The atom ordering in a torsion can be reversed to accommodate this requirement if necessary.

So it looks like impropers are symmetric, oops. But it's probably worth tracking which are harmonic and which are periodic, at least for ParmEd i/o.

[IAlibay]

So not relevant to the the discussion on impropers, but on a larger side note on AMBER topologies here (considering the whole ParmEd i/o), beyond impropers we have a wider issue that some of the fancier terms are not supported / read by MDAnalysis.

Of particular note here is that CMAP terms are not supported.
I have been slowly trying to put together a plan to extend my previous work on the reader to also support CHAMBER-style topologies, so that should be doable (hopefully proximity to Xmas will offer me time to tie that up and the netcdf readers/writer). However, ff19SB implements a different CMAP standard to the charmm approach (see ParmEd/ParmEd#1066 (comment)).

For now ParmEd doesn't support ff19SB-style CMAP terms, however we might need to implement a means at topology read to block/warn users that load say ff19SB parm7 files (I haven't tried this yet, but one of my students claims it just works), and then try to pass it to ParmEd, that not all the terms in the original parm7 are accounted for.

[lilyminium]

Closed due to not being real issue.