diff --git a/devtools/conda-envs/dev_env.yaml b/devtools/conda-envs/dev_env.yaml index 6560cb623..2e0ce25bc 100644 --- a/devtools/conda-envs/dev_env.yaml +++ b/devtools/conda-envs/dev_env.yaml @@ -7,19 +7,18 @@ dependencies: - python - pip - versioningit - - pip - numpy <2.3 - pydantic =2 # OpenFF stack - - openff-toolkit-base ~=0.16.6 + - openff-toolkit-base =0.17 - openff-units - ambertools =23 # Optional features - openmm =8.2 # smirnoff-plugins =2024 # de-forcefields # add back after smirnoff-plugins update - - openff-nagl ~=0.5 - - openff-nagl-models ~=0.3 + - openff-nagl + - openff-nagl-models - mbuild ~=0.18 - foyer =1 - gmso ~=0.12 diff --git a/devtools/conda-envs/docs_env.yaml b/devtools/conda-envs/docs_env.yaml index ae2ae91ae..0e66398c3 100644 --- a/devtools/conda-envs/docs_env.yaml +++ b/devtools/conda-envs/docs_env.yaml @@ -10,7 +10,7 @@ dependencies: - pip - numpy =1 - pydantic =2 - - openff-toolkit-base ~=0.16.6 + - openff-toolkit-base =0.17 - openmm =8.2 - mbuild - foyer =1 diff --git a/devtools/conda-envs/examples_env.yaml b/devtools/conda-envs/examples_env.yaml index d896f67e9..85f10b8e8 100644 --- a/devtools/conda-envs/examples_env.yaml +++ b/devtools/conda-envs/examples_env.yaml @@ -9,12 +9,12 @@ dependencies: - numpy <2.3 - pydantic =2 # OpenFF stack - - openff-toolkit-base ~=0.16.6 + - openff-toolkit-base =0.17 - openff-units - ambertools =23 # Optional features - - openff-nagl ~=0.5 - - openff-nagl-models ~=0.3 + - openff-nagl + - openff-nagl-models - mbuild - foyer - nglview diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index d541d7f43..e62a79de9 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -8,7 +8,7 @@ dependencies: - numpy <2.3 - pydantic =2 # OpenFF stack - - openff-toolkit-base ~=0.16.8 + - openff-toolkit-base =0.17 - openff-units =0.3 - ambertools =24 # Needs to be explicitly listed to not be dropped when AmberTools is removed diff --git a/devtools/conda-envs/test_not_py313.yaml b/devtools/conda-envs/test_not_py313.yaml index 800e9fef7..2edeffbe0 100644 --- a/devtools/conda-envs/test_not_py313.yaml +++ b/devtools/conda-envs/test_not_py313.yaml @@ -8,7 +8,7 @@ dependencies: - numpy <2.3 - pydantic =2 # OpenFF stack - - openff-toolkit-base ~=0.16.8 + - openff-toolkit-base =0.17 - openff-units =0.3 - ambertools # Needs to be explicitly listed to not be dropped when AmberTools is removed diff --git a/docs/releasehistory.md b/docs/releasehistory.md index da435df68..44bbb28ee 100644 --- a/docs/releasehistory.md +++ b/docs/releasehistory.md @@ -11,6 +11,12 @@ Dates are given in YYYY-MM-DD format. Please note that all releases prior to a version 1.0.0 are considered pre-releases and many API changes will come before a stable release. +## 0.4.5 - 2025-08-19 + +### New features + +* #1206 Support `` tags in SMIRNOFF force fields. + ## 0.4.4 - 2025-07-29 ### Behavior changes diff --git a/examples/host-guest/host_guest.ipynb b/examples/host-guest/host_guest.ipynb index 3b2afa110..4d17c6a8b 100644 --- a/examples/host-guest/host_guest.ipynb +++ b/examples/host-guest/host_guest.ipynb @@ -89,7 +89,7 @@ "source": [ "NAGLToolkitWrapper().assign_partial_charges(\n", " molecule=host,\n", - " partial_charge_method=\"openff-gnn-am1bcc-0.1.0-rc.1.pt\",\n", + " partial_charge_method=\"openff-gnn-am1bcc-0.1.0-rc.3.pt\",\n", ")\n", "\n", "host.partial_charges.round(3)" @@ -110,7 +110,7 @@ "metadata": {}, "outputs": [], "source": [ - "sage = ForceField(\"openff-2.1.0.offxml\")\n", + "sage = ForceField(\"openff-2.2.1.offxml\")\n", "\n", "out = sage.create_interchange(topology=docked_topology, charge_from_molecules=[host])" ] @@ -192,7 +192,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.11.13" } }, "nbformat": 4, diff --git a/openff/interchange/_tests/conftest.py b/openff/interchange/_tests/conftest.py index a8b505c92..2c3fcd1b3 100644 --- a/openff/interchange/_tests/conftest.py +++ b/openff/interchange/_tests/conftest.py @@ -53,7 +53,7 @@ def sage_with_bond_charge(sage): type="BondCharge", match="all_permutations", distance="0.8 * angstrom ** 1", - charge_increment1="0.0 * elementary_charge ** 1", + charge_increment1="0.123 * elementary_charge ** 1", charge_increment2="0.0 * elementary_charge ** 1", ), ) @@ -187,6 +187,19 @@ def sage_with_off_center_hydrogen(sage): return sage +@pytest.fixture +def sage_nagl(sage): + sage.get_parameter_handler( + "NAGLCharges", + { + "model_file": "openff-gnn-am1bcc-0.1.0-rc.3.pt", + "version": "0.3", + }, + ) + sage.deregister_parameter_handler("ToolkitAM1BCC") + return sage + + @pytest.fixture def _simple_force_field(): # TODO: Create a minimal force field for faster tests @@ -590,7 +603,6 @@ def hydrogen_cyanide_reversed(): def hexane_diol(): molecule = Molecule.from_smiles("OCCCCCCO") molecule.assign_partial_charges(partial_charge_method="gasteiger") - molecule.partial_charges.m return molecule diff --git a/openff/interchange/_tests/unit_tests/smirnoff/test_charge_assignment_logging.py b/openff/interchange/_tests/unit_tests/smirnoff/test_charge_assignment_logging.py index 3e8bedda3..112a95910 100644 --- a/openff/interchange/_tests/unit_tests/smirnoff/test_charge_assignment_logging.py +++ b/openff/interchange/_tests/unit_tests/smirnoff/test_charge_assignment_logging.py @@ -12,10 +12,12 @@ """ Hierarchy is -1. Match molceules with preset charges +1. Match molecules with preset charges 2. Match chemical environment against library charges -3. Match chemical environment against charge increments -4. Run AM1-BCC (or a variant) on remaining molecules +3. Run NAGLCharges (if present in the FF) +3. Run charge method in ChargeIncrementModel section (if present in the FF) and then + match chemical environment against charge increments +4. Run AM1-BCC (or a variant) on remaining molecules (if present in the FF) Test cases ---------- @@ -86,6 +88,19 @@ * Ions get library charges * Ligand gets charge increments +11. Sage with NAGL and ligand in vacuum +* Ligand gets NAGL charges + +12. Sage with NAGL and water molecules (library precedence) +* Water gets library charges (precedence over NAGL) + +13. Sage with NAGL mixed topology (ligand and water) +* Ligand gets NAGL charges +* Water gets library charges + +14. Sage with NAGL and preset charges on ligand +* Ligand gets preset charges (precedence over NAGL) + Other details * Specifics of charge method (AM1-BCC, AM1-BCC-ELF10, AM1-BCC via NAGL) * Molecules with preset charges can be similar but not exact enough @@ -118,6 +133,9 @@ def map_methods_to_atom_indices(caplog: pytest.LogCaptureFixture) -> dict[str, l elif message.startswith("Charge section ToolkitAM1BCC"): info[message.split(", ")[1].split(" ")[-1]].append(int(message.split("atom index ")[-1])) + elif message.startswith("Charge section NAGLCharges"): + info["NAGLChargesHandler"].append(int(message.split("atom index ")[-1])) + # without also pulling the virtual site - particle mapping (which is different for each engine) # it's hard to store more information than the orientation atoms that are affected by each # virtual site's charges @@ -149,7 +167,7 @@ def map_methods_to_atom_indices(caplog: pytest.LogCaptureFixture) -> dict[str, l @pytest.fixture -def sage_with_nagl(sage): +def sage_with_nagl_chargeincrements(sage): from openff.toolkit.typing.engines.smirnoff.parameters import ChargeIncrementModelHandler, ChargeIncrementType sage.register_parameter_handler( @@ -241,6 +259,10 @@ def ligand_and_water_and_ions(ligand, water_and_ions) -> Topology: 8.xSage with preset charges on water 9.xSage with (ligand) virtual site parameters 10.xAM1-with-custom-BCCs Sage with ligand and ions water +11.xSage with NAGL and ligand in vacuum +12.xSage with NAGL and water molecules (library precedence) +13.xSage with NAGL mixed topology (ligand and water) +14.xSage with NAGL and preset charges on ligand """ @@ -416,7 +438,7 @@ def test_case9(caplog, sage_with_bond_charge): assert info["orientation"] == [0, 1] -def test_case10(caplog, sage_with_nagl, ligand): +def test_case10(caplog, sage_with_nagl_chargeincrements, ligand): from openff.toolkit.utils.nagl_wrapper import NAGLToolkitWrapper from openff.toolkit.utils.rdkit_wrapper import RDKitToolkitWrapper from openff.toolkit.utils.toolkit_registry import ToolkitRegistry, toolkit_registry_manager @@ -430,7 +452,7 @@ def test_case10(caplog, sage_with_nagl, ligand): toolkit_precedence=[NAGLToolkitWrapper, RDKitToolkitWrapper], ), ): - sage_with_nagl.create_interchange(ligand.to_topology()) + sage_with_nagl_chargeincrements.create_interchange(ligand.to_topology()) info = map_methods_to_atom_indices(caplog) @@ -445,3 +467,70 @@ def test_case10(caplog, sage_with_nagl, ligand): # the standard AM1-BCC should not have ran assert AM1BCC_KEY not in info + + +def test_case11(caplog, sage_nagl, ligand): + """Test that NAGL charge assignment is properly logged.""" + pytest.importorskip("openff.nagl") + + with caplog.at_level(logging.INFO): + sage_nagl.create_interchange(ligand.to_topology()) + + info = map_methods_to_atom_indices(caplog) + + # Should log NAGL charges for all atoms + assert "NAGLChargesHandler" in info + assert info["NAGLChargesHandler"] == [*range(0, ligand.n_atoms)] + + +def test_case12(caplog, sage_nagl, water): + """Test logging when LibraryCharges takes precedence over NAGL.""" + pytest.importorskip("openff.nagl") + + topology = Topology.from_molecules([water, water]) + + with caplog.at_level(logging.INFO): + sage_nagl.create_interchange(topology) + + info = map_methods_to_atom_indices(caplog) + + # Water should get library charges, not NAGL + assert info["library"] == [*range(0, 6)] # 2 water molecules, 3 atoms each + assert "NAGLChargesHandler" not in info + + +def test_case13(caplog, sage_nagl, ligand, water): + """Test logging with mixed molecule types (some library, some NAGL).""" + pytest.importorskip("openff.nagl") + + topology = Topology.from_molecules([ligand, water]) + + with caplog.at_level(logging.INFO): + sage_nagl.create_interchange(topology) + + info = map_methods_to_atom_indices(caplog) + + # Ligand should get NAGL charges + assert info["NAGLChargesHandler"] == [*range(0, ligand.n_atoms)] + + # Water should get library charges + assert info["library"] == [*range(ligand.n_atoms, ligand.n_atoms + water.n_atoms)] + + +def test_case14(caplog, sage_nagl, ligand): + """Test logging when preset charges are used instead of NAGL.""" + pytest.importorskip("openff.nagl") + + ligand.assign_partial_charges("gasteiger") + + with caplog.at_level(logging.INFO): + sage_nagl.create_interchange( + ligand.to_topology(), + charge_from_molecules=[ligand], + ) + + info = map_methods_to_atom_indices(caplog) + + # Should use preset charges, not NAGL + assert info["preset"] == [*range(0, ligand.n_atoms)] + assert "NAGLChargesHandler" not in info diff --git a/openff/interchange/_tests/unit_tests/smirnoff/test_nonbonded.py b/openff/interchange/_tests/unit_tests/smirnoff/test_nonbonded.py index fe3af8bf1..1a03bbc29 100644 --- a/openff/interchange/_tests/unit_tests/smirnoff/test_nonbonded.py +++ b/openff/interchange/_tests/unit_tests/smirnoff/test_nonbonded.py @@ -1,14 +1,16 @@ import numpy import pytest -from openff.toolkit import Molecule, Quantity, Topology, unit +from openff.toolkit import ForceField, Molecule, Quantity, RDKitToolkitWrapper, Topology, unit from openff.toolkit.typing.engines.smirnoff import ( ChargeIncrementModelHandler, ElectrostaticsHandler, LibraryChargeHandler, + NAGLChargesHandler, ToolkitAM1BCCHandler, vdWHandler, ) -from openff.toolkit.utils.exceptions import SMIRNOFFVersionError +from openff.toolkit.utils.exceptions import MissingPackageError, SMIRNOFFVersionError +from openff.toolkit.utils.toolkit_registry import ToolkitRegistry, toolkit_registry_manager from packaging.version import Version from openff.interchange import Interchange @@ -135,6 +137,501 @@ def test_toolkit_am1bcc_uses_elf10_if_oe_is_available(self, sage, hexane_diol): assert not uses_elf10 numpy.testing.assert_allclose(partial_charges, assigned_charges) + def test_nagl_charge_assignment_matches_reference(self, sage_nagl, hexane_diol): + hexane_diol.assign_partial_charges("openff-gnn-am1bcc-0.1.0-rc.3.pt") + # Leave the ToolkitAM1BCC tag in openff-2.1.0 to ensure that the NAGLCharges handler takes precedence + + interchange = sage_nagl.create_interchange(topology=hexane_diol.to_topology()) + + assigned_charges_unitless = interchange["Electrostatics"].get_charge_array().m + + expected_charges = hexane_diol.partial_charges + assert expected_charges is not None + assert expected_charges.units == unit.elementary_charge + assert not all(charge == 0 for charge in expected_charges.magnitude) + expected_charges_unitless = [v.m for v in expected_charges] + numpy.testing.assert_allclose(expected_charges_unitless, assigned_charges_unitless) + + +class TestNAGLChargesErrorHandling: + """Test NAGLCharges error conditions.""" + + def test_nagl_charges_missing_toolkit_error(self, sage_nagl, hexane_diol): + """Test MissingPackageError when NAGL toolkit is not available. This should fail immediately instead of falling + back to ToolkitAM1BCC, since it doesn't know whether the molecule would have successfully + had charges assigned by NAGL if it were available.""" + + # Mock the toolkit registry to not have NAGL + # RDKit is needed for SMARTS matching. + with toolkit_registry_manager(ToolkitRegistry(toolkit_precedence=[RDKitToolkitWrapper])): + with pytest.raises(MissingPackageError, match="NAGL software isn't present"): + sage_nagl.create_interchange(topology=hexane_diol.to_topology()) + + # No error should be raised if using charge_from_molecules + sage_nagl.create_interchange( + topology=hexane_diol.to_topology(), + charge_from_molecules=[hexane_diol], + ) + + def test_nagl_charges_invalid_model_file(self, sage, hexane_diol): + """Test error handling for invalid model file paths. This should fail immediately instead of falling + back to ToolkitAM1BCC, since it doesn't know whether the molecule would have successfully + had charges assigned by this model it had been able to find it.""" + sage.get_parameter_handler( + "NAGLCharges", + { + "model_file": "nonexistent_model.pt", + "version": "0.3", + }, + ) + with pytest.raises(FileNotFoundError): + sage.create_interchange(topology=hexane_diol.to_topology()) + + sage["NAGLCharges"].model_file = "" + with pytest.raises(FileNotFoundError): + sage.create_interchange(topology=hexane_diol.to_topology()) + + sage["NAGLCharges"].model_file = None + with pytest.raises(FileNotFoundError): + sage.create_interchange(topology=hexane_diol.to_topology()) + + def test_nagl_charges_bad_hash(self, sage, hexane_diol, monkeypatch): + """Test error handling for a bad hash. This should fail immediately instead of falling + back to ToolkitAM1BCC, since it doesn't know whether the molecule would have successfully + had charges assigned by this model if the hash comparison hadn't failed.""" + from openff.nagl_models._dynamic_fetch import HashComparisonFailedException + + sage.get_parameter_handler( + "NAGLCharges", + { + "model_file": "openff-gnn-am1bcc-0.1.0-rc.3.pt", + "model_file_hash": "bad_hash", + "version": "0.3", + }, + ) + with pytest.raises(HashComparisonFailedException): + sage.create_interchange(topology=hexane_diol.to_topology()) + + def test_nagl_charges_bad_doi(self, sage, hexane_diol, monkeypatch): + """Test error handling for a bad DOI. This should fail immediately instead of falling + back to ToolkitAM1BCC, since it doesn't know whether the molecule would have successfully + had charges assigned by this model, since it's unfetchable.""" + from openff.nagl_models._dynamic_fetch import UnableToParseDOIException + + sage.get_parameter_handler( + "NAGLCharges", + { + "model_file": "nonexistent_model.pt", + "digital_object_identifier": "blah.foo/bar", + "version": "0.3", + }, + ) + with pytest.raises(UnableToParseDOIException): + sage.create_interchange(topology=hexane_diol.to_topology()) + + # For more information on why this test is skipped, see + # https://github.com/openforcefield/openff-interchange/pull/1206/commits/f99a10e17ad56235ba1f36ae35f6383a22ed840a#r2248864028 + @pytest.mark.xfail( + reason="charge assignment handler fallback behavior not yet implemented", + raises=ValueError, + ) + def test_nagl_charges_fallback_to_charge_increment_model(self, sage): + """Test that NAGL falls back to ChargeIncrementModel when molecule contains unsupported elements.""" + pytest.importorskip("openff.nagl") + + # Create a boron-containing molecule with nonzero formal charge + # BF4- anion - boron is not supported by current NAGL models + boron_molecule = Molecule.from_smiles("[B-]([F])([F])([F])[F]") + + # Verify formal charges are not all zero + formal_charges = [atom.formal_charge.m for atom in boron_molecule.atoms] + assert not all(charge == 0 for charge in formal_charges) + + # Create minimal force field with only the needed handlers + ff = ForceField() + + # Add Electrostatics handler + ff.register_parameter_handler( + ElectrostaticsHandler(version="0.4"), + ) + + # Add NAGLCharges handler + ff.register_parameter_handler( + NAGLChargesHandler( + version="0.3", + model_file="openff-gnn-am1bcc-0.1.0-rc.3.pt", + ), + ) + + # Add ChargeIncrementModel handler with formal_charge method and no increments + charge_increment_handler = ChargeIncrementModelHandler( + version="0.3", + partial_charge_method="formal_charge", + ) + ff.register_parameter_handler(charge_increment_handler) + + # Should succeed despite NAGL not supporting boron + interchange = ff.create_interchange(topology=boron_molecule.to_topology()) + + # Should have assigned charges to all atoms + assigned_charges = interchange["Electrostatics"].get_charge_array() + + # Assigned charges should match formal charges (fallback to ChargeIncrementModel) + expected_charges = [atom.formal_charge.m for atom in boron_molecule.atoms] + numpy.testing.assert_allclose(assigned_charges.m, expected_charges) + + # Net charge should match molecule's total formal charge + assert abs(sum(assigned_charges.m) - boron_molecule.total_charge.m) < 1e-10 + + @pytest.mark.xfail( + reason="charge assignment handler fallback behavior not yet implemented", + raises=ValueError, + ) + def test_nagl_charges_all_handlers_fail_comprehensive_error(self, sage): + """Test error reporting when all charge assignment methods fail.""" + pytest.importorskip("openff.nagl") + + # Create a uranium compound - not supported by any current charge assignment method + uranium_molecule = Molecule.from_smiles("[U+4]") + + # Create force field with multiple charge assignment handlers + ff = ForceField() + + # Add Electrostatics handler + ff.register_parameter_handler( + ElectrostaticsHandler(version="0.4"), + ) + + # Add NAGLCharges handler + ff.register_parameter_handler( + NAGLChargesHandler( + version="0.3", + model_file="openff-gnn-am1bcc-0.1.0-rc.3.pt", + ), + ) + + # Add ToolkitAM1BCC handler + ff.register_parameter_handler( + ToolkitAM1BCCHandler(version="0.3"), + ) + + # Add ChargeIncrementModel handler with gasteiger method + charge_increment_handler = ChargeIncrementModelHandler( + version="0.3", + partial_charge_method="mmff94", + ) + ff.register_parameter_handler(charge_increment_handler) + + # Should fail with comprehensive error message + with pytest.raises(RuntimeError) as excinfo: + ff.create_interchange(topology=uranium_molecule.to_topology()) + + error_message = str(excinfo.value) + + # Error should mention that no charges could be assigned + assert "could not be fully assigned charges" in error_message + + # Error should contain information about each handler's failure + assert "NAGLCharges" in error_message + assert "ToolkitAM1BCC" in error_message + assert "ChargeIncrementModel" in error_message + + # Should mention the exceptions raised by various handlers + assert "exceptions raised by various handlers" in error_message + + +class TestNAGLChargesPrecedence: + """Test NAGLCharges precedence in the hierarchy of charge assignment methods.""" + + def test_nagl_charges_precedence_over_am1bcc(self, sage_nagl, hexane_diol): + """Test that NAGLCharges takes precedence over ToolkitAM1BCC.""" + sage_nagl.get_parameter_handler("ToolkitAM1BCC", {"version": "0.3"}) + # Get reference charges from NAGL + hexane_diol.assign_partial_charges("openff-gnn-am1bcc-0.1.0-rc.3.pt") + nagl_charges = [c.m for c in hexane_diol.partial_charges] + + # Get reference charges from AM1BCC + hexane_diol.assign_partial_charges("am1bcc") + am1bcc_charges = [c.m for c in hexane_diol.partial_charges] + + # Ensure they're different + assert not numpy.allclose(nagl_charges, am1bcc_charges) + + interchange = sage_nagl.create_interchange(topology=hexane_diol.to_topology()) + assigned_charges = interchange["Electrostatics"].get_charge_array() + + # Should match NAGL charges, not AM1BCC + numpy.testing.assert_allclose(assigned_charges, nagl_charges) + + def test_library_charges_precedence_over_nagl(self, sage_nagl, methane): + """Test that LibraryCharges takes precedence over NAGLCharges.""" + + sage_nagl["LibraryCharges"].add_parameter( + { + "smirks": "[#6X4:1]-[#1:2]", + "charge1": -0.2 * unit.elementary_charge, + "charge2": 0.05 * unit.elementary_charge, + }, + ) + + interchange = sage_nagl.create_interchange(topology=methane.to_topology()) + assigned_charges = interchange["Electrostatics"].get_charge_array() + + # Should match library charges + expected_charges = [-0.2, 0.05, 0.05, 0.05, 0.05] + numpy.testing.assert_allclose(assigned_charges, expected_charges) + + def test_nagl_charges_precedence_over_charge_increments(self, sage_nagl, hexane_diol): + """Test that NAGLCharges takes precedence over ChargeIncrementModel as base charges.""" + + # Get reference charges from NAGL + hexane_diol.assign_partial_charges("openff-gnn-am1bcc-0.1.0-rc.3.pt") + nagl_charges = [c.m for c in hexane_diol.partial_charges] + + # Add ChargeIncrementModel handler (should provide base charges, not increments) + increment_handler = ChargeIncrementModelHandler( + version=0.3, + partial_charge_method="formal_charge", + ) + sage_nagl.register_parameter_handler(increment_handler) + + interchange = sage_nagl.create_interchange(topology=hexane_diol.to_topology()) + assigned_charges = interchange["Electrostatics"].get_charge_array() + + # Should match NAGL charges, not formal charges + numpy.testing.assert_allclose(assigned_charges, nagl_charges) + + +class TestNAGLChargesIntegration: + """Test NAGLCharges integration with other handlers.""" + + def test_nagl_charges_multi_molecule_topology(self, sage_nagl): + """Test NAGLCharges with multiple molecules in topology.""" + methane = Molecule.from_smiles("C") + ethane = Molecule.from_smiles("CC") + + topology = Topology.from_molecules([methane, ethane]) + + interchange = sage_nagl.create_interchange(topology=topology) + assigned_charges = interchange["Electrostatics"].get_charge_array() + + # Each molecule should have approximately zero net charge + methane_charge_sum = sum(assigned_charges[: methane.n_atoms]) + ethane_charge_sum = sum(assigned_charges[methane.n_atoms :]) + + assert abs(methane_charge_sum) < 1e-10 * unit.elementary_charge + assert abs(ethane_charge_sum) < 1e-10 * unit.elementary_charge + + def test_nagl_charges_with_virtual_sites(self, sage_with_bond_charge): + """Test NAGLCharges compatibility with virtual sites.""" + + # Create a molecule that would have virtual sites + molecule = Molecule.from_smiles("[Cl]CCO") + + # Add NAGLCharges to the force field + sage_with_bond_charge.get_parameter_handler( + "NAGLCharges", + { + "model_file": "openff-gnn-am1bcc-0.1.0-rc.3.pt", + "version": "0.3", + }, + ) + + # Should not raise an error + interchange = sage_with_bond_charge.create_interchange( + topology=molecule.to_topology(), + ) + + # Should have charges for real atoms + assigned_charges = interchange["Electrostatics"]._get_charges() + assert len(assigned_charges.values()) - 1 == molecule.n_atoms + + # Net charge should be approximately zero + all_particle_charge_sum = sum(assigned_charges.values()) + assert abs(all_particle_charge_sum) < 1e-10 * unit.elementary_charge + # Charge without the vsite should be nonzero + atom_charge_sum = sum([charge for tk, charge in assigned_charges.items() if tk.atom_indices is not None]) + assert abs(atom_charge_sum - (0.123 * unit.elementary_charge)) < 1e-10 * unit.elementary_charge + + def test_nagl_charges_force_field_creation_complete(self, hexane_diol): + """Test complete interchange creation with NAGLCharges.""" + + ff = ForceField("openff-2.1.0.offxml") + ff.get_parameter_handler( + "NAGLCharges", + { + "model_file": "openff-gnn-am1bcc-0.1.0-rc.3.pt", + "version": "0.3", + }, + ) + + # Should create complete interchange without errors + interchange = ff.create_interchange(topology=hexane_diol.to_topology()) + + # Should have all expected collections + expected_collections = ["Bonds", "Angles", "ProperTorsions", "ImproperTorsions", "vdW", "Electrostatics"] + for collection_name in expected_collections: + assert collection_name in interchange.collections + + # Electrostatics should have charges + charges = interchange["Electrostatics"].get_charge_array() + assert len(charges) == hexane_diol.n_atoms + + # Net charge should be approximately zero + total_charge = sum(charge.m for charge in charges) + assert abs(total_charge) < 1e-10 + + def test_nagl_charges_identical_molecules_same_charges(self): + """Test that identical molecules get identical charges from NAGLCharges.""" + + # Create topology with two identical molecules + molecule1 = Molecule.from_smiles("CCO") + molecule2 = Molecule.from_smiles("CCO") + topology = Topology.from_molecules([molecule1, molecule2]) + + ff = ForceField("openff-2.1.0.offxml") + ff.get_parameter_handler( + "NAGLCharges", + { + "model_file": "openff-gnn-am1bcc-0.1.0-rc.3.pt", + "version": "0.3", + }, + ) + + interchange = ff.create_interchange(topology=topology) + assigned_charges = interchange["Electrostatics"].get_charge_array() + + # First molecule charges + mol1_charges = assigned_charges[: molecule1.n_atoms] + # Second molecule charges + mol2_charges = assigned_charges[molecule1.n_atoms :] + + # Should be identical + numpy.testing.assert_allclose(mol1_charges, mol2_charges) + + def test_nagl_charges_with_charge_from_molecules(self, sage_nagl, hexane_diol): + """Test that charge_from_molecules takes precedence over NAGLCharges.""" + # Assign preset charges using a different method + hexane_diol.assign_partial_charges("gasteiger") + preset_charges = [c.m for c in hexane_diol.partial_charges] + + # Create interchange with charge_from_molecules - should use preset charges + interchange = sage_nagl.create_interchange( + topology=hexane_diol.to_topology(), + charge_from_molecules=[hexane_diol], + ) + + assigned_charges = interchange["Electrostatics"].get_charge_array() + + # Should match preset charges, not NAGL charges + numpy.testing.assert_allclose(assigned_charges.m, preset_charges) + + # Verify NAGL would give different charges + hexane_diol_copy = Molecule.from_smiles(hexane_diol.to_smiles()) + hexane_diol_copy.assign_partial_charges("openff-gnn-am1bcc-0.1.0-rc.3.pt") + nagl_charges = [c.m for c in hexane_diol_copy.partial_charges] + + # Preset and NAGL charges should be different + assert not numpy.allclose(preset_charges, nagl_charges, atol=1e-3) + + def test_nagl_charges_with_mixed_charge_sources(self, sage_nagl): + """Test NAGLCharges with some molecules having preset charges and others not.""" + # Create molecules + ethanol = Molecule.from_smiles("CCO") + methanol = Molecule.from_smiles("CO") + + # Assign preset charges to only one molecule + ethanol.assign_partial_charges("gasteiger") + preset_ethanol_charges = [c.m for c in ethanol.partial_charges] + + topology = Topology.from_molecules([ethanol, methanol]) + + # Create interchange with preset charges for ethanol only + interchange = sage_nagl.create_interchange( + topology=topology, + charge_from_molecules=[ethanol], + ) + + assigned_charges = interchange["Electrostatics"].get_charge_array() + + # First molecule (ethanol) should match preset charges + ethanol_charges = assigned_charges[: ethanol.n_atoms] + numpy.testing.assert_allclose(ethanol_charges.m, preset_ethanol_charges) + + # Second molecule (methanol) should get NAGL charges + methanol_charges = assigned_charges[ethanol.n_atoms :] + + # Get reference NAGL charges for methanol + methanol_copy = Molecule.from_smiles("CO") + methanol_copy.assign_partial_charges("openff-gnn-am1bcc-0.1.0-rc.3.pt") + nagl_methanol_charges = [c.m for c in methanol_copy.partial_charges] + + numpy.testing.assert_allclose(methanol_charges.m, nagl_methanol_charges) + + @pytest.mark.slow + def test_nagl_charges_large_molecule_performance(self, sage_nagl): + """Test that NAGL charge assignment completes in reasonable time for large molecules.""" + import time + + # Create a very large molecule + large_molecule = Molecule.from_smiles("C" * 200) # 200-carbon alkane chain + + start_time = time.time() + + # Should complete without error + interchange = sage_nagl.create_interchange(topology=large_molecule.to_topology()) + + end_time = time.time() + execution_time = end_time - start_time + + # Should complete within reasonable time (less than 30 seconds) + assert execution_time < 30.0, f"NAGL charge assignment took {execution_time:.2f}s, which is too long" + + # Net charge should be approximately zero + charges = interchange["Electrostatics"].get_charge_array() + total_charge = sum(charges.m) + assert abs(total_charge) < 1e-10 + + @pytest.mark.slow + def test_nagl_charges_multiple_large_molecules_performance(self, sage_nagl): + """Test performance with multiple large molecules in topology.""" + import time + + # Create multiple copies of medium-sized molecules + base_molecules = [ + Molecule.from_smiles("C" * 20), # 20-carbon chain + Molecule.from_smiles("C" * 25), # 25-carbon chain + Molecule.from_smiles("C" * 30), # 30-carbon chain + ] + + # Create 20 copies of each + molecules = [] + for _ in range(20): + for base_mol in base_molecules: + molecules.append(base_mol) + + topology = Topology.from_molecules(molecules) + + start_time = time.time() + + # Should complete without error + interchange = sage_nagl.create_interchange(topology=topology) + + end_time = time.time() + execution_time = end_time - start_time + + # Should complete within reasonable time + assert execution_time < 30.0, f"Multi-molecule NAGL assignment took {execution_time:.2f}s, which is too long" + + # Each molecule should have approximately zero net charge + charges = interchange["Electrostatics"].get_charge_array() + start_idx = 0 + for molecule in molecules: + mol_charges = charges[start_idx : start_idx + molecule.n_atoms] + mol_total_charge = sum(mol_charges.m) + assert abs(mol_total_charge) < 1e-10 + start_idx += molecule.n_atoms + @pytest.mark.skip( reason="Turn on if toolkit ever allows non-standard scale12/13/15", ) diff --git a/openff/interchange/smirnoff/_create.py b/openff/interchange/smirnoff/_create.py index 5a6da6d65..31cdb0b05 100644 --- a/openff/interchange/smirnoff/_create.py +++ b/openff/interchange/smirnoff/_create.py @@ -69,7 +69,7 @@ def _check_supported_handlers(force_field: ForceField): unsupported = list() for handler_name in force_field.registered_parameter_handlers: - if handler_name in {"ToolkitAM1BCC"}: + if handler_name in {"ToolkitAM1BCC", "NAGLCharges"}: continue if handler_name not in _SUPPORTED_PARAMETER_HANDLERS: unsupported.append(handler_name) @@ -348,6 +348,7 @@ def _electrostatics( force_field._parameter_handlers.get(name, None) for name in [ "Electrostatics", + "NAGLCharges", "ChargeIncrementModel", "ToolkitAM1BCC", "LibraryCharges", diff --git a/openff/interchange/smirnoff/_nonbonded.py b/openff/interchange/smirnoff/_nonbonded.py index e97d5f55e..92067d272 100644 --- a/openff/interchange/smirnoff/_nonbonded.py +++ b/openff/interchange/smirnoff/_nonbonded.py @@ -11,10 +11,12 @@ ChargeIncrementModelHandler, ElectrostaticsHandler, LibraryChargeHandler, + NAGLChargesHandler, ParameterHandler, ToolkitAM1BCCHandler, vdWHandler, ) +from openff.toolkit.utils.exceptions import MissingPackageError from pydantic import Field, PrivateAttr, computed_field from typing_extensions import Self @@ -48,6 +50,7 @@ ElectrostaticsHandlerType = Union[ ElectrostaticsHandler, + NAGLChargesHandler, ToolkitAM1BCCHandler, ChargeIncrementModelHandler, LibraryChargeHandler, @@ -248,7 +251,7 @@ class SMIRNOFFElectrostaticsCollection(ElectrostaticsCollection, SMIRNOFFCollect * global settings for the electrostatic interactions such as the cutoff distance and the intramolecular scale factors. * partial charges which have been assigned by a ``ToolkitAM1BCC``, - ``LibraryCharges``, or a ``ChargeIncrementModel`` parameter + ``LibraryCharges``, ``NAGLCharges``, or a ``ChargeIncrementModel`` parameter handler. * charge corrections applied by a ``ChargeIncrementHandler`` @@ -284,6 +287,7 @@ def allowed_parameter_handlers(cls): """Return a list of allowed types of ParameterHandler classes.""" return [ LibraryChargeHandler, + NAGLChargesHandler, ChargeIncrementModelHandler, ToolkitAM1BCCHandler, ElectrostaticsHandler, @@ -384,6 +388,7 @@ def _get_charges( if potential_key.associated_handler in ( "LibraryCharges", + "NAGLChargesHandler", "ToolkitAM1BCCHandler", "molecules_with_preset_charges", "ExternalSource", @@ -450,7 +455,7 @@ def parameter_handler_precedence(cls) -> list[str]: """ Return the order in which parameter handlers take precedence when computing charges. """ - return ["LibraryCharges", "ChargeIncrementModel", "ToolkitAM1BCC"] + return ["LibraryCharges", "NAGLCharges", "ChargeIncrementModel", "ToolkitAM1BCC"] @classmethod def create( @@ -510,10 +515,21 @@ def _compute_partial_charges( molecule: Molecule, mapped_smiles: str, method: str, + exception_types_to_raise: Iterable[Exception], + additional_args: tuple[tuple[str, str]], ) -> Quantity: """Call out to the toolkit's toolkit wrappers to generate partial charges.""" + from openff.toolkit.utils.toolkits import GLOBAL_TOOLKIT_REGISTRY + + additional_args = {i: j for i, j in additional_args} molecule = copy.deepcopy(molecule) - molecule.assign_partial_charges(method) + GLOBAL_TOOLKIT_REGISTRY.call( + "assign_partial_charges", + molecule=molecule, + partial_charge_method=method, + raise_exception_types=exception_types_to_raise, + **additional_args, + ) return molecule.partial_charges @@ -662,7 +678,7 @@ def _find_slot_matches( @classmethod def _find_charge_model_matches( cls, - parameter_handler: ToolkitAM1BCCHandler | ChargeIncrementModelHandler, + parameter_handler: ToolkitAM1BCCHandler | ChargeIncrementModelHandler | NAGLChargesHandler, unique_molecule: Molecule, ) -> tuple[ str, @@ -670,6 +686,8 @@ def _find_charge_model_matches( dict[PotentialKey, Potential], ]: """Construct a slot and potential map for a charge model based parameter handler.""" + from openff.nagl_models._dynamic_fetch import HashComparisonFailedException, UnableToParseDOIException + unique_molecule = copy.deepcopy(unique_molecule) reference_smiles = unique_molecule.to_smiles( isomeric=True, @@ -679,8 +697,29 @@ def _find_charge_model_matches( handler_name = parameter_handler.__class__.__name__ - if handler_name == "ChargeIncrementModelHandler": + if handler_name == "NAGLChargesHandler": + from openff.toolkit.utils.toolkits import GLOBAL_TOOLKIT_REGISTRY, NAGLToolkitWrapper + + partial_charge_method = parameter_handler.model_file + if NAGLToolkitWrapper not in {type(tk) for tk in GLOBAL_TOOLKIT_REGISTRY.registered_toolkits}: + raise MissingPackageError( + "The force field has a NAGLCharges section, but the NAGL software isn't " + "present in GLOBAL_TOOLKIT_REGISTRY", + ) + exception_types_to_raise = ( + FileNotFoundError, + MissingPackageError, + HashComparisonFailedException, + UnableToParseDOIException, + ) + additional_args = ( + ("doi", parameter_handler.digital_object_identifier), + ("file_hash", parameter_handler.model_file_hash), + ) + elif handler_name == "ChargeIncrementModelHandler": partial_charge_method = parameter_handler.partial_charge_method + exception_types_to_raise = tuple() + additional_args = tuple() elif handler_name == "ToolkitAM1BCCHandler": from openff.toolkit.utils.toolkits import GLOBAL_TOOLKIT_REGISTRY @@ -690,12 +729,13 @@ def _find_charge_model_matches( partial_charge_method = "am1bccelf10" else: partial_charge_method = "am1bcc" + exception_types_to_raise = tuple() + additional_args = tuple() else: raise InvalidParameterHandlerError( f"Encountered unknown handler of type {type(parameter_handler)} where only " - "ToolkitAM1BCCHandler or ChargeIncrementModelHandler are expected.", + "ToolkitAM1BCCHandler, NAGLChargesHandler, or ChargeIncrementModelHandler are expected.", ) - partial_charges = cls._compute_partial_charges( unique_molecule, unique_molecule.to_smiles( @@ -704,8 +744,9 @@ def _find_charge_model_matches( mapped=True, ), method=partial_charge_method, + exception_types_to_raise=exception_types_to_raise, + additional_args=additional_args, ) - matches = {} potentials = {} @@ -761,7 +802,7 @@ def _find_reference_matches( unique_molecule, ) - if handler_type in ["ToolkitAM1BCC", "ChargeIncrementModel"]: + if handler_type in ["ToolkitAM1BCC", "ChargeIncrementModel", "NAGLCharges"]: ( partial_charge_method, am1_matches, @@ -942,6 +983,12 @@ def store_matches( f"{new_key.extras['partial_charge_method']}, " f"applied to topology atom index {topology_atom_index}", ) + elif new_key.extras["handler"] == "NAGLChargesHandler": + logger.info( + "Charge section NAGLCharges, using NAGL model " + f"{new_key.extras['partial_charge_method']}, " + f"applied to topology atom index {topology_atom_index}", + ) elif new_key.extras["handler"] == "preset": logger.info(