Sometimes there are no angle terms in the forcefield for a given angle proposal

[pgrinaway]

The FFAllAngleGeometryEngine proposals work like this:

Assume we have a topological torsion connected like:

A--B--C--D

and we are trying to propose a position for atom D. It's not necessarily the case that the forcefield will have a term for the torsion defined by A-B-C-D, and this is not an issue in the GeometryEngine. However, it also assumes that the angle defined by B-C-D will always have a definition in the forcefield; this is usually but not necessarily true. When such an issue is encountered, the GeometryEngine fails. I can think of several ways to solve this:

• Propose from a reference distribution with a weak force constant
• Propose from a reference distribution parameterized by the angle generated by Omega and a weak force constant
• Something else?

I'm leaning toward the second option, but am open to other suggestions.

cc: @jchodera

[jchodera]

Can you give an example of when this fails? All angles should have angle terms in well-formed forcefields. Is this an improper?

[pgrinaway]

Can you give an example of when this fails? All angles should have angle terms in well-formed forcefields. Is this an improper?

OK, getting together an example now. Is this true? I was talking to some others who suggested they also ran into this issue.

[pgrinaway]

Is this an improper?

The torsion is always proper since I specifically look for topological torsions, not torsion force terms.

[jchodera]

We would have to see some examples to diagnose further, but my understanding is that missing angles can only be the result of a forcefield typing failure.

[jchodera]

The only thing I can think of causing this is that the frcmod files generated by parmchk2 are either missing angles in some cases or the frcmod files are not being read in, which would be a significant issue.

[pgrinaway]

We would have to see some examples to diagnose further, but my understanding is that missing angles can only be the result of a forcefield typing failure.

This lists one result for me:

smiles= "[H]c1c(c(c(c(c1[H])[H])OC(=O)[H])[H])[H]"
from openeye import oechem, oeiupac
import parmed
mol = oechem.OEMol()
oechem.OESmilesToMol(mol, smiles)
iupac_name = oeiupac.OECreateIUPACName(mol)
from io import StringIO

def createSystemFromIUPAC(iupac_name):
    """
    Create an openmm system out of an oemol
    Parameters
    ----------
    iupac_name : str
        IUPAC name
    Returns
    -------
    molecule : openeye.OEMol
        OEMol molecule
    system : openmm.System object
        OpenMM system
    positions : [n,3] np.array of floats
        Positions
    topology : openmm.app.Topology object
        Topology
    """

    # Create OEMol
    from perses.tests.utils import createOEMolFromIUPAC, extractPositionsFromOEMOL
    molecule = createOEMolFromIUPAC(iupac_name)

    # Generate a topology.
    from openmoltools.forcefield_generators import generateTopologyFromOEMol
    topology = generateTopologyFromOEMol(molecule)

    # Initialize a forcefield with GAFF.
    # TODO: Fix path for `gaff.xml` since it is not yet distributed with OpenMM
    from simtk.openmm.app import ForceField, HBonds
    from perses.tests.utils import get_data_filename
    gaff_xml_filename = get_data_filename('data/gaff.xml')
    forcefield = ForceField(gaff_xml_filename)

    # Generate template and parameters.
    from openmoltools.forcefield_generators import generateResidueTemplate
    [template, ffxml] = generateResidueTemplate(molecule)

    # Register the template.
    forcefield.registerResidueTemplate(template)

    # Add the parameters.
    forcefield.loadFile(StringIO(ffxml))

    # Create the system.
    system = forcefield.createSystem(topology, removeCMMotion=False, constraints=HBonds)

    # Extract positions
    positions = extractPositionsFromOEMOL(molecule)

    return (molecule, system, positions, topology)

mol, sys, pos, top = createSystemFromIUPAC(iupac_name)

angle_force = sys.getForces()[1]
bond_force = sys.getForces()[0]
angle_list = []
for idx in range(angle_force.getNumAngles()):
    [a1, a2, a3, theta, k] = angle_force.getAngleParameters(idx); angle_list.append([a1, a2, a3])

bond_list = []

for idx in range(bond_force.getNumBonds()):
    [a1, a2, r, k] = bond_force.getBondParameters(idx)
    bond_list.append([a1, a2])

s = parmed.openmm.load_topology(top, sys)

for atom in s.atoms:
    angles_to_try = []
    bond_partners = atom.bond_partners
    for bond_partner in bond_partners:
        for double_partner in bond_partner.bond_partners:
            angles_to_try.append([atom, bond_partner, double_partner])

for angle in angles_to_try:
    atom1_angles = set(angle[0].angles)
    atom2_angles = set(angle[1].angles)
    atom3_angles = set(angle[2].angles)

    common_angles = atom1_angles.intersection(atom2_angles, atom3_angles)

    if len(common_angles) == 0:
        print(angle)

[<Atom H6 [14]; In phenyl formate 0>, <Atom C7 [6]; In phenyl formate 0>, <Atom O1 [8]; In phenyl formate 0>]

[pgrinaway]

It doesn't seem like there's an angle for type (h5, c, os) in gaff.xml or in the generated ffxml

[jchodera]

@pgrinaway : If your leap session gives you parameters for this angle, can you post the leap input files and script necessary to reproduce and inspect the prmtop?

[jchodera]

I suspect this will be addressed by the networkx-based torsion selection scheme discussed in #508