-
Notifications
You must be signed in to change notification settings - Fork 22
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Handle some edge cases with preset charges #1070
base: develop
Are you sure you want to change the base?
Changes from all commits
9400e6c
68a390f
884b221
31edc75
99801b0
3aa0337
4099a1e
25b144e
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
@@ -1,5 +1,34 @@ | ||||||||||||
# Sharp edges | ||||||||||||
|
||||||||||||
## Quirks of charge assignment | ||||||||||||
|
||||||||||||
### Charge assignment hierarchy | ||||||||||||
|
||||||||||||
Interchange, following the [SMIRNOFF specification](https://openforcefield.github.io/standards/standards/smirnoff/#partial-charge-and-electrostatics-models), assigns charges to (heavy) atoms with the following priority: | ||||||||||||
|
||||||||||||
1. **Preset charges**: Look for molecule matches in the `charge_from_molecules` argument | ||||||||||||
2. **Library charges**: Look for chemical environment matches in library charges | ||||||||||||
3. **Charge increments**: Look for chemical environment matches in charge increments | ||||||||||||
4. **AM1-BCC**: Try to run some variant of AM1-BCC (presumably this is were graph charges go) | ||||||||||||
|
||||||||||||
If charges are successfully assigned using a method, no lower-priority methods in this list are attempted. For example: | ||||||||||||
|
||||||||||||
* A force field with library charge parameters for peptides (i.e. a biopolymer force field, covering all appropriate residues) will NOT try to call AM1-BCC on a biopolymers. | ||||||||||||
* If a ligand is successfully assigned preset charges, chemical environment matching of library charges and charge increments will be skipped, as will AM1-BCC. | ||||||||||||
* If a variant of AM1-BCC (i.e. using something other than AM1 and/or using custom BCCs) is encoded in a `<ChargeIncrementModel>` section, other AM1-BCC implementations will not be called. | ||||||||||||
* If preset charges are not provided, a force field like OpenFF's Parsley or Sage lines without (many) library charges or charge increments will attempt AM1-BCC on all molecules (except water and monoatomic ions). | ||||||||||||
|
||||||||||||
After all of these steps are complete and all heavy atoms given partial charges, virtual sites are assigned charges using the values of `charge_increment`s in the virtual site parameters. Because virtual site charge are only described by the force field, using preset charges with virtual sites is discouraged. | ||||||||||||
|
||||||||||||
### Preset charges | ||||||||||||
|
||||||||||||
The following restrictions are in place when using preset charges: | ||||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
Just think this could use a little extra context. |
||||||||||||
|
||||||||||||
* All molecules in the the `charge_from_molecules` list must be non-isomorphic with each other. | ||||||||||||
* All molecules in the the `charge_from_molecules` list must have partial charges. | ||||||||||||
Comment on lines
+27
to
+28
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
I think this one-to-many relationship is worth making explicit. |
||||||||||||
|
||||||||||||
Using preset charges with virtual sites is discouraged as it can provide surprising results. | ||||||||||||
|
||||||||||||
## Quirks of `from_openmm` | ||||||||||||
|
||||||||||||
### Modified masses are ignored | ||||||||||||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -7,6 +7,7 @@ | |
from openff.interchange._tests import get_test_file_path | ||
from openff.interchange.constants import kcal_mol, kcal_mol_a2 | ||
from openff.interchange.exceptions import ( | ||
PresetChargesError, | ||
UnassignedAngleError, | ||
UnassignedBondError, | ||
UnassignedTorsionError, | ||
|
@@ -21,6 +22,7 @@ | |
_upconvert_vdw_handler, | ||
library_charge_from_molecule, | ||
) | ||
from openff.interchange.warnings import PresetChargesAndVirtualSitesWarning | ||
|
||
|
||
def _get_interpolated_bond_k(bond_handler) -> float: | ||
|
@@ -160,7 +162,8 @@ def test_library_charges_from_molecule(): | |
assert library_charges.charge == [*mol.partial_charges] | ||
|
||
|
||
class TestChargeFromMolecules: | ||
class TestPresetCharges: | ||
@pytest.mark.slow | ||
def test_charge_from_molecules_basic(self, sage): | ||
molecule = Molecule.from_smiles("CCO") | ||
molecule.assign_partial_charges(partial_charge_method="am1bcc") | ||
|
@@ -231,6 +234,98 @@ def test_charges_from_molecule_reordered( | |
|
||
assert numpy.allclose(expected_charges, found_charges) | ||
|
||
def test_error_when_any_missing_partial_charges(self, sage): | ||
topology = Topology.from_molecules( | ||
[ | ||
Molecule.from_smiles("C"), | ||
Molecule.from_smiles("CCO"), | ||
], | ||
) | ||
Comment on lines
+238
to
+243
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Should this topology include a molecule with partial charges, to ensure that it's testing if there are any missing charges? |
||
|
||
with pytest.raises( | ||
PresetChargesError, | ||
match="All.*must have partial charges", | ||
): | ||
sage.create_interchange( | ||
topology, | ||
charge_from_molecules=[Molecule.from_smiles("C")], | ||
) | ||
|
||
def test_error_multiple_isomorphic_molecules(self, sage, ethanol, reversed_ethanol): | ||
# these ethanol fixtures *do* have charges, if they didn't have charges it would be | ||
# ambiguous which error should be raised | ||
topology = Topology.from_molecules([ethanol, reversed_ethanol, ethanol]) | ||
|
||
with pytest.raises( | ||
PresetChargesError, | ||
match="All molecules.*isomorphic.*unique", | ||
): | ||
sage.create_interchange( | ||
topology, | ||
charge_from_molecules=[ethanol, reversed_ethanol], | ||
) | ||
|
||
def test_virtual_site_charge_increments_applied_after_preset_charges_water( | ||
self, | ||
tip4p, | ||
water, | ||
): | ||
""" | ||
Test that charge increments to/from virtual sites are still applied after preset charges are set. | ||
|
||
This is funky user input, since preset charges override library charges, which are important for water. | ||
""" | ||
water.partial_charges = Quantity( | ||
[-2.0, 1.0, 1.0], | ||
"elementary_charge", | ||
) | ||
|
||
with pytest.warns(PresetChargesAndVirtualSitesWarning): | ||
# This dict is probably ordered O H H VS | ||
charges = tip4p.create_interchange( | ||
water.to_topology(), | ||
charge_from_molecules=[water], | ||
)["Electrostatics"].charges | ||
|
||
# tip4p_fb.offxml is meant to result in charges of -1.0517362213526 (VS) 0 (O) and 0.5258681106763 (H) | ||
# the force field has 0.0 for all library charges (skipping AM1-BCC), so preset charges screw this up. | ||
# the resulting charges, if using preset charges of [-2, 1, 1], should be the result of adding together | ||
# [-2, 1, 1, 0] (preset charges) and | ||
# [0, 1.5258681106763001, 1.5258681106763001, -1.0517362213526] (charge increments) | ||
|
||
numpy.testing.assert_allclose( | ||
[charge.m for charge in charges.values()], | ||
[-2.0, 1.5258681106763001, 1.5258681106763001, -1.0517362213526], | ||
) | ||
|
||
def test_virtual_site_charge_increments_applied_after_preset_charges_ligand( | ||
self, | ||
sage_with_sigma_hole, | ||
): | ||
ligand = Molecule.from_mapped_smiles("[C:1]([Cl:2])([Cl:3])([Cl:4])[Cl:5]") | ||
|
||
ligand.partial_charges = Quantity( | ||
[1.2, -0.3, -0.3, -0.3, -0.3], | ||
"elementary_charge", | ||
) | ||
|
||
with pytest.warns(PresetChargesAndVirtualSitesWarning): | ||
# This dict is probably ordered C Cl Cl Cl Cl VS VS VS VS | ||
charges = sage_with_sigma_hole.create_interchange( | ||
ligand.to_topology(), | ||
charge_from_molecules=[ligand], | ||
)["Electrostatics"].charges | ||
|
||
# this fake sigma hole parameter shifts 0.1 e from the carbon and 0.2 e from the chlorine, so the | ||
# resulting charges should be like | ||
# [1.2, -0.3, -0.3, -0.3, -0.3] + | ||
# [0.4, 0.2, 0.2, 0.2, 0.2, -0.3, -0.3, -0.3, -0.3] | ||
|
||
numpy.testing.assert_allclose( | ||
[charge.m for charge in charges.values()], | ||
[1.6, -0.1, -0.1, -0.1, -0.1, -0.3, -0.3, -0.3, -0.3], | ||
) | ||
|
||
|
||
@skip_if_missing("openmm") | ||
class TestPartialBondOrdersFromMolecules: | ||
|
Original file line number | Diff line number | Diff line change | ||||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
@@ -109,7 +109,9 @@ def from_smirnoff( | |||||||||||||||||||||
positions are taken from the molecules in topology, if present on all molecules. | ||||||||||||||||||||||
charge_from_molecules : `List[openff.toolkit.molecule.Molecule]`, optional | ||||||||||||||||||||||
If specified, partial charges will be taken from the given molecules | ||||||||||||||||||||||
instead of being determined by the force field. | ||||||||||||||||||||||
instead of being determined by the force field. All molecules in this list | ||||||||||||||||||||||
must have partial charges assigned and must not be isomorphic with any other | ||||||||||||||||||||||
molecules in the list. | ||||||||||||||||||||||
Comment on lines
111
to
+114
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
I think this might help clear up the misconception that I sometimes see hinted at that this only affects the same molecule in the topology - like if I have four copies of a molecule, and pass literally one of those objects to |
||||||||||||||||||||||
partial_bond_orders_from_molecules : List[openff.toolkit.molecule.Molecule], optional | ||||||||||||||||||||||
If specified, partial bond orders will be taken from the given molecules | ||||||||||||||||||||||
instead of being determined by the force field. | ||||||||||||||||||||||
|
@@ -144,7 +146,7 @@ def from_smirnoff( | |||||||||||||||||||||
topology=topology, | ||||||||||||||||||||||
box=box, | ||||||||||||||||||||||
positions=positions, | ||||||||||||||||||||||
charge_from_molecules=charge_from_molecules, | ||||||||||||||||||||||
molecules_with_preset_charges=charge_from_molecules, | ||||||||||||||||||||||
partial_bond_orders_from_molecules=partial_bond_orders_from_molecules, | ||||||||||||||||||||||
allow_nonintegral_charges=allow_nonintegral_charges, | ||||||||||||||||||||||
) | ||||||||||||||||||||||
|
Original file line number | Diff line number | Diff line change | ||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
@@ -1,3 +1,5 @@ | ||||||||||||||||
import warnings | ||||||||||||||||
|
||||||||||||||||
from openff.toolkit import ForceField, Molecule, Quantity, Topology | ||||||||||||||||
from openff.toolkit.typing.engines.smirnoff import ParameterHandler | ||||||||||||||||
from openff.toolkit.typing.engines.smirnoff.plugins import load_handler_plugins | ||||||||||||||||
|
@@ -8,6 +10,7 @@ | |||||||||||||||
from openff.interchange.components.toolkit import _check_electrostatics_handlers | ||||||||||||||||
from openff.interchange.exceptions import ( | ||||||||||||||||
MissingParameterHandlerError, | ||||||||||||||||
PresetChargesError, | ||||||||||||||||
SMIRNOFFHandlersNotImplementedError, | ||||||||||||||||
) | ||||||||||||||||
from openff.interchange.plugins import load_smirnoff_plugins | ||||||||||||||||
|
@@ -25,6 +28,7 @@ | |||||||||||||||
SMIRNOFFProperTorsionCollection, | ||||||||||||||||
) | ||||||||||||||||
from openff.interchange.smirnoff._virtual_sites import SMIRNOFFVirtualSiteCollection | ||||||||||||||||
from openff.interchange.warnings import PresetChargesAndVirtualSitesWarning | ||||||||||||||||
|
||||||||||||||||
_SUPPORTED_PARAMETER_HANDLERS: set[str] = { | ||||||||||||||||
"Constraints", | ||||||||||||||||
|
@@ -93,17 +97,61 @@ def validate_topology(value): | |||||||||||||||
) | ||||||||||||||||
|
||||||||||||||||
|
||||||||||||||||
def _preprocess_preset_charges( | ||||||||||||||||
molecules_with_preset_charges: list[Molecule] | None, | ||||||||||||||||
) -> list[Molecule] | None: | ||||||||||||||||
""" | ||||||||||||||||
Pre-process the molecules_with_preset_charges argument. | ||||||||||||||||
|
||||||||||||||||
If molecules_with_preset_charges is None, return None. | ||||||||||||||||
|
||||||||||||||||
If molecules_with_preset_charges is list[Molecule], ensure that | ||||||||||||||||
|
||||||||||||||||
1. The input is a list of Molecules | ||||||||||||||||
2. Each molecule has assign partial charges | ||||||||||||||||
3. Ensure no molecules are isomorphic with another in the list | ||||||||||||||||
|
||||||||||||||||
""" | ||||||||||||||||
if molecules_with_preset_charges is None: | ||||||||||||||||
return None | ||||||||||||||||
|
||||||||||||||||
molecule_set = {molecule.to_smiles() for molecule in molecules_with_preset_charges} | ||||||||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
|
||||||||||||||||
|
||||||||||||||||
if len(molecule_set) != len(molecules_with_preset_charges): | ||||||||||||||||
raise PresetChargesError( | ||||||||||||||||
"All molecules in the molecules_with_preset_charges list must be isomorphically unique from each other", | ||||||||||||||||
) | ||||||||||||||||
|
||||||||||||||||
for molecule in molecules_with_preset_charges: | ||||||||||||||||
if molecule.partial_charges is None: | ||||||||||||||||
raise PresetChargesError( | ||||||||||||||||
"All molecules in the molecules_with_preset_charges list must have partial charges assigned.", | ||||||||||||||||
) | ||||||||||||||||
|
||||||||||||||||
return molecules_with_preset_charges | ||||||||||||||||
|
||||||||||||||||
|
||||||||||||||||
def _create_interchange( | ||||||||||||||||
force_field: ForceField, | ||||||||||||||||
topology: Topology | list[Molecule], | ||||||||||||||||
box: Quantity | None = None, | ||||||||||||||||
positions: Quantity | None = None, | ||||||||||||||||
charge_from_molecules: list[Molecule] | None = None, | ||||||||||||||||
molecules_with_preset_charges: list[Molecule] | None = None, | ||||||||||||||||
partial_bond_orders_from_molecules: list[Molecule] | None = None, | ||||||||||||||||
allow_nonintegral_charges: bool = False, | ||||||||||||||||
) -> Interchange: | ||||||||||||||||
molecules_with_preset_charges = _preprocess_preset_charges(molecules_with_preset_charges) | ||||||||||||||||
|
||||||||||||||||
_check_supported_handlers(force_field) | ||||||||||||||||
|
||||||||||||||||
if molecules_with_preset_charges is not None and "VirtualSites" in force_field.registered_parameter_handlers: | ||||||||||||||||
warnings.warn( | ||||||||||||||||
"Preset charges were provided (via `charge_from_molecules`) alongside a force field that includes " | ||||||||||||||||
"virtual site parameters. Note that virtual sites will be applied charges from the force field and " | ||||||||||||||||
"cannot be given preset charges.", | ||||||||||||||||
Comment on lines
+149
to
+151
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
Do you think adding something like this to the warning will be useful for people? It seems like the sort of thing that would surprise me but maybe it's just noise. |
||||||||||||||||
PresetChargesAndVirtualSitesWarning, | ||||||||||||||||
) | ||||||||||||||||
|
||||||||||||||||
# interchange = Interchange(topology=topology) | ||||||||||||||||
# or maybe | ||||||||||||||||
interchange = Interchange(topology=validate_topology(topology)) | ||||||||||||||||
|
@@ -138,7 +186,7 @@ def _create_interchange( | |||||||||||||||
interchange, | ||||||||||||||||
force_field, | ||||||||||||||||
interchange.topology, | ||||||||||||||||
charge_from_molecules, | ||||||||||||||||
molecules_with_preset_charges, | ||||||||||||||||
allow_nonintegral_charges, | ||||||||||||||||
) | ||||||||||||||||
_plugins(interchange, force_field, interchange.topology) | ||||||||||||||||
|
@@ -274,7 +322,7 @@ def _electrostatics( | |||||||||||||||
interchange: Interchange, | ||||||||||||||||
force_field: ForceField, | ||||||||||||||||
topology: Topology, | ||||||||||||||||
charge_from_molecules: list[Molecule] | None = None, | ||||||||||||||||
molecules_with_preset_charges: list[Molecule] | None = None, | ||||||||||||||||
allow_nonintegral_charges: bool = False, | ||||||||||||||||
): | ||||||||||||||||
if "Electrostatics" not in force_field.registered_parameter_handlers: | ||||||||||||||||
|
@@ -304,7 +352,7 @@ def _electrostatics( | |||||||||||||||
if handler is not None | ||||||||||||||||
], | ||||||||||||||||
topology=topology, | ||||||||||||||||
charge_from_molecules=charge_from_molecules, | ||||||||||||||||
molecules_with_preset_charges=molecules_with_preset_charges, | ||||||||||||||||
allow_nonintegral_charges=allow_nonintegral_charges, | ||||||||||||||||
), | ||||||||||||||||
}, | ||||||||||||||||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The spec uses the same language for
ChargeIncrementModel
s and thecharge_increment
s that are applied to move charges to virtual sites, so I think it's important to keep this distinction very clear. I also think it's clearer to specify the section of the force field that is being applied, instead of just repeating the name.