Unexpected dihedral parameters from Open Force Field Toolkit

[martimunicoy]

Description

We noticed that some of the parameters assigned to torsions by the Open Force Field Toolkit are unexpected. The current implementation of offpele expects the following values for torsional parameters:

• phase constant: either 0.0 or 180.0 degrees
• periodicity: one of (1, 2, 3, 4, 6)

See issue #33 for more details.

While the requirement in the periodicity value is trivial, the phase constant does really need to be either 0.0 or 180.0 degrees. That is because of PELE's current dihedral equation (explained here), which is:

Since, there is not a phase constant term, we need to use the prefactor which can equal 1 or -1. The sign determines whether the shape of the torsion behaves as in Open Force Field with a phase constant equal to 0.0 or 180.0 degrees. Thus, in case that the Open Force Field Toolkit works with phase values other than 0.0 and 180.0 degrees we will need to change the equation of dihedrals in PELE.

[martimunicoy]

Some cases where we have obtained unexpected parameters in torsions:

Case 1

from openforcefield.topology import Molecule, Topology
from openforcefield.typing.engines.smirnoff import ForceField

mol = Molecule.from_smiles('[H]c1c(sc([n+]1[H])[C@@]2(C(C(=O)N2[H])([H])[H])C([H])([H])[H])[H]')
ff = ForceField('openff_unconstrained-1.2.1.offxml')
topology = Topology.from_molecules([mol])
parameters = ff.label_molecules(topology)[0]

• Molecule: [H]c1c(sc([n+]1[H])[C@@]2(C(C(=O)N2[H])([H])[H])C([H])([H])[H])[H]
  
• Unexpected propers:
  • t26, phase=270.0 deg, periodicity=4
  • t26, phase=270.0 deg, periodicity=2
  • t26, phase=270.0 deg, periodicity=1

Case 2

from openforcefield.topology import Molecule, Topology
from openforcefield.typing.engines.smirnoff import ForceField

mol = Molecule.from_smiles('[H]c1c(onc1C2(C(C2([H])[H])([H])[H])C(=O)O[H])[H]')
ff = ForceField('openff_unconstrained-1.2.1.offxml')
topology = Topology.from_molecules([mol])
parameters = ff.label_molecules(topology)[0]

• Molecule: [H]c1c(onc1C2(C(C2([H])[H])([H])[H])C(=O)O[H])[H]
  
• Unexpected propers:
  • t42, phase=320.0 deg, periodicity=2
  • t42, phase=320.0 deg, periodicity=2

Case 3

from openforcefield.topology import Molecule, Topology
from openforcefield.typing.engines.smirnoff import ForceField

mol = Molecule.from_smiles('[H]c1c(c(n(n1)S(=O)(=O)C([H])([H])C([H])([H])[H])[H])O[H]')
ff = ForceField('openff_unconstrained-1.2.1.offxml')
topology = Topology.from_molecules([mol])
parameters = ff.label_molecules(topology)[0]

• Molecule: [H]c1c(c(n(n1)S(=O)(=O)C([H])([H])C([H])([H])[H])[H])O[H]
  
• Unexpected propers:
  • t142, phase=90.0 deg, periodicity=3

[j-wags]

Thanks for the excellent writeup, @martimunicoy.

I've checked in the original FFs and those phase values are indeed intentional. I apologize -- I was wrong about this previously. I'm happy to see it's a small number of terms, but I had thought that there were none at all :-/

Would you be up to chat about this soon? In the next year or two, OpenFF is anticipating adding virtual sites (off center charges) to our FFs. It may be worth looking at having PELE energy evaluations in OpenMM, since we work with their development team, and will ensure that it can handle any terms in the force fields we release.

[martimunicoy]

Thanks for your comment, @j-wags.

I guess we will have to modify the current implementation of the torsion equation in PELE. I have already started looking into a workaround. Besides, once you release a FF with virtual sites we will have to implement them into PELE too. Although it seems you do not plan to release a FF with virtual sites anytime soon, we can already start looking at it.

We can talk about it this current week, if it suits you (either Wednesday, Thursday or Friday).

Thank you again!:)

[j-wags]

Great. I'll coordinate the meeting on Slack :-)