PeriodicTorsionGenerator behavior with GAFF- frcmod ffxml use

[rafwiewiora]

Tagging @jchodera @swails @peastman

I'm not sure if we had talked about this before, I hope not given how much it confused me...

I was making a test system for GAFF and I chose imatinib, mol2: https://github.com/choderalab/openmoltools/blob/master/openmoltools/chemicals/imatinib/imatinib.mol2

So we have a:

• converted gaff.dat --> gaff.xml
• frcmod made for imatinib with parmch2, converted to ffxml - imatinib_frcmod.xml. There are no double calls to parameters between gaff.xml and imatinib_frcmod.xml, except that some dihedrals using wildcards are also re-specified without wildcards in the frcmod. (so I guess these are double calls after all)
• imatinib.lib made with LeAP from the mol2, converted to xml - imatinib.xml

and consider two calls to ForceField:

A. ff = ForceField('gaff.xml', 'imatinib_frcmod.xml', 'imatinib.xml')

b. ff = ForceField('imatinib_frcmod.xml', 'gaff.xml', 'imatinib.xml')

These two give different energies, particularly on the Impropers: A:

[('HarmonicBondForce', 26.14324505045579),
 ('HarmonicAngleForce', 22.300047364727916),
 ('PeriodicTorsionForce', 36.33697380921836),
 ('PeriodicTorsionForce', 0.4123695363923616),
 ('NonbondedForce', -852.4329694509506),
 ('CMMotionRemover', 0.0)]

b:

[('HarmonicBondForce', 26.14324505045579),
 ('HarmonicAngleForce', 22.300047364727916),
 ('PeriodicTorsionForce', 36.3369079713571),
 ('PeriodicTorsionForce', 0.14815193784295388),
 ('NonbondedForce', -852.4329694509506),
 ('CMMotionRemover', 0.0)]

there are also different number of torsions in the systems: (sys_a.getForce(x).getNumTorsions()), A: 207, B: 183.

So what I understand and think is happening:

• we have two <PeriodicTorsionForce> xml children, so ForceField will create two independent PeriodicTorsionGenerator 's.
• on calling createSystem() both of these generators will call createForce (https://github.com/pandegroup/openmm/blob/master/wrappers/python/simtk/openmm/app/forcefield.py#L1330) and do assignments for all torsions in the system data independently
• so the problem here seemed to be the overlapping dihedrals that use wildcards from gaff.xml and the same dihedrals but specified with all 4 types from frcmod? This would need to be solved by doing a 'join' on generators prior to creating forces, or maybe easier on the xml trees.

But what I find weird is that the situation with two generators seems additive, so the order of files in ForceField shouldn't matter - but the numbers of torsions and energies in the systems are different. What am I missing here? What adds to my confusion is the fact that if I do ff = ForceField('gaff.xml', 'imatinib.xml') (no frcmod) the energies come out the same as in case B:

[('HarmonicBondForce', 26.14324505045579),
 ('HarmonicAngleForce', 22.300047364727916),
 ('PeriodicTorsionForce', 36.3369079713571),
 ('PeriodicTorsionForce', 0.14815193784295388),
 ('NonbondedForce', -852.4329694509506),
 ('CMMotionRemover', 0.0)]

[rafwiewiora]

The files are here:

gaff.xml: https://gist.github.com/rafwiewiora/c48d6ac828263ae19642 imatinib_frcmod.xml: https://gist.github.com/rafwiewiora/0d5cc642ea73d5c12934 imatinib.xml: https://gist.github.com/rafwiewiora/8ddf558b76c4f7d7cfca

[jchodera]

Pinging @peastman as well.

I hadn't realized each <PeriodicTorsionForce> would create a new PeriodicTorsionGenerator.

• @peastman: Is this the intended behavior?
• Is there some way we could provide for <PeriodicTorsionForce> groups to be joined across different files?
• Are there opportunities where this would be undesirable, and we would want to prevent this from happening?

Perhaps some additional attribute in the XML tag could specify the group that the forces should be merged across different files into?

[rafwiewiora]

Ok, I get now why ff = ForceField('imatinib_frcmod.xml', 'gaff.xml', 'imatinib.xml') has the same result energy-wise as ff = ForceField('gaff.xml', 'imatinib.xml'). The torsions in the frcmod, if called before gaff.xml, which defines atom types, are ignored because the types are not defined yet.

After we decide how to solve this, can I please go through this code and add the much needed warnings everywhere? I'm getting stuck for hours over and over again, because these things just happen and say nothing.

So it's just the wild cards - not wild cards duplication in two different generators. I'm just gonna propose we take the https://github.com/pandegroup/openmm/pull/1373 that has been hanging for a moment and go further with: don't duplicate generators and put everything in one per force to solve the wild card problem. This should all be compatible with the overriding for templates we decided on - I'm hoping to code this up this week too.

EDIT: I should clarify that the wild cards - not wild cards issue still clearly remains, because the energies returned by AMBER and converted ffxml-OpenMM do not pass validation. But I guess that was already obvious from seeing how the ForceField code would behave. Just wanted to point this out, because I only talked about the different energies of different ForceField calls in the first comment.

[peastman]

This is a good example of why it's a bad idea to have the results depend on the order in which you list files! I agree with your recommendation. We should create only a single generator (some of them already do that, such as NonbondedForceGenerator), and we should process the input files in a way that doesn't make the results depend on the order they were listed in.

[rafwiewiora]

Closing this as problem is resolved in https://github.com/pandegroup/openmm/pull/1373, just waiting for merging.