Atom ordering in AMBER improper torsions does not match LEaP

[jchodera]

@kyleabeauchamp and I are running into some weird issues with GAFF-to-ffxml conversion which seems to be due to the atom ordering in improper torsions assigned by ForceField not precisely matching the order they would be assigned by AmberTools LEaP. These differences result in small (but nonzero) energy differences in most cases, but make it very difficult to validate the conversion and may account for larger differences under different circumstances.

Closer inspection of AmberTools LEaP code suggests that the difference may be primarily due to the fact that LEaP seems to alphabetically sort the improper torsions before applying the torsions, supposedly to "speed up" the search:

amber12/AmberTools/src/leap/src/leap/parmSet.c:251

/*                                                                                                                                                                                                                                                                                                                             
 *      zParmSetOrderImproperAtoms                                                                                                                                                                                                                                                                                             
 *                                                                                                                                                                                                                                                                                                                             
 *      Author: Christian Schafmeister (1991)                                                                                                                                                                                                                                                                                  
 *                                                                                                                                                                                                                                                                                                                             
 *      Order the improper atom _types_ into alphabetical order                                                                                                                                                                                                                                                                
 *      to speed up searches.                                                                                                                                                                                                                                                                                                  
 *                                                                                                                                                                                                                                                                                                                             
 *      The proper order is for the wild cards to be first                                                                                                                                                                                                                                                                     
 *      and the rest of the types to be in alphabetical order.                                                                                                                                                                                                                                                                 
 *      Only the first, second, and fourth atoms are affected,                                                                                                                                                                                                                                                                 
 *      the third atom is ALWAYS the central atom.   
 */                                                                                                                                                                                                                                                                                                                            

I wonder if doing this at the ffxml conversion step would be sufficient to eliminate this discrepancy, and could potentially allow you to remove some of the weird AMBER-specific juju in the assignment of impropers in forcefield.py.

[jchodera]

I think this is still an outstanding problem causing some discrepancies between OpenMM and LEaP-parameterized molecules.

[peastman]

I don't believe it's possible to fix this. Last time I checked, the way Amber ordered improper torsions was not fully specified and not even internally consistent. For example, you could have the same residue appear twice in a protein, and in one residue a torsion would be ordered one way, and in the other residue it would be ordered a different way.

[jchodera]

@swails may have insight into this, but I thought it just involved some alphabetical sorting of torsions that could be recapitulated. It may not be worthwhile to match, however.

[swails]

Amber improper torsions are sadly done in a way that is not well-defined. The central atom comes third, and the other three are in alphabetical order by atom type. When two atoms involved in the improper torsion have the same type, the ordering becomes undefined, and indeed there are instances where leap will assign two different copies of the same ligand with different orders of impropers that are still consistent according to their sorting rules.

I have plans to go into tleap and order by atom type, then order degenerate atoms by atom name, but I haven't done this yet. atom name should have been done from the beginning, since that choice does not suffer from undefined behavior (all atom names in a residue must be unique), and if you change the atom types as part of refactoring a force field, you might unknowingly change the ordering of the impropers.

At the end of the day, these torsions are not rigorously derived, anyway -- they just exist to keep things planar. And they'll do that even if the order isn't exactly the same between programs. So I don't think it's a big deal.

[jchodera]

But didn't we find there were significant differences when running a set of OpenEye-prepped FDA-approved drugs through what was formerly gaff2xml?
choderalab/openmoltools#24

[jchodera]

Fun idea: Maybe we should test how much this impacts physical quantities like hydration free energies?

[jchodera]

I think this as been addressed by #16 and #193 above.