Errors with get charges in generateResidueTemplete

[gregoryross]

I'm trying to setup a system in openmm that contains only (non-peptide) small molecules and am running into problems with generating the parameters, particular with using Modeller to add solvent.

My code snippet is the following:

pdb = app.PDBFile('nanoparticle.pdb')
forcefield = app.ForceField('gaff.xml','tip3p.xml')
forcefield.registerTemplateGenerator(gaffTemplateGenerator)
modeller = app.Modeller(pdb.topology, pdb.positions)
modeller.addSolvent(forcefield)
system = forcefield.createSystem(modeller.topology,nonbondedMethod=app.PME, nonbondedCutoff =1.0*unit.nanometer, constraints=app.HBonds)

Running the last command throws up the error below. Any suggestions? I've attached the input PDB file as a .txt file so github can load it. Thanks!

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-8-1c51bb5ba4ee> in <module>()
----> 1 system = forcefield.createSystem(modeller.topology,nonbondedMethod=app.PME, nonbondedCutoff=1.0*unit.nanometer, constraints=app.HBonds)

/Users/rossg/miniconda2/lib/python2.7/site-packages/simtk/openmm/app/forcefield.pyc in createSystem(self, topology, nonbondedMethod, nonbondedCutoff, constraints, rigidWater, removeCMMotion, hydrogenMass, **args)
    771                     # No existing templates match.  Try any registered residue template generators.
    772                     for generator in self._templateGenerators:
--> 773                         if generator(self, res):
    774                             # This generator has registered a new residue template that should match.
    775                             [template, matches] = self._getResidueTemplateMatches(res, bondedToAtom)

/Users/rossg/miniconda2/lib/python2.7/site-packages/openmoltools/forcefield_generators.pyc in gaffTemplateGenerator(forcefield, residue, structure)
    483 
    484     # Generate template and parameters.
--> 485     [template, ffxml] = generateResidueTemplate(molecule)
    486 
    487     # Register the template.

/Users/rossg/miniconda2/lib/python2.7/site-packages/openmoltools/forcefield_generators.pyc in generateResidueTemplate(molecule, residue_atoms)
    260 
    261     # Generate canonical AM1-BCC charges and a reference conformation.
--> 262     molecule = get_charges(molecule, strictStereo=False, keep_confs=1)
    263 
    264     # DEBUG: This may be necessary.

/Users/rossg/miniconda2/lib/python2.7/site-packages/openmoltools/openeye.pyc in get_charges(molecule, max_confs, strictStereo, normalize, keep_confs)
     66 
     67     if not status:
---> 68         raise(RuntimeError("OEAssignPartialCharges returned error code %d" % status))
     69 
     70     #Determine conformations to return

RuntimeError: OEAssignPartialCharges returned error code 0

nanoparticle.txt

[jchodera]

Looks like charging failed. Usually this is because the molecule was misperceived because the topology was missing some bonds, but it looks like your PDB has the appropriate CONECT records. I'll see if I can have this dump a mol2 file when it raises an exception so that we can more easily debug what is going on.

[gregoryross]

Thanks. The error does not occur when modeller isn't called. For instance, the snippet below doesn't throw up any errors.

forcefield = app.ForceField('gaff.xml','tip3p.xml')
forcefield.registerTemplateGenerator(gaffTemplateGenerator)
pdb = app.PDBFile("nanoparticle.pdb")
system = forcefield.createSystem(pdb.topology,nonbondedMethod=app.PME, nonbondedCutoff=1.0*unit.nanometer, constraints=app.HBonds)

[jchodera]

Working on this now!

[jchodera]

I got some more useful debug info out:

mski1776:openmoltools.choderalab choderaj$ python test.py
Warning: SelectElfDiverseConfs: elfPop.NumConfs 1 <= elfLimit 1
Warning: BCC model does not support molecules containing Na
Warning:  Getting BCCs for sodium(1+) was unsuccessful; stopping the charging process for this molecule
Warning:  Getting BCCs for sodium(1+) was unsuccessful; stopping the charging process for this molecule
Warning: Formal charge(1) is not equal to sum of partial charges(0.00)
     for sodium(1+).
Traceback (most recent call last):
  File "test.py", line 10, in <module>
    system = forcefield.createSystem(modeller.topology,nonbondedMethod=app.PME, nonbondedCutoff =1.0*unit.nanometer, constraints=app.HBonds)
  File "/Users/choderaj/miniconda/lib/python2.7/site-packages/simtk/openmm/app/forcefield.py", line 773, in createSystem
    if generator(self, res):
  File "/Users/choderaj/github/choderalab/openmoltools/openmoltools.choderalab/openmoltools/forcefield_generators.py", line 491, in gaffTemplateGenerator
    [template, ffxml] = generateResidueTemplate(molecule)
  File "/Users/choderaj/github/choderalab/openmoltools/openmoltools.choderalab/openmoltools/forcefield_generators.py", line 263, in generateResidueTemplate
    molecule = get_charges(molecule, strictStereo=False, keep_confs=1)
  File "/Users/choderaj/github/choderalab/openmoltools/openmoltools.choderalab/openmoltools/openeye.py", line 68, in get_charges
    raise(RuntimeError("OEAssignPartialCharges returned error code %d" % status))
RuntimeError: OEAssignPartialCharges returned error code 0

Looks like the problem is sodium. Can you include an xml file with ions in the app.ForceField call?

[jchodera]

I can confirm that changing

forcefield = app.ForceField('gaff.xml','tip3p.xml')

to

forcefield = app.ForceField('gaff.xml','tip3p.xml','amber99sbildn.xml')

eliminates the error because Na parameters are provided in amber99sbildn.xml.

[gregoryross]

Great, thanks for looking into this!