`compute_partial_charges_am1bcc` seems Unable to assign charges for some ligands

Question

`compute_partial_charges_am1bcc` seems Unable to assign charges for some ligands

hannahbrucemacdonald opened this issue 5 years ago · 12 comments

hannahbrucemacdonald commented 5 years ago

Describe the bug
I am seeing Exception: Unable to assign charges when trying to parameterise some small molecules. This seems similar to #346, however the molecules fail for both oequacpac.OEAM1BCCELF10Charges and OEAM1BCCCharges

To Reproduce
python minimal_example.py
Of this set, there are 4 molecules that are failing (problematic_ids)

min_example.zip

Output

Traceback (most recent call last):
File "minimial_example.py", line 10, in
molecules[i].compute_partial_charges_am1bcc()
File "/home/brucemah/miniconda3/envs/omm_source/lib/python3.6/site-packages/openforcefield/topology/molecule.py", line 2014, in compute_partial_charges_am1bcc
self
File "/home/brucemah/miniconda3/envs/omm_source/lib/python3.6/site-packages/openforcefield/utils/toolkits.py", line 3030, in call
return method(*args, **kwargs)
File "/home/brucemah/miniconda3/envs/omm_source/lib/python3.6/site-packages/openforcefield/utils/toolkits.py", line 1278, in compute_partial_charges_am1bcc
raise Exception('Unable to assign charges')
Exception: Unable to assign charges

Computing environment (please complete the following information):

Operating system
Output of running conda list

Additional context

Answer 1 · 2020-01-21T19:55:33.000Z

@hannahbrucemacdonald Thanks for the report and example. I'm able to reproduce this error on my machine. I'll keep you updated as I work on this.

Answer 2 · 2020-01-22T04:40:57.000Z

@j-wags : Let's mark this as "high priority" since it's blocking us from finishing the free energy benchmarks of openff-1.0.0.

Answer 3 · 2020-01-22T22:21:32.000Z

@j-wags thanks for looking into this!

I think @jchodera figured it out - the problem goes away if molecule.generate_conformers()

He's going to open a PR into https://github.com/openmm/openmm-forcefields that should fix this.

Answer 4 · 2020-01-22T22:46:41.000Z

The bigger issue here is that some molecules in that tarball succeed when calling just

molecule.compute_partial_charges_am1bcc()

but others fail, and succeed only if calling

molecule.generate_conformers(n_conformers=10)
molecule.compute_partial_charges_am1bcc()

I'd love to figure out what's going on with this, and which one of these idioms is the "right" way to assign charges.

Answer 5 · 2020-01-23T00:09:35.000Z

Hm, I had figured this had something to do with the molecules' connectivity, and not the geometry. That's a bit surprising, and I'll keep looking into this issue.

@hannahbrucemacdonald Normally SMIRNOFF calculates partial charges based on several conformers it generates behind-the-scenes for a molecule. This makes the output of the charge calculation deterministic (since OE omega will generate the same conformer for the same graph representation of a molecule). Is your workflow intending to look specifically at the charges applied to the input conformers in the attached file, or is it OK to regenerate the conformers (our normal workflow)?

Regardless, I'm going to debug the quacpac failure. I'd just like to figure out if generating new conformers can get you running again in the short-term.

Answer 6 · 2020-01-23T00:13:35.000Z

I think I was under the impression that Molecule.compute_partial_charges_am1bcc() implemented our AM1-BCC ELF10 best practices, but it appears that is not the case!

The key issue was that openmm-forcefields just called compute_partial_charges_am1bcc() for GAFF, but uses the openforcefield.typing.engines.smirnoff.ForceField.create_openmm_system for SMIRNOFF (which generates conformers first).

It sounds like compute_partial_charges_am1bcc() may just use the existing conformer(s) defined in the molecule, so if a bad geometry is defined, that appears to cause quacpac to fail.

We should figure out what we actually want compute_partial_charges_am1bcc() to do---whether it is our best-practices charging method, or just a "charge this conformer" method.

