assert_dihedral_params incorrect behavior for layered dihedrals (e.g., GAFF)

[rsdefever]

Describe the behavior you would like added to Foyer

Currently, the _check_dihedrals function of forcefield.py verifies that the total number of dihedrals identified matches the number of parameterized dihedrals. This is incorrect behavior if multiple dihedrals potentials are placed on the same physical dihedral (i.e., the same i-j-k-l atoms). Foyer exits with a Parameters have not been assigned to all proper dihedrals. Total system dihedrals: XX, parameterized dihedrals: YY. message, even if YY > XX.

Describe the solution you'd like

Foyer should check that for each physical dihedral i-j-k-l atoms, that there is at least one parameterized dihedral.

Describe alternatives you've considered

We could just ensure that the number of parameterized dihedrals is greater than or equal to the number of system dihedrals. However, this solution may fail to catch unparameterized dihedrals in the case where some physical dihedrals had multiple potentials while other physical dihedrals had zero potentials.

[ahy3nz]

structure.join_dihedrals()

If dihedrals have multiple parameters layered on them, they parameters applying to the same dihedral get grouped together into a single DihedralType. I think this might help address the issue where the number of parametrized dihedrals exceeds the number of system dihedrals.

[mattwthompson]

I think this might help address the issue where the number of parametrized dihedrals exceeds the number of system dihedrals.

Is this necessarily a case we want to avoid, though? The three _check_ functions only look at equality. It is a problem in ParmEd world if there are more dihedral types than dihedrals?

[ahy3nz]

I might need to test this out further, but let's say you had a system with 3 chemically-identical dihedrals, but each dihedraltype has 4 PeriodicTorsions to it. When this gets initially converted into parmed, structure.dihedral_types will have 4 entities, reflecting the 4 different periodic torsion parameters even though they all belong to the same dihedral "key". If we call structure.join_dihedrals(), the 4 DihedralTypes will get converted into 1 DihedralTypeList object, and structure.dihedral_types should have len=1

[rsdefever]

@ahy3nz nice catch. A quick test suggests this is indeed the behavior.

This would affect the writers as they would need to loop through the DihedralTypeList to extract dihedral parameters.

Does it make sense to apply structure.join_dihedrals() to the final parameterized structure or should we copy the data off to a temporary temp_structure and then do temp_structure.join_dihedrals() when we check for the correct number of parameterized dihedrals in _check_dihedrals?

[ahy3nz]

structure.join_dihedrals() would only affect periodic torsions (I think), and the only place it would really affect us down-the-line in the MoSDeF code would be in the lammpswriter, but we already call structure.join_dihedrals() in the lammpswriter.

Anyways, point is, I think it's okay to call join_dihedrals on the final, parametrized structure that we return - I don't think we're hurting anyone with this, and there's the added benefit of not having to do a copy

[rsdefever]

An update on this issue:

It appears that if the force constant for a dihedral is 0.0, then foyer (really openmm) doesn't add the dihedral to the parameterized list. An example of this is GAFF where there is a dihedral defined as: class1="" class2="c3" class3="ca" class4="". The force constant is 0.0. When I run a relevant molecule through foyer, it add no parameters for a ca-c3-ca-ca dihedral.

Unfortunately, this means that the number of parameterized dihedrals will sometimes be less than the number of physical dihedrals.

[mattwthompson]

I recall Andrew or somebody pointing this out some time ago, but I don't recall our workaround. I don't think there was a clear workaround while using OpenMM, except that we intend to be forgiving with that assert (compared to the case of missing bonds or angles, which default to be an error). Unfortunately, this doesn't really help you do your bookkeepping with GAFF or any other force field with zero-force dihedrals.