`compute_partial_charges_am1bcc` seems Unable to assign charges for some ligands #492

hannahbrucemacdonald · 2020-01-21T19:43:16Z

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

The text was updated successfully, but these errors were encountered:

j-wags · 2020-01-21T19:55:33Z

@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.

jchodera · 2020-01-22T04:40:57Z

@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.

hannahbrucemacdonald · 2020-01-22T22:21:32Z

@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.

jchodera · 2020-01-22T22:46:41Z

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.

j-wags · 2020-01-23T00:09:35Z

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.

jchodera · 2020-01-23T00:13:35Z

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.

j-wags · 2020-01-23T00:25:03Z

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)?

j-wags · 2020-01-23T02:23:45Z

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]

jchodera · 2020-01-23T04:46:38Z

@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.

jchodera · 2020-01-25T04:18:44Z

@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.

jchodera · 2020-02-04T00:28:07Z

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).

j-wags · 2020-02-13T19:15:42Z

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">

j-wags added bug external dependencies labels Jan 21, 2020

jchodera mentioned this issue Jan 22, 2020

Fix GAFF AM1-BCC charging failures openmm/openmmforcefields#95

Merged

j-wags added the priority:high label Jan 23, 2020

j-wags mentioned this issue Jan 23, 2020

Make OpenEyeToolkitWrapper failures more verbose #495

Open

jchodera removed the priority:high label Jan 25, 2020

j-wags mentioned this issue Feb 13, 2020

Rerun omega without stereo requirement if initial conf gen failed #517

Merged

4 tasks

j-wags mentioned this issue Feb 21, 2020

Add ChargeIncrementModelHandler #471

Merged

14 tasks

j-wags closed this as completed in #471 Jun 19, 2020

j-wags mentioned this issue Jun 26, 2020

Error raised for create_openmm_system() (AM1-BCC charges) from molecules with net charge != 0 #627

Closed

mattwthompson mentioned this issue Nov 29, 2023

Do not generate conformers that will not be used openmm/openmmforcefields#315

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

`compute_partial_charges_am1bcc` seems Unable to assign charges for some ligands #492

`compute_partial_charges_am1bcc` seems Unable to assign charges for some ligands #492

hannahbrucemacdonald commented Jan 21, 2020

j-wags commented Jan 21, 2020

jchodera commented Jan 22, 2020

hannahbrucemacdonald commented Jan 22, 2020

jchodera commented Jan 22, 2020 •

edited

Loading

j-wags commented Jan 23, 2020

jchodera commented Jan 23, 2020

j-wags commented Jan 23, 2020

j-wags commented Jan 23, 2020 •

edited

Loading

jchodera commented Jan 23, 2020

jchodera commented Jan 25, 2020

jchodera commented Feb 4, 2020

j-wags commented Feb 13, 2020 •

edited

Loading

compute_partial_charges_am1bcc seems Unable to assign charges for some ligands #492

compute_partial_charges_am1bcc seems Unable to assign charges for some ligands #492

Comments

hannahbrucemacdonald commented Jan 21, 2020

Name Version Build Channel

j-wags commented Jan 21, 2020

jchodera commented Jan 22, 2020

hannahbrucemacdonald commented Jan 22, 2020

jchodera commented Jan 22, 2020 • edited Loading

j-wags commented Jan 23, 2020

jchodera commented Jan 23, 2020

j-wags commented Jan 23, 2020

j-wags commented Jan 23, 2020 • edited Loading

jchodera commented Jan 23, 2020

jchodera commented Jan 25, 2020

jchodera commented Feb 4, 2020

j-wags commented Feb 13, 2020 • edited Loading

`compute_partial_charges_am1bcc` seems Unable to assign charges for some ligands #492

`compute_partial_charges_am1bcc` seems Unable to assign charges for some ligands #492

jchodera commented Jan 22, 2020 •

edited

Loading

j-wags commented Jan 23, 2020 •

edited

Loading

j-wags commented Feb 13, 2020 •

edited

Loading