Answer 7 · 2020-01-23T00:25:03.000Z

I think I was under the impression that Molecule.compute_partial_charges_am1bcc() implemented our AM1-BCC ELF10 best practices, but it appears that is not the case!

Ah, that makes sense. You're right -- It's not clear from the name whether compute_partial_charges_am1bcc() implements the whole workflow under the hood or not.

When I initially put it together, I had the same question. I decided that it's best to keep the toolkit wrappers as modular as possible, so generate_conformers and compute_partial_charges became separate functions. I'd be OK to put in a "whole shebang" method for charge calculation, with a name that reflects its function -- and it could even mirror the ChargeIncrementModel SMIRNOFF spec section -- so maybe generate_charges(partial_charge_method, n_conformers)?

Answer 8 · 2020-01-23T02:23:45.000Z

Ok, I've confirmed that it's a problem with OpenEye Quacpac rejecting the conformers due to a short interatomic distance. I think we can consider this issue to be resolved if you agree, @hannahbrucemacdonald.

I'm opening a new issue to make OpenEyeToolkitWrapper failures more verbose now that I know how to wrangle OE output streams -- having the distance error print out initially would be much friendlier to users than the current opaque "Unable to generate charges" message.

from openforcefield.topology import Molecule
from openeye import oechem, oeomega, oequacpac
molecules = Molecule.from_file('/Users/jwagner/Downloads/PTP1B_ligands_shifted.sdf')

problematic_ids = [3,7,14,21]

for idx in range(len(molecules)):
    mol = molecules[idx]
    print(idx, mol.total_charge, mol.to_smiles())
    oemol = mol.to_openeye()
    copymol = oechem.OEMol(oemol)

    errfs = oechem.oeosstream()
    oechem.OEThrow.SetOutputStream(errfs)
    oechem.OEThrow.Clear()
    quacpac_status = oequacpac.OEAssignCharges(oemol, oequacpac.OEAM1BCCELF10Charges())
    oechem.OEThrow.SetOutputStream(oechem.oeerr)#restoring to original state
    print(errfs.str().decode("UTF-8"))

