Charge generation is failing

[nividic]

Hi there,

I was trying to use generate an openmm system starting from a given molecule topology. I'm using the openforcefield v 0.3.0 on Mac Osx. Unfortunately, I was unsuccessful. The problem is related to the charge assignment. Here is the snippet code to reproduce the issue:

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

oemol = oechem.OEMol()

with oechem.oemolistream("2799.oeb") as ifs:
    oechem.OEReadMolecule(ifs, oemol)

# This is a temporary bug fix for the bug #344 
oechem.OE3DToAtomStereo(oemol)
oechem.OE3DToBondStereo(oemol)

molecule = Molecule(oemol)

forcefield = ForceField('test_forcefields/smirnoff99Frosst.offxml')

topology = Topology.from_molecules(molecule)
system = forcefield.create_openmm_system(topology)

I'm getting the following error:

Warning: SelectElfPop: issue with removing trans COOH conformers
Warning: ElfPrepareMol: SelectElfPop failed
Warning: OEELFCharges: ELF preparation failed on mol DrugBank_2799.
Traceback (most recent call last):
File "conf.py", line 19, in
system = forcefield.create_openmm_system(topology)
File "/Users/gcalabro/local/miniconda3/envs/smirnoff_new/lib/python3.6/site-packages/openforcefield/typing/engines/smirnoff/forcefield.py", line 997, in create_openmm_system
parameter_handler.create_force(system, topology, **kwargs)
File "/Users/gcalabro/local/miniconda3/envs/smirnoff_new/lib/python3.6/site-packages/openforcefield/typing/engines/smirnoff/parameters.py", line 2255, in create_force
temp_mol.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry)
File "/Users/gcalabro/local/miniconda3/envs/smirnoff_new/lib/python3.6/site-packages/openforcefield/topology/molecule.py", line 1991, in compute_partial_charges_am1bcc
self
File "/Users/gcalabro/local/miniconda3/envs/smirnoff_new/lib/python3.6/site-packages/openforcefield/utils/toolkits.py", line 2964, in call
return method(*args, **kwargs)
File "/Users/gcalabro/local/miniconda3/envs/smirnoff_new/lib/python3.6/site-packages/openforcefield/utils/toolkits.py", line 1237, in compute_partial_charges_am1bcc
raise Exception('Unable to assign charges')
Exception: Unable to assign charges

I've attached the offending zipped molecule:

2799.oeb.zip

[davidlmobley]

Naively this sounds like an OpenEye issue. Have you checked if you can use OpenEye's AM1BCC charging engine to charge this molecule? In other words, is this really an OpenFF issue or an OpenEye issue?

[nividic]

@davidlmobley This is probably an OpenEye problem

[jchodera]

I think this might be an openforcefield bug. If you run the pure OpenEye code, it works:

from openeye import oechem, oequacpac

oemol = oechem.OEMol()

with oechem.oemolistream("2799.oeb") as ifs:
    oechem.OEReadMolecule(ifs, oemol)

# This is a temporary bug fix for the bug #344 
oechem.OE3DToAtomStereo(oemol)
oechem.OE3DToBondStereo(oemol)

# Assign ELF10 charges
oequacpac.OEAssignCharges(oemol, oequacpac.OEAM1BCCELF10Charges())

(@nividic : Note that the need for explicit stereochemistry is a feature, not a bug! If you don't know what molecule you want to parameterize, we can't guess it!)

But if you do this in OFF, we are unable to assign charges:

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

molecule = Molecule(oemol)

forcefield = ForceField('test_forcefields/smirnoff99Frosst.offxml')

topology = Topology.from_molecules(molecule)
system = forcefield.create_openmm_system(topology)

resulting in the error:

Exception                                 Traceback (most recent call last)
<ipython-input-2-7253629acf33> in <module>
      8 
      9 topology = Topology.from_molecules(molecule)
---> 10 system = forcefield.create_openmm_system(topology)

~/miniconda/lib/python3.7/site-packages/openforcefield/typing/engines/smirnoff/forcefield.py in create_openmm_system(self, topology, **kwargs)
   1126         # Add forces and parameters to the System
   1127         for parameter_handler in parameter_handlers:
