Looking for the right force field

[tarungog]

I need to rapidly minimize conformations of organic molecules that are anywhere from 250-500 atoms in size (lignin polymers). I was using MMFF in RDKit before, but I would like to access the GPU acceleration. I think this still fits under the GAFF small molecule regime, so I'm wondering if someone would point me to the right force field that might get me the results I want. When I used the one from the example, I got an error b.c. the system is expecting a protein.

molOrig = Chem.MolFromSmiles("CCCCCCCCCC")
# AllChem.EmbedMultipleConfs(molOrig)
# Chem.rdmolops.AssignAtomChiralTagsFromStructure(molOrig)
from openforcefield.topology import Molecule
ofmol = Molecule.from_rdkit(molOrig)
ofmol.generate_conformers(n_conformers=10)
pos = ofmol.conformers[0]

forcefield_kwargs = { 'constraints' : app.HBonds, 'rigidWater' : True, 'removeCMMotion' : False, 'hydrogenMass' : 4*unit.amu }
# Initialize a SystemGenerator using GAFF
from openmmforcefields.generators import SystemGenerator
system_generator = SystemGenerator(forcefields=['amber/ff14SB.xml', 'amber/tip3p_standard.xml'], small_molecule_forcefield='gaff-2.11', forcefield_kwargs=forcefield_kwargs, cache='db.json')
# Create an OpenMM System from an Open Force Field toolkit Topology object
system = system_generator.create_system(ofmol.to_topology().to_openmm())

Any suggestions where I could minimize such a molecule rapidly would be appreciated. It takes around 1 second in RDKit with MMFF, so any improvement on that would be great.

[jchodera]

Try the modified code:

from rdkit import Chem
molOrig = Chem.MolFromSmiles("CCCCCCCCCC")
molOrig = Chem.AddHs(molOrig)
# AllChem.EmbedMultipleConfs(molOrig)
# Chem.rdmolops.AssignAtomChiralTagsFromStructure(molOrig)
from openforcefield.topology import Molecule
ofmol = Molecule.from_rdkit(molOrig)
#generate_unique_atom_names(ofmol)
ofmol.generate_conformers(n_conformers=10)
pos = ofmol.conformers[0]
ofmol.name = 'molecule'

from simtk.openmm import app
from simtk import unit
forcefield_kwargs = { 'constraints' : app.HBonds, 'rigidWater' : True, 'removeCMMotion' : False, 'hydrogenMass' : 4*unit.amu }
# Initialize a SystemGenerator using GAFF
from openmmforcefields.generators import SystemGenerator
system_generator = SystemGenerator(small_molecule_forcefield='gaff-2.11', molecules=[ofmol], forcefield_kwargs=forcefield_kwargs, cache='db.json')
# Create an OpenMM System from an Open Force Field toolkit Topology object
system = system_generator.create_system(ofmol.to_topology().to_openmm())

There were two issues I spotted:

1. You specified an implicit hydrogen molecule, which (for some reason) was not converted into an explicit hydrogen openforcefield.topology.Molecule when converted to an openforcefield Molecule object.
2. You neglected to pass the ofmol to the molecules=[...] argument of SystemGenerator.

[tarungog]

@jchodera Thank you for your input! works like a charm for that decane. But when I try more complicated molecules, I keep coming across this error:

If I instead embed the conformer in RDKit, assign chirality, and then create the ofmol, I get this separate error:

Any ideas?

[jchodera]

You need to either specify the stereochemistry of your molecule (otherwise, we don't know which molecule to parameterize!) or use the Molecule.from_rdkit(rdmol, allow_unspecified_stereo=True) flag if you just want us to pick a random stereoisomer to parameterize!

[tarungog]

@jchodera if i provide an exact conformation from the RDKit side, the stereochemistry is fully parametrized, correct?

AllChem.EmbedMolecule(molOrig)
Chem.rdmolops.AssignAtomChiralTagsFromStructure(molOrig)

Then why does that second stack trace show up? I can't have a random stereoisomer, as I need to minimize an exact conformer provided from the RDKit side.

Edit: Yeah, in fact, after using the allow_undefined_stereo I am right back at this error:

[jchodera]

@jchodera if i provide an exact conformation from the RDKit side, the stereochemistry is fully parametrized, correct?

It should!

Yeah, in fact, after using the allow_undefined_stereo I am right back at this error:

This looks like an issue with antechamber failing to charge your molecule. If you write out the molecule directly to SDF or PDB, are you able to charge it with antechamber directly?

antechamber -fi mdl -i input.sdf -fo mol2 -o output-gaff.mol2 -c bcc

I suspect this is what is failing now.

I can check if you upload an SDF.

[tarungog]

@jchodera Uploaded my sdf to this link: https://bashupload.com/wp9l4/out.sdf

There are no formal charges on the molecule, so I'm confused what's going wrong.