Fix remaining issues with full drugs test

[kyleabeauchamp]

So it appears that there are still some issues with the all drugs test. I have confirmed that they issues appear to be:

1. Primarily energy differences
2. Accompanied by differences in forces (e.g. not just a trivial energy offset)
3. Are probably caused by some of the "nastier" force components (e.g. non-obvious origin)
4. Are independent of whether we use the OChem or the MDTraj loading path.

[jchodera]

I suggest we proceed by collecting a list of all molecules with energy differences and sorting by number of atoms, taking the simplest cases (fewest atoms) first.

[kyleabeauchamp]

So the smallest failure suggests that this is related to known issues with the ordering of atoms in torsions. There is essentially a single torsion term that disagrees:

(2, 11, 4, 5, 2, Quantity(value=3.14159265359, unit=radian), Quantity(value=4.6024, unit=kilojoule/mole))

(11, 2, 4, 5, 2, Quantity(value=3.141594, unit=radian), Quantity(value=4.6024, unit=kilojoule/mole))

We can possibly avoid raising errors on these issues, but perhaps we should check with Jason about what he thinks. My concern is that there is a nonzero force difference between the two permutations:

In [4]: delta
Out[4]: 
Quantity(value=array([[  7.80569666e-05,  -2.94378036e-05,  -3.59438541e-04],
       [  6.44469299e-05,   6.29361630e-04,  -1.82332906e-03],
       [  6.47828517e+00,  -7.58939674e+00,  -4.05393155e+01],
       [  7.72669684e-04,   1.59426371e-04,   1.21480494e-04],
       [ -3.40229215e+01,   4.48256524e+01,   4.01034047e+01],
       [  2.50570050e+01,  -3.92174250e+01,  -1.08051860e+01],
       [ -1.83442232e-03,  -5.47431499e-03,   3.94541665e-04],
       [  1.80442957e-03,   1.90104896e-03,   3.01372908e-03],
       [  5.48325661e-03,  -3.79568526e-03,  -4.02752879e-04],
       [ -5.97974115e-04,   9.56827477e-04,  -5.10031847e-04],
       [  2.48774144e+00,   1.99057251e+00,   1.12428893e+01],
       [ -4.46687181e-04,  -1.34580806e-03,  -6.11047927e-04],
       [  1.34711418e-03,   7.38158402e-05,  -6.19783394e-04]]), unit=kilojoule/(nanometer*mole))

Note the four entries with values that are ~1E4 larger than the others...

[kyleabeauchamp]

It's worth noting that this permutation is happening for torsions that involve 4 distinct atoms--previously I thought this permutation issue was limited to torsions with 3 atoms...

[jchodera]

This is not good. I believe our previous discussions found that OpenMM's ForceField class handles torsion assignments in a manner that differs from AmberTools LEaP. I had tried to convince Peter that this was problematic, but we didn't have the "killer evidence" at that point. But this seems to be it.

I think the only way to fix this and still integrate with the ForceField ffxml framework is to make sure ForceField assigns torsions in exactly the same way LEaP does.

[jchodera]

See this thread: https://github.com/SimTk/openmm/issues/220

[jchodera]

Are the examples above of improper torsions, or proper dihedrals?

[kyleabeauchamp]

Improper

[jchodera]

We may need help from @swails to decode what LEaP is actually doing.

[kyleabeauchamp]

I'll email him to watch this thread...

[swails]

So tleap's handling of impropers is admittedly strange and currently has a glaring problem for efforts like gaff2xml. Improper torsions in Amber are treated with the same potential as proper torsions, but they are not assigned the same way.

Improper torsions are stored (and matched) inside tleap by fixing the 3rd atom as defined in the frcmod or parm.dat files, and then arranging the remaining 3 atoms in alphabetical order by atom type. So an improper defined between atom types c3, c1, cc, and ca in the parm.dat file will be rearranged into an improper of the form c1-c3-cc-ca (since cc is fixed and 1 < a).

This is all well and good until you have an improper in which two of the atoms that are not fixed share the same type. For instance, ca-ca-cc-c1 will become c1-ca-cc-ca, but how the two ca's map to the atoms that make up a particular torsion is not well-defined. For a "true" improper (like the kind CHARMM has), I don't think that matters. Since Amber treats it like a proper, though, it does matter, and you can get slightly different energies and forces depending on what (arbitrary) ordering is picked.

I talked with Dave Case a while ago about this very topic, and I think what we agreed upon was that atom names should be used to specify the order of any atoms that share a common type. I have not yet implemented that in tleap, though. I will try to implement that tomorrow for the sake of making comparisons with gaff2xml easier.

In practice, this ordering thing is not a big deal. Regardless of the ordering, the improper torsion does its job (it maintains planarity). So in the context of a gaff2xml-like project, this is a "serious" flaw. In the context of a real-world simulation, it can be safely neglected.

[jchodera]

Thanks, @swails!

Can't we just add an alphabetizing step at some point inside gaff2xml or ForceField to ensure we get the same results as LEaP? Even if there is little practical difference, energy difference between programs worry many people...

[kyleabeauchamp]

For now, it looks like I can slightly increase the tolerance for the torsional energies and merge my PR, and get travis running.

On Wed, Jun 4, 2014 at 11:17 PM, John Chodera notifications@github.com wrote:

Thanks, @swails https://github.com/swails!

Can't we just add an alphabetizing step at some point inside gaff2xml or ForceField to ensure we get the same results as LEaP? Even if there is little practical difference, energy difference between programs worry many people...