-> 1128             parameter_handler.create_force(system, topology, **kwargs)
   1129 
   1130         # Let force Handlers do postprocessing

~/miniconda/lib/python3.7/site-packages/openforcefield/typing/engines/smirnoff/parameters.py in create_force(self, system, topology, **kwargs)
   2445                 #temp_mol.compute_partial_charges(quantum_chemical_method=self._quantum_chemical_method,
   2446                 #                                 partial_charge_method=self._partial_charge_method)
-> 2447                 temp_mol.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry)
   2448 
   2449             # Assign charges to relevant atoms

~/miniconda/lib/python3.7/site-packages/openforcefield/topology/molecule.py in compute_partial_charges_am1bcc(self, toolkit_registry)
   1989             charges = toolkit_registry.call(
   1990                       'compute_partial_charges_am1bcc',
-> 1991                       self
   1992             )
   1993         elif isinstance(toolkit_registry, ToolkitWrapper):

~/miniconda/lib/python3.7/site-packages/openforcefield/utils/toolkits.py in call(self, method_name, *args, **kwargs)
   2962                 method = getattr(toolkit, method_name)
   2963                 try:
-> 2964                     return method(*args, **kwargs)
   2965                 except NotImplementedError:
   2966                     pass

~/miniconda/lib/python3.7/site-packages/openforcefield/utils/toolkits.py in compute_partial_charges_am1bcc(self, molecule)
   1235 
   1236         if result is False:
-> 1237             raise Exception('Unable to assign charges')
   1238 
   1239         # Extract and return charges

Exception: Unable to assign charges

[jchodera]

It appears neither the OpenEyeToolkitWrapper

from openforcefield.utils.toolkits import ToolkitRegistry, OpenEyeToolkitWrapper, RDKitToolkitWrapper, AmberToolsToolkitWrapper
    
toolkit_registry = ToolkitRegistry()
for toolkit in [OpenEyeToolkitWrapper,]:
    toolkit_registry.register_toolkit(toolkit)
 
system = forcefield.create_openmm_system(topology, toolkit_registry=toolkit_registry)

nor the AmberToolsToolkitWrapper

from openforcefield.utils.toolkits import ToolkitRegistry, OpenEyeToolkitWrapper, RDKitToolkitWrapper, AmberToolsToolkitWrapper
    
toolkit_registry = ToolkitRegistry()
for toolkit in [AmberToolsToolkitWrapper, RDKitToolkitWrapper]:
    toolkit_registry.register_toolkit(toolkit)
 
system = forcefield.create_openmm_system(topology, toolkit_registry=toolkit_registry)

works, both failing with identical error messages that give no information about the cause:

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-10-c01906baa177> in <module>
      5     toolkit_registry.register_toolkit(toolkit)
      6 
----> 7 system = forcefield.create_openmm_system(topology, toolkit_registry=toolkit_registry)

~/miniconda/lib/python3.7/site-packages/openforcefield/typing/engines/smirnoff/forcefield.py in create_openmm_system(self, topology, **kwargs)
   1126         # Add forces and parameters to the System
   1127         for parameter_handler in parameter_handlers:
-> 1128             parameter_handler.create_force(system, topology, **kwargs)
   1129 
   1130         # Let force Handlers do postprocessing

~/miniconda/lib/python3.7/site-packages/openforcefield/typing/engines/smirnoff/parameters.py in create_force(self, system, topology, **kwargs)
   2445                 #temp_mol.compute_partial_charges(quantum_chemical_method=self._quantum_chemical_method,
   2446                 #                                 partial_charge_method=self._partial_charge_method)
-> 2447                 temp_mol.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry)
   2448 
   2449             # Assign charges to relevant atoms

~/miniconda/lib/python3.7/site-packages/openforcefield/topology/molecule.py in compute_partial_charges_am1bcc(self, toolkit_registry)
   1989             charges = toolkit_registry.call(
   1990                       'compute_partial_charges_am1bcc',
-> 1991                       self
   1992             )
   1993         elif isinstance(toolkit_registry, ToolkitWrapper):

