Implement GAFF/AM1-BCC residue template generator for simtk.openmm.app.ForceField

[jchodera]

This is the beginning of the implementation of a GAFF/AM1-BCC residue template generator for small molecules for use with the new simtk.openmm.app.ForceField ability to hook in residue template generators.

Example usage:

from simtk.openmm.app import ForceField, NoCutoff
# Load my PDB file containing a ligand
pdb = PDBFile(pdb_filename)
# Create a ForceField object, loading GAFF parameters.
forcefield = ForceField('amber99sbildn.xml', 'tip3p.xml', 'gaff.xml')
# Register the GAFF residue template generator.
from openmoltools.forcefield_generators import gaffTemplateGenerator
forcefield.registerTemplateGenerator(gaffTemplateGenerator)
# Create an OpenMM system, calling antechamber as needed to parameterize small molecules.
system = forcefield.createSystem(pdb.topology, nonbondedMethod=NoCutoff)

Note that this does not yet work with residues that have external bonds (like unnatural amino acids or covalently-attached ligands), but that feature is coming soon.

There are several hurdles to be overcome:

• Test and debug bond order perception
• Extra parameters generated by parmchk2 need to be incorporated
• Find a solution for molecules for multiple chiral centers

[jchodera]

Note that this feature also currently requires OpenEye (for writing Tripos mol2 files) and ParmEd (for reading GAFF-format mol2 files).

[jchodera]

This nearly works, but we need to fix an issue with scee/scnb in ParmEd:
ParmEd/ParmEd#568

[jchodera]

cc: choderalab/perses#67

It looks like we might be able to go either from an OpenMM Topology Residue (using Antechamber's bondtype bond perception) or an OEMol molecule.

Both the pre-generation of residue templates:

forcefield = ForceField('amber99sb.xml', 'tip3p.xml', 'gaff.xml')
[template, ffxml] = generateResidueTemplate(molecule)
forcefield.registerResidueTemplate(template)
forcefield.loadFile(StringIO(ffxml))
system = forcefield.createSystem(topology)

and the residue template generator plugin

# Initialize ForceField
forcefield = ForceField('amber99sb.xml', 'tip3p.xml', 'gaff.xml')
from openmoltools.forcefield_generators import gaffTemplateGenerator
forcefield.registerTemplateGenerator(gaffTemplateGenerator)
# Now, when you need to create a new System...
system = forcefield.createSystem(topology)

should work.

I've also added a convenience function to go from a Residue to an OEMol:

molecule = generateOEMolFromTopologyResidue(residue)

though it only works with residues that are separate molecules and not polymeric molecules connected to other residues.

generateResidueTemplate() also supports an optional argument residue_atoms where you can provide a list of atoms that form the polymeric repeating unit of the residue. This would allow us to pass an OEMol molecule for, say, a modified amino acid with capping groups, along with a list of the atoms that are not capping groups, and it will parameterize the residue and make sure the charges are integral. This feature hasn't been tested yet, but it should be relatively easy to get working when we need it.

[jchodera]

OK, this feature passes tests locally. Merging to run full test suite.