— Reply to this email directly or view it on GitHub https://github.com/choderalab/gaff2xml/issues/24#issuecomment-45177114.

[jchodera]

For now, it looks like I can slightly increase the tolerance for the torsional energies and merge my PR, and get travis running.

That's what we did last time, and we still found problems. I'd rather we find a real solution.

[kyleabeauchamp]

I agree, but the "real" solution will involve coordination that's outside the scope of the current PR--we might be waiting for an AmberTools patch. In the meantime, we need to make progress on the other outstanding issues.

[jchodera]

I think the solution here is just to introduce an atom sorting step to match tleap behavior. Wouldn't that do the trick?

If we wanted the fix to be local to just gaff2xml-derived ffxml files, we could even do the sorting step just in gaff2xml...

[jchodera]

On second thought, if @swails is actually interested in fixing tleap for an upcoming release, we should probably forge ahead rather than trying to match its behavior...

[swails]

Just to clarify, the "fix" I'm talking about does not remove the need to sort atoms 1, 2, and 4 of impropers in alphabetical order by atom type. The "fix" merely aims to solve instances in which that sorting is ambiguous by sorting degenerate atom types by atom names.

Imagine a phenol group where you want an improper between the O and 3 of the aromatic carbons to maintain planarity. You'll get the 4 atom types as something like ca-ca-cb-oh, where the atom names might be C3-C5-C4-O.

Because C3 and C5 have the same atom type, tleap can do one of two things -- it can define that improper as C3-C5-C4-O or it can define the improper as C5-C3-C4-O, and both correspond to a "properly sorted" improper.

In a practical example, this bit us when we were trying to do a TI calculation and tleap assigned the same improper different orders in the two end states. Because the Amber implementation requires conserved parameters to be exactly identical, this was a problem. (This is actually when I discovered the problem in the first place).

All I was suggesting here was to go ahead and sort degenerate types by atom names when necessary because that is the fix that will match tleap eventually (after I implement that fix). I'm not sure if that fix will be turned into an update or not (I'll have to consult Dave about that), but I can certainly push it to ambermini.

[jchodera]

Got it. It does sound like it would be optimal to put a fix in to ambermini and have gaff2xml duplicate this behavior and match this behavior with gaff2xml/ForceField to validate our implementation.

I do worry that atom names may not be optimal. It would fix the TI issue you had, but molecules with the same types and charges but different atom names will have different energies, which feels like it has the potential to lead to confusion down the line.

It's probably the best solution, but is there another unambiguous way to assign impropers? Amber avoids ambiguities with multiple dihedral by assigning all dihedrals that match along a bond, then dividing by the number of dihedrals. Could something like that work for impropers?

[swails]

I do worry that atom names may not be optimal. It would fix the TI issue you had, but molecules with the same types and charges but different atom names will have different energies, which feels like it has the potential to lead to confusion down the line.

I don't think it will really cause that much (if any) confusion. Don't forget that for the same ligand, people are going to use the same library files (and therefore must have the same names for the same atoms or LEaP will fail). If you have a different number of atoms or a different connectivity, single-point force field energies are not comparable, anyway. Furthermore, enantiomers can (and should) use the same residue library template (just make sure the input PDB file has the necessary chirality defined or create a new residue from the mirror image of the original residue template).

It's probably the best solution, but is there another unambiguous way to assign impropers?

I don't know. I'm curious about how CHARMM defines impropers when writing the PSF file. I don't have access to the CHARMM source code and the psfgen source code in vmd severely lacks comments (at least with regards to how impropers are found and defined). There are only 2 ways that I can think of that will guarantee a unique ordering: sorting by the ordering of the atoms in the template and sorting by atom name (assuming that, unlike atom types, numbers and names cannot be degenerate). Anything else is necessarily arbitrary.

Amber avoids ambiguities with multiple dihedral by assigning all dihedrals that match along a bond, then dividing by the number of dihedrals. Could something like that work for impropers?

I disagree here. Proper dihedrals are always well-defined by the 3 bonds making up the torsion. The torsion is defined around the second bond, which is between atoms 2 and 3. Atom 1 is bonded to 2, and 4 to 3. The only two orders that work for this definition are 1-2-3-4 and 4-3-2-1 -- and those two orderings give rigorously identical torsion angles. Improper torsions do not have this serial bond pattern -- one of the 4 atoms is connected to all other atoms (and none of the other 3 atoms are connected to any other). Only the central atom can be well-defined with respect to the other atoms

The division factor you refer to is a stylistic choice that is applied exclusively to dihedrals defined with wild-cards. The energy is what you would get with a torsion scan of a set of molecules containing that wild-card dihedral, and the division factor is simply how many of those dihedral terms you expect to have. This scaling factor can be divided into the force constant and set to 1 everywhere (it has no impact on ordering or determining uniqueness).

[swails]

So I locally patched my version of tleap to order degenerate atom types by the atom name, only to find out that tleap already sorts these degenerate types by the ordering of the atoms in the unit.

It also appears that the atom ordering in the topology file always matches the atom ordering in the library file defining the unit (so tleap will reorder atoms from an input PDB file to match the ordering specified in either the ligand mol2 file or the standard residue libraries). As a result, improper torsions are already unambiguously ordered by tleap (just sort by atom type, and sort degenerate atom types by atom indexes inside the mol2 file).

I believe I have an example where tleap does not behave this way, but that is a bug that needs to be fixed. If you guys have more examples where this is true, please forward them to me so I can begin to put together some kind of bug report.

However, since tleap already defines a working convention, there is no way I'll be able to push another one on them (I just have to get the current convention working rigorously if it isn't already).