~/miniconda/lib/python3.7/site-packages/openforcefield/utils/toolkits.py in call(self, method_name, *args, **kwargs)
   2962                 method = getattr(toolkit, method_name)
   2963                 try:
-> 2964                     return method(*args, **kwargs)
   2965                 except NotImplementedError:
   2966                     pass

~/miniconda/lib/python3.7/site-packages/openforcefield/utils/toolkits.py in compute_partial_charges_am1bcc(self, molecule)
   1235 
   1236         if result is False:
-> 1237             raise Exception('Unable to assign charges')
   1238 
   1239         # Extract and return charges

Exception: Unable to assign charges

from openforcefield.utils.toolkits import ToolkitRegistry, OpenEyeToolkitWrapper, RDKitToolkitWrapper, AmberToolsToolkitWrapper
    
toolkit_registry = ToolkitRegistry()
for toolkit in [AmberToolsToolkitWrapper, RDKitToolkitWrapper]:
    toolkit_registry.register_toolkit(toolkit)
 
system = forcefield.create_openmm_system(topology, toolkit_registry=toolkit_registry)
from openforcefield.utils.toolkits import ToolkitRegistry, OpenEyeToolkitWrapper, RDKitToolkitWrapper, AmberToolsToolkitWrapper
    
toolkit_registry = ToolkitRegistry()
for toolkit in [AmberToolsToolkitWrapper, RDKitToolkitWrapper]:
    toolkit_registry.register_toolkit(toolkit)
 
system = forcefield.create_openmm_system(topology, toolkit_registry=toolkit_registry)
---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-11-6dd9f1491efc> in <module>
      5     toolkit_registry.register_toolkit(toolkit)
      6 
----> 7 system = forcefield.create_openmm_system(topology, toolkit_registry=toolkit_registry)

~/miniconda/lib/python3.7/site-packages/openforcefield/typing/engines/smirnoff/forcefield.py in create_openmm_system(self, topology, **kwargs)
   1126         # Add forces and parameters to the System
   1127         for parameter_handler in parameter_handlers:
-> 1128             parameter_handler.create_force(system, topology, **kwargs)
   1129 
   1130         # Let force Handlers do postprocessing

~/miniconda/lib/python3.7/site-packages/openforcefield/typing/engines/smirnoff/parameters.py in create_force(self, system, topology, **kwargs)
   2445                 #temp_mol.compute_partial_charges(quantum_chemical_method=self._quantum_chemical_method,
   2446                 #                                 partial_charge_method=self._partial_charge_method)
-> 2447                 temp_mol.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry)
   2448 
   2449             # Assign charges to relevant atoms

~/miniconda/lib/python3.7/site-packages/openforcefield/topology/molecule.py in compute_partial_charges_am1bcc(self, toolkit_registry)
   1989             charges = toolkit_registry.call(
   1990                       'compute_partial_charges_am1bcc',
-> 1991                       self
   1992             )
   1993         elif isinstance(toolkit_registry, ToolkitWrapper):

~/miniconda/lib/python3.7/site-packages/openforcefield/utils/toolkits.py in call(self, method_name, *args, **kwargs)
   2962                 method = getattr(toolkit, method_name)
   2963                 try:
-> 2964                     return method(*args, **kwargs)
   2965                 except NotImplementedError:
   2966                     pass

~/miniconda/lib/python3.7/site-packages/openforcefield/utils/toolkits.py in compute_partial_charges_am1bcc(self, molecule)
   1235 
   1236         if result is False:
-> 1237             raise Exception('Unable to assign charges')
   1238 
   1239         # Extract and return charges

Exception: Unable to assign charges

[j-wags]

@jchodera, the current default constructor for ToolkitRegistry will attempt to initialize OE, Amber, and RDK toolkits. The command to initialize an empty ToolkitRegistry is ToolkitRegistry(toolkit_precedence=[]). Revisiting this months later, this is really unintuitive, and I think we should change it (ToolkitRegistry() should initialize an empty registry). Opening a new issue.