0 -2 [H]c1c(c(c(c(c1[H])[H])C(=O)N([H])c2c(c(c(c(c2[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H])[H])[H])[H])[H]

1 -2 [H]c1c(c(c(c(c1[H])N([H])S(=O)(=O)c2c(c(c(c(c2[H])[H])F)[H])[H])[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

2 -2 [H]c1c(c(c(c(c1[H])N([H])C([H])([H])C2(C(C(C(C(C2([H])[H])([H])[H])([H])[H])([H])[H])([H])[H])[H])[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

3 -2 [H]c1c(c(c(c(c1[H])N([H])C2(C(C(N(C(C2([H])[H])([H])[H])C(=O)OC([H])([H])[H])([H])[H])([H])[H])[H])[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]
Warning: Distance between atom 29 and 38 is below 0.7A
Please, correct the input
Warning: OESEmpPartialCharges: Stopping the charging process for molecule 23480
Warning: OESEmpBCCPartialCharges: Stopping the charging process at the second stage for mol 23480

4 -2 [H]c1c(c(c(c(c1[H])N([H])C(=O)N([H])C([H])(C([H])([H])[H])C([H])([H])[H])[H])C2=C(C(=C(S2)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

5 -2 [H]c1c(c(c(c(c1[H])[H])C([H])([H])S(=O)(=O)N2C(C(C(C(C2([H])[H])([H])[H])([H])N([H])c3c(c(c(c(c3[H])C4=C(C(=C(S4)C(=O)[O-])OC([H])([H])C(=O)[O-])Cl)[H])[H])[H])([H])[H])([H])[H])[H])[H]

6 -2 [H]c1c(c(c(c(c1[H])[H])C([H])([H])N([H])c2c(c(c(c(c2[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H])[H])[H])[H])[H]

7 -2 [H]c1c(c(c(c(c1[H])N([H])C2(C(C(N(C(C2([H])[H])([H])[H])S(=O)(=O)C([H])([H])[H])([H])[H])([H])[H])[H])[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]
Warning: Distance between atom 29 and 38 is below 0.7A
Please, correct the input
Warning: OESEmpPartialCharges: Stopping the charging process for molecule 23482
Warning: OESEmpBCCPartialCharges: Stopping the charging process at the second stage for mol 23482

8 -2 [H]c1c(c(c(c(c1[H])[H])S(=O)(=O)N2C(C(C(C(C2([H])[H])([H])[H])([H])N([H])c3c(c(c(c(c3[H])C4=C(C(=C(S4)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H])[H])[H])([H])[H])([H])[H])[H])[H]

9 -2 [H]c1c(c(c(c(c1[H])[H])C([H])([H])Oc2c(c(c(c(c2[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H])[H])[H])[H])[H]

10 -2 [H]c1c(c(c(c(c1[H])N([H])C2(C(C(C(C(C2([H])[H])([H])[H])([H])[H])([H])[H])([H])[H])[H])[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

11 -2 [H]c1c(c(c(c(c1[H])N([H])C(=O)C([H])([H])C([H])([H])[H])[H])C2=C(C(=C(S2)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

12 -2 [H]c1c(c(c(c(c1[H])N([H])C2(C(C(C(C2([H])[H])([H])[H])([H])[H])([H])[H])[H])[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

13 -2 [H]c1c(c(c(c(c1[H])[H])C([H])([H])S(=O)(=O)N2C(C(C(C(C2([H])[H])([H])[H])([H])N([H])c3c(c(c(c(c3[H])C4=C(C(=C(S4)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H])[H])[H])([H])[H])([H])[H])[H])[H]

14 -2 [H]c1c(c(c(c(c1[H])N([H])C2(C(C(OC(C2([H])[H])([H])[H])([H])[H])([H])[H])[H])[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]
Warning: Distance between atom 16 and 33 is below 0.7A
Please, correct the input
Warning: OESEmpPartialCharges: Stopping the charging process for molecule 23477
Warning: OESEmpBCCPartialCharges: Stopping the charging process at the second stage for mol 23477

15 -2 [H]c1c(c(c(c(c1[H])N([H])C(=O)OC([H])([H])[H])[H])C2=C(C(=C(S2)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

16 -2 [H]c1c(c(c(c(c1[H])[H])C([H])([H])S(=O)(=O)N2C(C(C(C(C2([H])[H])([H])[H])([H])N([H])c3c(c(c(c(c3[H])C4=C(C(=C(S4)C(=O)[O-])OC([H])([H])C(=O)[O-])C([H])([H])[H])[H])[H])[H])([H])[H])([H])[H])[H])[H]

17 -2 [H]c1c(c(c(c(c1[H])OC([H])([H])[H])[H])C2=C(C(=C(S2)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

18 -2 [H]c1c(c(c(c(c1[H])N([H])[H])[H])C2=C(C(=C(S2)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

19 -2 [H]c1c(c(c(c(c1[H])N([H])C2(C(C(C(C(C(C2([H])[H])([H])[H])([H])[H])([H])[H])([H])[H])([H])[H])[H])[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

20 -2 [H]c1c(c(c(c(c1[H])[H])C([H])([H])S(=O)(=O)N2C(C(C(C(C2([H])[H])([H])[H])([H])N([H])c3c(c(c(c(c3[H])C4=C(C(=C(S4)C(=O)[O-])OC([H])([H])C(=O)[O-])[H])[H])[H])[H])([H])[H])([H])[H])[H])[H]

21 -2 [H]c1c(c(c(c(c1[H])N([H])C2(C(C(N(C(C2([H])[H])([H])[H])C(=O)C([H])([H])[H])([H])[H])([H])[H])[H])[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]
Warning: Distance between atom 28 and 37 is below 0.7A
Please, correct the input
Warning: OESEmpPartialCharges: Stopping the charging process for molecule 23479
Warning: OESEmpBCCPartialCharges: Stopping the charging process at the second stage for mol 23479

22 -2 [H]c1c(c(c(c(c1[H])N([H])C2(C(C(C(C(C2([H])[H])(C([H])([H])[H])C([H])([H])[H])([H])[H])(C([H])([H])[H])C([H])([H])[H])([H])[H])[H])[H])C3=C(C(=C(S3)C(=O)[O-])OC([H])([H])C(=O)[O-])Br)[H]

Answer 9 · 2020-01-23T04:46:38.000Z

@j-wags : It would also be great to discuss what the API should do. I don't think it's good design to have molecule.compute_partial_charges_am1bcc() either

use a conformer if it's defined
use some of multiple conformers (how many?) if they're defined
fail if it doesn't like a conformer
fail if no conformers are defined

and to require that, to achieve the established best practice for assigning charges, we have to overwrite the conformers in the molecule with

molecule.generate_conformers(n_conformers=10)
molecule.compute_partial_charges_am1bcc()

This is just bad design, because it (1) doesn't produce a well-determined result, (2) violates the principle of minimum surprise in causing (or using) unintended side-effects, and (3) is verbose and confusing compared to a single API point that just does the best practice.

Answer 10 · 2020-01-25T04:18:44.000Z

@j-wags: We've resolved the urgency in openmm/openmmforcefields#95, but it would still be useful to revisit what API design we really want here.

Answer 11 · 2020-02-04T00:28:07.000Z

One other consideration here: we want to be sure we are steering people toward clear, unambiguous behavior to make it very easy to use current best practices, while still keeping optional flexibility so we can experiment as needed (with a bit more effort to select non-default behavior).

Answer 12 · 2020-02-13T19:15:42.000Z

Per discussion between @jchodera and myself on Monday --

Big picture is that we'll do best practices, as defined by:

compute_partial_charges_am1bcc, assign_partial_charges, and assign_fractional_bond_orders will disregard existing conformers unless an override is provided.
create_openmm_system will disregard existing partial charges and fractional bond orders unless an override is provided
We define best practice for the am1bccelf10 method as "generate a lot of conformers, potentially by using a low RMS cutoff, since ~~we need to start with at least 50 to get multiple conformers due to the lowest-2% step in the ELF10 process"~~ (I talked to Chris Bayly and, for molecules with fewer than 50 conformers, the fraction considered is >2%)

This issue will be closed when the following changes are made:

forcefield.create_openmm_system
- new kwarg: fractional_bond_orders_from_molecules: iterable of Molecule
molecule.generate_conformers
- new kwarg: rms_cutoff --> float or Quantity (distance dimension). If float, assume Angstrom
Molecule/OpenEyeToolkitWrapper/AmberToolsToolkitWrapper.compute_partial_charges_am1bcc
- Raise DeprecationWarning -- Tell users to use new assign_partial_charges(partial_charge_method="am1bcc") in the future
- behavior change: disregard existing confs by default
- new kwarg: use_conformers = conformer or [iterable_of_conformers]
- Update releasenotes: mention behavior change
new methods: Molecule/OpenEyeToolkitWrapper/AmberToolsToolkitWrapper.assign_partial_charges
- behavior: disregard existing confs by default. Make new conformers with a small RMS cutoff so ELF10 has several to choose from
- kwarg: partial_charge_method --> str. Name of a partial charge method that is validated by underlying toolkit. ValueError raised if the method can't be handled.
- kwarg: use_conformers --> conformer or [iterable_of_conformers]
- kwarg: strict_n_conformers --> bool, saying whether to raise a ValueError if the the number of input conformers is incorrect (versus a warning). This provides a pathway for an error if someone produces a logically inconsistent OFFXML with, for example <ChargeIncrementModel method="graph-based-method" num_conformers="5">

Name Version Build Channel