[j-wags]

The OE-backend failure reproduces for me when I generate multiple conformers before the AM1BCCELF10 calculation [1]. The molecule successfully gets charges assigned by RDK/Ambertools [2].

It's not actually clear that the calculation fails per se. But the command result = oequacpac.OEAssignCharges(oemol, oequacpac.OEAM1BCCELF10Charges()) raises some warnings, and consequently has the calculation return False, which we interpret as failure.

One possible solution would be to have some sort of fallback-behavior in this case. Maybe if AM1BCCELF10 fails, we try again, but with vanilla AM1BCC.

[1]

from openforcefield.utils.toolkits import RDKitToolkitWrapper, AmberToolsToolkitWrapper, OpenEyeToolkitWrapper, ToolkitRegistry
from openforcefield.topology import Molecule
from openforcefield.topology import Topology
from openforcefield.typing.engines.smirnoff import ForceField
from openeye import oechem

tr = ToolkitRegistry(toolkit_precedence=[OpenEyeToolkitWrapper])

oemol = oechem.OEMol()

with oechem.oemolistream("2799.oeb") as ifs:
    oechem.OEReadMolecule(ifs, oemol)

# This is a temporary bug fix for the bug #344 
oechem.OE3DToAtomStereo(oemol)
oechem.OE3DToBondStereo(oemol)

molecule = Molecule(oemol)

molecule.generate_conformers(n_conformers=10, toolkit_registry=tr)

molecule.compute_partial_charges_am1bcc(toolkit_registry=tr)

Warning: Failed to generate conformers for molecule DrugBank_2799
Warning: SelectElfPop: issue with removing trans COOH conformers
Warning: ElfPrepareMol: SelectElfPop failed
...
~/projects/OpenForceField/openforcefield/openforcefield/utils/toolkits.py in compute_partial_charges_am1bcc(self, molecule)
1244
1245 if result is False:
-> 1246 raise Exception('Unable to assign charges')
1247
1248 # Extract and return charges

Exception: Unable to assign charges

[2]

Same as [1], but replace tr=... line with the following:

tr = ToolkitRegistry(toolkit_precedence=[AmberToolsToolkitWrapper, RDKitToolkitWrapper])

[jchodera]

Warning: Failed to generate conformers for molecule DrugBank_2799
Warning: SelectElfPop: issue with removing trans COOH conformers
Warning: ElfPrepareMol: SelectElfPop failed

@j-wags : Can you come up with the minimal test case (not involving OFF tools) to reproduce this warning? It seems useful for @nvidic to relay those issues internally.

[j-wags]

@jchodera Yes -- A minimal reproducing case using only OE tools is:

from openeye import oechem
from openeye import oequacpac
from openeye import oeomega

oemol = oechem.OEMol()

with oechem.oemolistream("2799.oeb") as ifs:
    oechem.OEReadMolecule(ifs, oemol)
    
omega = oeomega.OEOmega()
omega.SetSampleHydrogens(True)
omega_status = omega(oemol)
quacpac_status = oequacpac.OEAssignCharges(oemol, oequacpac.OEAM1BCCELF10Charges())

If SetSampleHydrogens is False, then we don't get this problem. It may be the case that we shouldn't expect quacpac to handle variable protonation int he input confs. I'd be interested in your feedback on this, @nividic.

[nividic]

@j-wags @jchodera Let me talk about this with someone here at OpenEye. I'll be back to you asap

[cbayly13]

@j-wags @jchodera @nividic @davidlmobley We know what the problem is and a fix is coming. The problem is ELF10Charges (or any AM1-BCC charge set) must eliminate trans-COOH configurations from charging, and when it does that in this case there are no conformers left. It is fishy that this happens when we are asking for more sampling and it doesn't happen without... we will look into it and fix a bug if there is one (likely in my view). In the meantime, the fix James Haigh is working on will check the error message from the failed oequacpac.OEAM1BCCELF10Charges() call, and if it is this one it will re-run omega with SetSampleHydrogens set to False. Crude but effective.