diff --git a/The-SMIRNOFF-force-field-format.md b/The-SMIRNOFF-force-field-format.md index dd00659e3..48d4cfc3a 100644 --- a/The-SMIRNOFF-force-field-format.md +++ b/The-SMIRNOFF-force-field-format.md @@ -225,30 +225,27 @@ For example, to ensure water molecules are assigned partial charges for [TIP3P]( ### ``: Small molecule and fragment charges -.. warning:: This functionality is not yet implemented and will appear in a future version of the toolkit. This area of the SMIRNOFF spec is under further consideration. Please see `Issue 208 on the Open Force Field Toolkit issue tracker `_. - In keeping with the AMBER force field philosophy, especially as implemented in small molecule force fields such as [GAFF](http://ambermd.org/antechamber/gaff.html), [GAFF2](https://mulan.swmed.edu/group/gaff.php), and [parm@Frosst](http://www.ccl.net/cca/data/parm_at_Frosst/), partial charges for small molecules are usually assigned using a quantum chemical method (usually a semiempirical method such as [AM1](https://en.wikipedia.org/wiki/Austin_Model_1)) and a [partial charge determination scheme](https://en.wikipedia.org/wiki/Partial_charge) (such as [CM2](http://doi.org/10.1021/jp972682r) or [RESP](http://doi.org/10.1021/ja00074a030)), then subsequently corrected via charge increment rules, as in the highly successful [AM1-BCC](https://dx.doi.org/10.1002/jcc.10128) approach. Here is an example: ```XML - + - - - - - + + + + + ``` The sum of formal charges for the molecule or fragment will be used to determine the total charge the molecule or fragment will possess. `` provides several optional attributes to control its behavior: -* The `number_of_conformers` attribute (default: `"10"`) is used to specify how many conformers will be generated for the molecule (or capped fragment) prior to charging. -* The `quantum_chemical_method` attribute (default: `"AM1"`) is used to specify the quantum chemical method applied to the molecule or capped fragment. -* The `partial_charge_method` attribute (default: `"CM2"`) is used to specify how uncorrected partial charges are to be generated from the quantum chemical wavefunction. Later additions will add restrained electrostatic potential fitting (RESP) capabilities. +* The `number_of_conformers` attribute (default: `"1"`) is used to specify how many conformers will be generated for the molecule (or capped fragment) prior to charging. +* The `partial_charge_method` attribute (default: `"AM1-Mulliken"`) is used to specify how uncorrected partial charges are to be generated. Later additions will add restrained electrostatic potential fitting (RESP) capabilities. The `` tags specify how the quantum chemical derived charges are to be corrected to produce the final charges. -The `chargeincrement#` attributes specify how much the charge on the associated tagged atom index (replacing `#`) should be modified. +The `charge_increment#` attributes specify how much the charge on the associated tagged atom index (replacing `#`) should be modified. The sum of charge increments should equal zero. Note that atoms for which library charges have already been applied are excluded from charging via ``. @@ -654,7 +651,7 @@ We will support the following different types or geometries of off-center charge - `TrivalentLonePair`: This is suitable for planar or tetrahedral nitrogen lone pairs; a charge site `S` lies above the central atom (e.g. nitrogen, blue) a distance `d` along the vector perpendicular to the plane of the three connected atoms (2,3,4). With positive values of `d` the site lies above the nitrogen and with negative values it lies below the nitrogen. ![Trivalent lone pair virtual site](figures/vsite_trivalent.jpg) -Each virtual site receives charge which is transferred from the desired atoms specified in the SMIRKS pattern via a `chargeincrement#` parameter, e.g., if `chargeincrement1=+0.1*elementary_charge` then the virtual site will receive a charge of -0.1 and the atom labeled `1` will have its charge adjusted upwards by +0.1. +Each virtual site receives charge which is transferred from the desired atoms specified in the SMIRKS pattern via a `charge_increment#` parameter, e.g., if `charge_increment1=+0.1*elementary_charge` then the virtual site will receive a charge of -0.1 and the atom labeled `1` will have its charge adjusted upwards by +0.1. N may index any indexed atom. Increments which are left unspecified default to zero. Additionally, each virtual site can bear Lennard-Jones parameters, specified by `sigma` and `epsilon` or `rmin_half` and `epsilon`. @@ -665,23 +662,23 @@ In the SMIRNOFF format, these are encoded as: - + - + - + - + - - + + ``` diff --git a/docs/releasehistory.rst b/docs/releasehistory.rst index f1f828bc1..fb9355929 100644 --- a/docs/releasehistory.rst +++ b/docs/releasehistory.rst @@ -10,14 +10,53 @@ Releases follow the ``major.minor.micro`` scheme recommended by `PEP440 `_: In order to provide the same results for the same chemical species, regardless of input - conformation, fractional bond order calculation methods now default to ignore input conformers - and generate a new conformer of the molecule before running semiempirical calculations. + conformation, ``assign_partial_charges``, ``compute_partial_charges_am1bcc``, and + ``assign_fractional_bond_orders`` methods now default to ignore input conformers + and generate new conformer(s) of the molecule before running semiempirical calculations. Users can override this behavior by specifying the keyword argument - ``use_conformers=molecule.conformers`` + ``use_conformers=molecule.conformers``. - `PR #281 `_: Closes `Issue #250 `_ by adding support for partial charge I/O in SDF. The partial charges are stored as a property in the @@ -58,10 +97,38 @@ Behavior changed it has failed to assign charges. - `PR #597 `_: Energy-minimized sample systems with Parsley 1.1.0. - +- `PR #558 `_: The + :py:class`Topology ` + particle indexing system now orders :py:class`TopologyVirtualSites ` + after all atoms. +- `PR #469 `_: + When running :py:meth:`Topology.to_openmm `, unique atom names + are generated if the provided atom names are not unique (overriding any existing atom names). This + uniqueness extends only to atoms in the same molecule. To disable this behavior, set the kwarg + ``ensure_unique_atom_names=False``. +- `PR #472 `_: + The :py:meth:`Molecule.__eq__ ` now uses the new + :py:meth:`Molecule.are_isomorphic ` to perform the + similarity checking. +- `PR #472 `_: + The :py:meth:`Topology.from_openmm ` and + :py:meth:`Topology.add_molecule ` now use the + :py:meth:`Molecule.are_isomorphic ` to match + molecules. +- `PR #544 `_: Raises + ``NotImplementedError`` when calling + :py:meth:`ParameterHandler.get_parameter `, + which is not yet implemented, but would previously silently return ``None``. +- `PR #551 `_: Implemented the + :py:meth:`ParameterHandler.get_parameter ` function. API-breaking changes """""""""""""""""""" +- `PR #471 `_: Closes + `Issue #465 `_. + ``atom.formal_charge`` and ``molecule.total_charge`` now return ``simtk.unit.Quantity`` objects + instead of integers. To preserve backward compatibility, the setter for ``atom.formal_charge`` + can accept either a ``simtk.unit.Quantity`` or an integer. - `PR #601 `_: Removes almost all of the previous :py:class:`ChemicalEnvironment ` @@ -99,6 +166,28 @@ API-breaking changes New features """""""""""" +- `PR #471 `_: Closes + `Issue #208 `_ + by implementing support for the + ``ChargeIncrementModel`` tag in the `SMIRNOFF specification `_. +- `PR #471 `_: Implements + ``Molecule.compute_partial_charges``, which calls one of the newly-implemented + ``OpenEyeToolkitWrapper.compute_partial_charges``, and + ``AmberToolsToolkitWrapper.compute_partial_charges``. ``strict_n_conformers`` is a + optional boolean parameter indicating whether an ``IncorrectNumConformersError`` should be raised if an invalid + number of conformers is supplied during partial charge calculation. For example, if two conformers are + supplied, but ``partial_charge_method="AM1BCC"`` is also set, then there is no clear use for + the second conformer. The previous behavior in this case was to raise a warning, and to preserve that + behavior, ``strict_n_conformers`` defaults to a value of ``False``. +- `PR #471 `_: Adds + keyword argument ``raise_exception_types`` (default: ``[Exception]``) to + :py:meth:`ToolkitRegistry.call . + The default value will provide the previous OpenFF Toolkit behavior, which is that the first ToolkitWrapper + that can provide the requested method is called, and it either returns on success or raises an exception. This new + keyword argument allows the ToolkitRegistry to *ignore* certain exceptions, but treat others as fatal. + If ``raise_exception_types = []``, the ToolkitRegistry will attempt to call each ToolkitWrapper that provides the + requested method and if none succeeds, a single ``ValueError`` will be raised, with text listing the + errors that were raised by each ToolkitWrapper. - `PR #601 `_: Adds :py:meth:`RDKitToolkitWrapper.get_tagged_smarts_connectivity ` and @@ -158,7 +247,7 @@ New features py:meth:`Molecule.enumerate_stereoisomers `, py:meth:`Molecule.enumerate_protomers ` .. warning:: - Enumerate protomoers is currently only available through the OpenEye toolkit. + Enumerate protomers is currently only available through the OpenEye toolkit. - `PR #573 `_: Adds ``quacpac`` error output to ``quacpac`` failure in ``Molecule.compute_partial_charges_am1bcc``. - `PR #560 `_: Added visualization method to the the Molecule class. @@ -302,6 +391,7 @@ Bugfixes :py:meth:`Molecule.to_file ` can now correctly write multi-model PDB files when using the RDKit backend toolkit. + Examples added """""""""""""" - `PR #591 `_ and @@ -314,6 +404,7 @@ Examples added a QCArchive entry level record and calculate the energy using RDKit through QCEngine. + 0.6.0 - Library Charges ----------------------- diff --git a/docs/utils.rst b/docs/utils.rst index bf66ff3f1..f85cc4232 100644 --- a/docs/utils.rst +++ b/docs/utils.rst @@ -63,6 +63,7 @@ Miscellaneous utility functions. inherit_docstrings all_subclasses + temporary_cd get_data_file_path convert_0_1_smirnoff_to_0_2 convert_0_2_smirnoff_to_0_3 diff --git a/openforcefield/tests/test_forcefield.py b/openforcefield/tests/test_forcefield.py index e1cd487a8..c3419c1e4 100644 --- a/openforcefield/tests/test_forcefield.py +++ b/openforcefield/tests/test_forcefield.py @@ -24,11 +24,12 @@ import pytest from tempfile import NamedTemporaryFile -from openforcefield.utils.toolkits import OpenEyeToolkitWrapper, RDKitToolkitWrapper, AmberToolsToolkitWrapper, ToolkitRegistry +from openforcefield.utils.toolkits import (OpenEyeToolkitWrapper, RDKitToolkitWrapper, AmberToolsToolkitWrapper, + ToolkitRegistry, ChargeMethodUnavailableError) from openforcefield.utils import get_data_file_path from openforcefield.topology import Molecule, Topology -from openforcefield.typing.engines.smirnoff import ForceField, ParameterHandler, IncompatibleParameterError, SMIRNOFFSpecError -from openforcefield.typing.engines.smirnoff import XMLParameterIOHandler +from openforcefield.typing.engines.smirnoff import (ForceField, IncompatibleParameterError, SMIRNOFFSpecError, + XMLParameterIOHandler, ParameterHandler) #====================================================================== @@ -227,6 +228,25 @@ ''' +xml_spec_docs_charge_increment_model_xml = ''' + + + + + + + + + + +''' + +xml_charge_increment_model_formal_charges = ''' + + + +''' + xml_ff_torsion_bo = ''' @@ -344,7 +364,6 @@ def create_reversed_ethanol(): ethanol.partial_charges = charges return ethanol - def create_benzene_no_aromatic(): """ Creates an openforcefield.topology.Molecule representation of benzene through the API with aromatic bonds @@ -377,7 +396,6 @@ def create_benzene_no_aromatic(): benzene.add_bond(5, 11, 1, False) # C5 - C11 return benzene - def create_acetaldehyde(): """ Creates an openforcefield.topology.Molecule representation of acetaldehyde through the API @@ -400,6 +418,27 @@ def create_acetaldehyde(): acetaldehyde.partial_charges = charges return acetaldehyde +def create_acetate(): + """ + Creates an openforcefield.topology.Molecule representation of + acetate without the use of a cheminformatics toolkit + """ + # Create an acetate molecule without using a toolkit + acetate = Molecule() + acetate.add_atom(6, 0, False) # C0 + acetate.add_atom(6, 0, False) # C1 + acetate.add_atom(8, 0, False) # O2 + acetate.add_atom(8, -1, False) # O3 + acetate.add_atom(1, 0, False) # H4 + acetate.add_atom(1, 0, False) # H5 + acetate.add_atom(1, 0, False) # H6 + acetate.add_bond(0, 1, 1, False) # C0 - C1 + acetate.add_bond(1, 2, 2, False) # C1 = O2 + acetate.add_bond(1, 3, 1, False) # C1 - O3[-1] + acetate.add_bond(0, 4, 1, False) # C0 - H4 + acetate.add_bond(0, 5, 1, False) # C0 - H5 + acetate.add_bond(0, 6, 1, False) # C0 - H6 + return acetate def create_cyclohexane(): """ @@ -475,6 +514,59 @@ def create_cyclohexane(): 'omm_force': openmm.NonbondedForce.NoCutoff, 'exception': None, 'exception_match': ''}, ] +partial_charge_method_resolution_matrix = [ + {'toolkit': AmberToolsToolkitWrapper, + 'partial_charge_method': 'AM1-Mulliken', + 'exception': None, + 'exception_match': '' + }, + {'toolkit': AmberToolsToolkitWrapper, + 'partial_charge_method': 'Gasteiger', + 'exception': None, + 'exception_match': '' + }, + {'toolkit': AmberToolsToolkitWrapper, + 'partial_charge_method': 'Madeup-ChargeMethod', + 'exception': ChargeMethodUnavailableError, + 'exception_match': '' + }, + {'toolkit': OpenEyeToolkitWrapper, + 'partial_charge_method': 'AM1-Mulliken', + 'exception': None, + 'exception_match': '' + }, + {'toolkit': OpenEyeToolkitWrapper, + 'partial_charge_method': 'Gasteiger', + 'exception': None, + 'exception_match': '' + }, + {'toolkit': OpenEyeToolkitWrapper, + 'partial_charge_method': 'MMFF94', + 'exception': None, + 'exception_match': '' + }, + {'toolkit': OpenEyeToolkitWrapper, + 'partial_charge_method': 'am1bcc', + 'exception': None, + 'exception_match': '' + }, + {'toolkit': OpenEyeToolkitWrapper, + 'partial_charge_method': 'am1bccnosymspt', + 'exception': None, + 'exception_match': '' + }, + {'toolkit': OpenEyeToolkitWrapper, + 'partial_charge_method': 'am1bccelf10', + 'exception': None, + 'exception_match': '' + }, + {'toolkit': OpenEyeToolkitWrapper, + 'partial_charge_method': 'Madeup-ChargeMethod', + 'exception': ChargeMethodUnavailableError, + 'exception_match': '' + } + ] + #============================================================================================= # TESTS @@ -1289,6 +1381,214 @@ def test_parse_library_charges_from_spec_docs(self): ff = ForceField(xml_spec_docs_ala_library_charges_xml) ff = ForceField(xml_spec_docs_tip3p_library_charges_xml) + def test_parse_charge_increment_model_from_spec_docs(self): + """Ensure that the examples for librarycharges in the SMIRNOFF spec page are still valid""" + # TODO: This test is practically useless while the XML strings are hard-coded at the top of this file. + # We should implement something like doctests for the XML snippets on the SMIRNOFF spec page. + ff = ForceField(xml_spec_docs_charge_increment_model_xml) + + def test_charge_increment_model_forward_and_reverse_ethanol(self): + """Test application of ChargeIncrements to the same molecule with different orderings in the topology""" + test_charge_increment_model_ff = ''' + + + + + + + ''' + file_path = get_data_file_path('test_forcefields/smirnoff99Frosst.offxml') + ff = ForceField(file_path, test_charge_increment_model_ff) + del ff._parameter_handlers['ToolkitAM1BCC'] + top = Topology.from_molecules([create_ethanol(), create_reversed_ethanol()]) + sys = ff.create_openmm_system(top) + nonbonded_force = [force for force in sys.getForces() if isinstance(force, openmm.NonbondedForce)][0] + expected_charges = [0.2, -0.15, -0.05, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., -0.05, -0.15, 0.2] * unit.elementary_charge + for idx, expected_charge in enumerate(expected_charges): + charge, _, _ = nonbonded_force.getParticleParameters(idx) + assert abs(charge - expected_charge) < 1.e-6 * unit.elementary_charge + + def test_charge_increment_model_initialize_with_no_elements(self): + """Ensure that we can initialize a ForceField object from an OFFXML with a ChargeIncrementModel header, but no + ChargeIncrement elements""" + ff = ForceField(xml_charge_increment_model_formal_charges) + + def test_charge_increment_model_net_charge(self): + """Test application of charge increments on a molecule with a net charge""" + from simtk import unit + test_charge_increment_model_ff = ''' + + + + + + + ''' + file_path = get_data_file_path('test_forcefields/smirnoff99Frosst.offxml') + ff = ForceField(file_path, test_charge_increment_model_ff) + del ff._parameter_handlers['ToolkitAM1BCC'] + + acetate = create_acetate() + top = acetate.to_topology() + sys = ff.create_openmm_system(top) + nonbonded_force = [force for force in sys.getForces() if isinstance(force, openmm.NonbondedForce)][0] + expected_charges = [0, 0.15, -0.2 , -0.95, 0, 0, 0] * unit.elementary_charge + for idx, expected_charge in enumerate(expected_charges): + charge, _, _ = nonbonded_force.getParticleParameters(idx) + assert abs(charge - expected_charge) < 1.e-6 * unit.elementary_charge + + + def test_charge_increment_model_deduplicate_symmetric_matches(self): + """Test that chargeincrementmodelhandler deduplicates symmetric matches""" + from simtk import unit + + ethanol = create_ethanol() + top = ethanol.to_topology() + + # Test a charge increment that matches all C-H bonds at once + # (this should be applied once: C0-H3-H4-H5) + test_charge_increment_model_ff = ''' + + + + + + ''' + file_path = get_data_file_path('test_forcefields/smirnoff99Frosst.offxml') + ff = ForceField(file_path, test_charge_increment_model_ff) + del ff._parameter_handlers['ToolkitAM1BCC'] + + + sys = ff.create_openmm_system(top) + nonbonded_force = [force for force in sys.getForces() if isinstance(force, openmm.NonbondedForce)][0] + expected_charges = [0.3, 0, 0, -0.1, -0.1, -0.1, 0., 0., 0.] * unit.elementary_charge + for idx, expected_charge in enumerate(expected_charges): + charge, _, _ = nonbonded_force.getParticleParameters(idx) + assert abs(charge - expected_charge) < 1.e-6 * unit.elementary_charge + + # Test a charge increment that matches two C-H bonds at a time + # (this should be applied 3 times: C0-H3-H4, C0-H3-H5, C0-H4-H5) + test_charge_increment_model_ff = ''' + + + + + + ''' + file_path = get_data_file_path('test_forcefields/smirnoff99Frosst.offxml') + ff = ForceField(file_path, test_charge_increment_model_ff) + del ff._parameter_handlers['ToolkitAM1BCC'] + + + sys = ff.create_openmm_system(top) + nonbonded_force = [force for force in sys.getForces() if isinstance(force, openmm.NonbondedForce)][0] + expected_charges = [0.3, 0, 0, -0.1, -0.1, -0.1, 0., 0., 0.] * unit.elementary_charge + for idx, expected_charge in enumerate(expected_charges): + charge, _, _ = nonbonded_force.getParticleParameters(idx) + assert abs(charge - expected_charge) < 1.e-6 * unit.elementary_charge + + + # Test a charge increment that matches ONE C-H bond at a time + # (this should be applied three times: C0-H3, C0-H4, C0-H5) + test_charge_increment_model_ff = ''' + + + + + + ''' + file_path = get_data_file_path('test_forcefields/smirnoff99Frosst.offxml') + ff = ForceField(file_path, test_charge_increment_model_ff) + del ff._parameter_handlers['ToolkitAM1BCC'] + + sys = ff.create_openmm_system(top) + nonbonded_force = [force for force in sys.getForces() if isinstance(force, openmm.NonbondedForce)][0] + expected_charges = [0.3, 0, 0, -0.1, -0.1, -0.1, 0., 0., 0.] * unit.elementary_charge + for idx, expected_charge in enumerate(expected_charges): + charge, _, _ = nonbonded_force.getParticleParameters(idx) + assert abs(charge - expected_charge) < 1.e-6 * unit.elementary_charge + + def test_charge_increment_model_completely_overlapping_matches_override(self): + """Ensure that DIFFERENT chargeincrements override one another if they apply to the + same atoms, regardless of order""" + from simtk import unit + test_charge_increment_model_ff = ''' + + + + + + + ''' + file_path = get_data_file_path('test_forcefields/smirnoff99Frosst.offxml') + ff = ForceField(file_path, test_charge_increment_model_ff) + del ff._parameter_handlers['ToolkitAM1BCC'] + + ethanol = create_ethanol() + top = ethanol.to_topology() + sys = ff.create_openmm_system(top) + nonbonded_force = [force for force in sys.getForces() if isinstance(force, openmm.NonbondedForce)][0] + expected_charges = [0.3, 0, 0, -0.1, -0.1, -0.1, 0., 0., 0.] * unit.elementary_charge + for idx, expected_charge in enumerate(expected_charges): + charge, _, _ = nonbonded_force.getParticleParameters(idx) + assert abs(charge - expected_charge) < 1.e-6 * unit.elementary_charge + + + def test_charge_increment_model_partially_overlapping_matches_both_apply(self): + """Ensure that DIFFERENT chargeincrements BOTH get applied if they match + a partially-overlapping set of atoms""" + from simtk import unit + test_charge_increment_model_ff = ''' + + + + + + + ''' + file_path = get_data_file_path('test_forcefields/smirnoff99Frosst.offxml') + ff = ForceField(file_path, test_charge_increment_model_ff) + del ff._parameter_handlers['ToolkitAM1BCC'] + + ethanol = create_ethanol() + top = ethanol.to_topology() + sys = ff.create_openmm_system(top) + nonbonded_force = [force for force in sys.getForces() if isinstance(force, openmm.NonbondedForce)][0] + expected_charges = [0.35, -0.05, 0, -0.1, -0.1, -0.1, 0., 0., 0.] * unit.elementary_charge + for idx, expected_charge in enumerate(expected_charges): + charge, _, _ = nonbonded_force.getParticleParameters(idx) + assert abs(charge - expected_charge) < 1.e-6 * unit.elementary_charge + + @pytest.mark.parametrize("inputs", partial_charge_method_resolution_matrix) + def test_partial_charge_resolution(self, inputs): + """Check that the proper partial charge methods are available, and that unavailable partial charge methods + raise an exception. + """ + toolkit_wrapper_class = inputs['toolkit'] + if not(toolkit_wrapper_class.is_available()): + pytest.skip(f"{toolkit_wrapper_class} is not available.") + toolkit_wrapper = toolkit_wrapper_class() + partial_charge_method = inputs['partial_charge_method'] + expected_exception = inputs['exception'] + expected_exception_match = inputs['exception_match'] + ethanol = create_ethanol() + ethanol.generate_conformers() + if expected_exception is None: + ethanol.assign_partial_charges(partial_charge_method=partial_charge_method, + toolkit_registry=toolkit_wrapper) + abs_charge_sum = 0. * unit.elementary_charge + + # Ensure that nonzero charges were assigned + for pc in ethanol.partial_charges: + abs_charge_sum += abs(pc) + assert abs_charge_sum > 0.5 * unit.elementary_charge + + else: + with pytest.raises(expected_exception, match=expected_exception_match) as excinfo: + ethanol.assign_partial_charges(partial_charge_method=partial_charge_method, + toolkit_registry=toolkit_wrapper) + def test_library_charge_hierarchy(self): """Test assigning charges to one water molecule using library charges, where two LCs match and the assignment is determined by order they are added to the force field""" @@ -1329,8 +1629,8 @@ def test_library_charges_to_two_waters(self): q, sigma, epsilon = nonbondedForce.getParticleParameters(particle_index) assert q == expected_charge - def test_library_charges_to_two_ethanols_different_atom_ordering(self): - """Test assigning charges to two ethanols with different atom orderings""" + def test_library_charges_to_three_ethanols_different_atom_ordering(self): + """Test assigning charges to three ethanols with different atom orderings""" from simtk.openmm import NonbondedForce # Define a library charge parameter for ethanol (C1-C2-O3) where C1 has charge -0.2, and its Hs have -0.02, @@ -1351,18 +1651,29 @@ def test_library_charges_to_two_ethanols_different_atom_ordering(self): # H6 - C1 - C3 - O2 - H4 # | | # H7 H9 + # + # create_reversed_ethanol() + # H5 H2 + # | | + # H4 - C8 - C7 - O6 - H0 + # | | + # H3 H1 + molecules = [Molecule.from_file(get_data_file_path('molecules/ethanol.sdf')), - Molecule.from_file(get_data_file_path('molecules/ethanol_reordered.sdf'))] + Molecule.from_file(get_data_file_path('molecules/ethanol_reordered.sdf')), + create_reversed_ethanol()] top = Topology.from_molecules(molecules) omm_system = ff.create_openmm_system(top) nonbondedForce = [f for f in omm_system.getForces() if type(f) == NonbondedForce][0] - expected_charges = [-0.2, -0.1, 0.3, 0.08, -0.02, -0.02, -0.02, -0.01, -0.01, -0.2, - 0.3, -0.1, 0.08, -0.02, -0.02, -0.02, -0.01, -0.01] * unit.elementary_charge + expected_charges = [-0.2, -0.1, 0.3, 0.08, -0.02, -0.02, -0.02, -0.01, -0.01, + -0.2, 0.3, -0.1, 0.08, -0.02, -0.02, -0.02, -0.01, -0.01, + 0.08, -0.01, -0.01, -0.02, -0.02, -0.02, 0.3, -0.1, -0.2] * unit.elementary_charge for particle_index, expected_charge in enumerate(expected_charges): q, sigma, epsilon = nonbondedForce.getParticleParameters(particle_index) assert q == expected_charge + @pytest.mark.parametrize('monatomic_ion,formal_charge', generate_monatomic_ions()) def test_library_charges_monatomic_ions(self, monatomic_ion, formal_charge): """Test assigning library charges to each of the monatomic ions in openff-1.1.0.xml""" @@ -1384,13 +1695,23 @@ def test_charge_method_hierarchy(self): ff = ForceField('test_forcefields/smirnoff99Frosst.offxml', xml_CH_zeroes_library_charges_xml, - 'test_forcefields/tip3p.offxml' + 'test_forcefields/tip3p.offxml', + xml_charge_increment_model_formal_charges ) + # Cyclohexane will be assigned nonzero charges based on `charge_from_molecules` kwarg cyclohexane = Molecule.from_file(get_data_file_path(os.path.join('systems', 'monomers','cyclohexane.sdf'))) + # Butanol will be assigned zero-valued charges based on `charge_from_molecules` kwarg butanol = Molecule.from_file(get_data_file_path(os.path.join('systems', 'monomers', 'butanol.sdf'))) + # Propane will be assigned charges from AM1BCC propane = Molecule.from_file(get_data_file_path(os.path.join('systems', 'monomers', 'propane.sdf'))) + # Water will be assigned TIP3P librarycharges water = Molecule.from_file(get_data_file_path(os.path.join('systems', 'monomers', 'water.sdf'))) + # Ethanol should be caught by ToolkitAM1BCC ethanol = Molecule.from_file(get_data_file_path(os.path.join('systems', 'monomers', 'ethanol.sdf'))) + # iodide should fail to be assigned charges by AM1-BCC, and instead be caught by ChargeIncrementHandler + # (and assigned formal charge by dummy charge method) + iodide = Molecule.from_smiles('[I-1]') + # Assign dummy partial charges to cyclohexane, which we expect to find in the final system since it # is included in the charge_from_molecules kwarg to create_openmm_system cyclohexane.partial_charges = np.array([-0.2, -0.2, -0.2, -0.2, -0.2, -0.2, @@ -1415,7 +1736,8 @@ def test_charge_method_hierarchy(self): butanol, # charge_from_molecules kwarg propane, # library charges water, # library charges - ethanol] # AM1-BCC + ethanol, # AM1-BCC + iodide] # charge increment model (formal charge) top = Topology.from_molecules(molecules) omm_system = ff.create_openmm_system(top, charge_from_molecules=[cyclohexane, butanol]) existing = [f for f in omm_system.getForces() if type(f) == NonbondedForce] @@ -1437,16 +1759,21 @@ def test_charge_method_hierarchy(self): # water (3 atoms) should have the following charges from tip3p.offxml -0.834, 0.417, 0.417] * unit.elementary_charge - # Ensure that the first three molecules have exactly the charges we intended + # Ensure that the first four molecules have exactly the charges we intended for particle_index, expected_charge in enumerate(expected_charges): q, sigma, epsilon = nonbondedForce.getParticleParameters(particle_index) assert q == expected_charge # Ensure the last molecule (ethanol) had _some_ nonzero charge assigned by an AM1BCC implementation - for particle_index in range(len(expected_charges), top.n_topology_atoms): + for particle_index in range(len(expected_charges), top.n_topology_atoms-1): q, sigma, epsilon = nonbondedForce.getParticleParameters(particle_index) assert q != 0 * unit.elementary_charge + # Ensure that iodine has a charge of -1, specified by charge increment model charge_method="formal charge" + q, sigma, epsilon = nonbondedForce.getParticleParameters(top.n_topology_atoms-1) + assert q == -1. * unit.elementary_charge + + def test_assign_charges_to_molecule_in_parts_using_multiple_library_charges(self): """Test assigning charges to parts of a molecule using two library charge lines. Note that these LibraryCharge SMIRKS have partial overlap, so this also tests that the hierarchy is correctly obeyed.""" @@ -2348,7 +2675,7 @@ def test_electrostatics_options(self): """ molecules_file_path = get_data_file_path('molecules/laromustine_tripos.mol2') molecule = openforcefield.topology.Molecule.from_file(molecules_file_path) - forcefield = ForceField([smirnoff99Frosst_offxml_file_path, chargeincrement_offxml_file_path]) + forcefield = ForceField([smirnoff99Frosst_offxml_file_path, charge_increment_offxml_file_path]) for method in ['PME', 'reaction-field', 'Coulomb']: # Change electrostatics method forcefield.forces['Electrostatics'].method = method @@ -2359,7 +2686,7 @@ def test_electrostatics_options(self): # AMBER-parameterized system to OFF-parameterized systems @pytest.mark.skip(reason='Needs to be updated for 0.2.0 syntax') -def test_chargeincrement(self): +def test_charge_increment(self): """Test parameter assignment using smirnoff99Frosst on laromustine with ChargeIncrementModel. """ molecules_file_path = get_data_file_path('molecules/laromustine_tripos.mol2') diff --git a/openforcefield/tests/test_molecule.py b/openforcefield/tests/test_molecule.py index 318f9381f..5ea9aaad9 100644 --- a/openforcefield/tests/test_molecule.py +++ b/openforcefield/tests/test_molecule.py @@ -227,16 +227,21 @@ def test_atom_constructor(self): # Create a non-aromatic carbon atom atom1 = Atom(6, 0, False) assert atom1.atomic_number == 6 - assert atom1.formal_charge == 0 + assert atom1.formal_charge == 0 * unit.elementary_charge # Create a chiral carbon atom atom2 = Atom(6, 0, False, stereochemistry='R', name='CT') assert atom1.stereochemistry != atom2.stereochemistry + # Ensure that formal charge can also be set as a Quantity + atom1 = Atom(6, 1*unit.elementary_charge, False) + assert atom1.formal_charge == 1 * unit.elementary_charge + + def test_atom_properties(self): """Test that atom properties are correctly populated and gettable""" from simtk.openmm.app import element - formal_charge = 0 + formal_charge = 0 * unit.elementary_charge is_aromatic = False # Attempt to create all elements supported by OpenMM elements = [getattr(element, name) for name in dir(element) if (type(getattr(element, name)) == element.Element)] @@ -1547,8 +1552,10 @@ def test_torsions(self, molecule): @pytest.mark.parametrize('molecule', mini_drug_bank()) def test_total_charge(self, molecule): """Test total charge""" - total_charge = sum([atom.formal_charge for atom in molecule.atoms]) - assert total_charge == molecule.total_charge + charge_sum = 0 * unit.elementary_charge + for atom in molecule.atoms: + charge_sum += atom.formal_charge + assert charge_sum == molecule.total_charge # ---------------------------------------------------- # Test magic methods. diff --git a/openforcefield/tests/test_parameters.py b/openforcefield/tests/test_parameters.py index e575545dc..6986dba96 100644 --- a/openforcefield/tests/test_parameters.py +++ b/openforcefield/tests/test_parameters.py @@ -23,7 +23,7 @@ ParameterAttribute, IndexedParameterAttribute, ParameterList, ParameterType, BondHandler, ParameterHandler, ProperTorsionHandler, ImproperTorsionHandler, LibraryChargeHandler, GBSAHandler, SMIRNOFFSpecError, - _ParameterAttributeHandler + _ParameterAttributeHandler, ChargeIncrementModelHandler, IncompatibleParameterError ) from openforcefield.utils import detach_units, IncompatibleUnitError from openforcefield.utils.collections import ValidatedList @@ -1207,6 +1207,101 @@ def test_create_library_charge_handler(self): """Test creation of an empty LibraryChargeHandler""" handler = LibraryChargeHandler(skip_version_check=True) + def test_library_charge_type_wrong_num_charges(self): + """Ensure that an error is raised if a LibraryChargeType is initialized with a different number of + tagged atoms and charges""" + lc_type = LibraryChargeHandler.LibraryChargeType(smirks='[#6:1]-[#7:2]', + charge1=0.1 * unit.elementary_charge, + charge2=-0.1 * unit.elementary_charge) + + lc_type = LibraryChargeHandler.LibraryChargeType(smirks='[#6:1]-[#7:2]-[#6]', + charge1=0.1 * unit.elementary_charge, + charge2=-0.1 * unit.elementary_charge) + + with pytest.raises(SMIRNOFFSpecError, match="initialized with unequal number of tagged atoms and charges") as excinfo: + lc_type = LibraryChargeHandler.LibraryChargeType(smirks='[#6:1]-[#7:2]', + charge1=0.05 * unit.elementary_charge, + charge2=0.05 * unit.elementary_charge, + charge3=-0.1 * unit.elementary_charge) + + with pytest.raises(SMIRNOFFSpecError, match="initialized with unequal number of tagged atoms and charges") as excinfo: + lc_type = LibraryChargeHandler.LibraryChargeType(smirks='[#6:1]-[#7:2]-[#6]', + charge1=0.05 * unit.elementary_charge, + charge2=0.05 * unit.elementary_charge, + charge3=-0.1 * unit.elementary_charge) + + with pytest.raises(SMIRNOFFSpecError, match="initialized with unequal number of tagged atoms and charges") as excinfo: + lc_type = LibraryChargeHandler.LibraryChargeType(smirks='[#6:1]-[#7:2]-[#6]', + charge1=0.05 * unit.elementary_charge) + +class TestChargeIncrementModelHandler: + def test_create_charge_increment_model_handler(self): + """Test creation of ChargeIncrementModelHandlers""" + handler = ChargeIncrementModelHandler(skip_version_check=True) + assert handler.number_of_conformers == 1 + assert handler.partial_charge_method == 'AM1-Mulliken' + handler = ChargeIncrementModelHandler(skip_version_check=True, number_of_conformers=10) + handler = ChargeIncrementModelHandler(skip_version_check=True, number_of_conformers=1) + handler = ChargeIncrementModelHandler(skip_version_check=True, number_of_conformers="10") + handler = ChargeIncrementModelHandler(skip_version_check=True, number_of_conformers=0) + handler = ChargeIncrementModelHandler(skip_version_check=True, number_of_conformers="0") + with pytest.raises(TypeError) as excinfo: + handler = ChargeIncrementModelHandler(skip_version_check=True, number_of_conformers=None) + with pytest.raises(SMIRNOFFSpecError) as excinfo: + handler = ChargeIncrementModelHandler(skip_version_check=True, n_conformers=[10]) + handler = ChargeIncrementModelHandler(skip_version_check=True, partial_charge_method="AM1-Mulliken") + handler = ChargeIncrementModelHandler(skip_version_check=True, partial_charge_method="Gasteiger") + handler = ChargeIncrementModelHandler(skip_version_check=True, partial_charge_method=None) + + def test_charge_increment_model_handler_getters_setters(self): + """Test ChargeIncrementModelHandler getters and setters""" + handler = ChargeIncrementModelHandler(skip_version_check=True) + assert handler.number_of_conformers == 1 + assert handler.partial_charge_method == 'AM1-Mulliken' + handler.number_of_conformers = 2 + assert handler.number_of_conformers == 2 + handler.number_of_conformers = "3" + assert handler.number_of_conformers == 3 + with pytest.raises(ValueError) as excinfo: + handler.number_of_conformers = "string that can't be cast to int" + + def test_charge_increment_model_handlers_are_compatible(self): + """Test creation of ChargeIncrementModelHandlers""" + handler1 = ChargeIncrementModelHandler(skip_version_check=True) + handler2 = ChargeIncrementModelHandler(skip_version_check=True) + handler1.check_handler_compatibility(handler2) + + handler3 = ChargeIncrementModelHandler(skip_version_check=True, number_of_conformers='9') + with pytest.raises(IncompatibleParameterError) as excinfo: + handler1.check_handler_compatibility(handler3) + + def test_charge_increment_type_wrong_num_increments(self): + """Ensure that an error is raised if a ChargeIncrementType is initialized with a different number of + tagged atoms and chargeincrements""" + ci_type = ChargeIncrementModelHandler.ChargeIncrementType(smirks='[#6:1]-[#7:2]', + charge_increment1=0.1 * unit.elementary_charge, + charge_increment2=-0.1 * unit.elementary_charge) + + ci_type = ChargeIncrementModelHandler.ChargeIncrementType(smirks='[#6:1]-[#7:2]-[#6]', + charge_increment1=0.1 * unit.elementary_charge, + charge_increment2=-0.1 * unit.elementary_charge) + + with pytest.raises(SMIRNOFFSpecError, match="initialized with unequal number of tagged atoms and charge increments") as excinfo: + ci_type = ChargeIncrementModelHandler.ChargeIncrementType(smirks='[#6:1]-[#7:2]', + charge_increment1=0.05 * unit.elementary_charge, + charge_increment2=0.05 * unit.elementary_charge, + charge_increment3=-0.1 * unit.elementary_charge) + + with pytest.raises(SMIRNOFFSpecError, match="initialized with unequal number of tagged atoms and charge increments") as excinfo: + ci_type = ChargeIncrementModelHandler.ChargeIncrementType(smirks='[#6:1]-[#7:2]-[#6]', + charge_increment1=0.05 * unit.elementary_charge, + charge_increment2=0.05 * unit.elementary_charge, + charge_increment3=-0.1 * unit.elementary_charge) + + with pytest.raises(SMIRNOFFSpecError, match="initialized with unequal number of tagged atoms and charge increments") as excinfo: + ci_type = ChargeIncrementModelHandler.ChargeIncrementType(smirks='[#6:1]-[#7:2]-[#6]', + charge_increment1=0.05 * unit.elementary_charge) + class TestGBSAHandler: def test_create_default_gbsahandler(self): @@ -1263,7 +1358,6 @@ def test_gbsahandlers_are_compatible(self): """ Test the check_handler_compatibility function of GBSAHandler """ - from openforcefield.typing.engines.smirnoff import IncompatibleParameterError from simtk import unit gbsa_handler_1 = GBSAHandler(skip_version_check=True) gbsa_handler_2 = GBSAHandler(skip_version_check=True) diff --git a/openforcefield/tests/test_toolkits.py b/openforcefield/tests/test_toolkits.py index e8be9b0b1..25f905124 100644 --- a/openforcefield/tests/test_toolkits.py +++ b/openforcefield/tests/test_toolkits.py @@ -21,12 +21,15 @@ import pytest from openforcefield.utils.toolkits import (OpenEyeToolkitWrapper, RDKitToolkitWrapper, - AmberToolsToolkitWrapper, ToolkitRegistry, - GAFFAtomTypeWarning, UndefinedStereochemistryError) + AmberToolsToolkitWrapper, BuiltInToolkitWrapper, ToolkitRegistry, + ToolkitWrapper, + GAFFAtomTypeWarning, UndefinedStereochemistryError, + ChargeMethodUnavailableError, IncorrectNumConformersError, + IncorrectNumConformersWarning) from openforcefield.utils import get_data_file_path from openforcefield.topology.molecule import Molecule from openforcefield.tests.test_forcefield import create_ethanol, create_cyclohexane, create_acetaldehyde, \ - create_reversed_ethanol + create_reversed_ethanol, create_acetate #============================================================================================= # FIXTURES @@ -657,77 +660,202 @@ def test_generate_conformers(self): assert molecule.n_conformers != 0 assert not(molecule.conformers[0] == (0.*unit.angstrom)).all() - # TODO: Make this test more robust - @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') - def test_compute_partial_charges(self): - """Test OpenEyeToolkitWrapper compute_partial_charges()""" + def test_generate_multiple_conformers(self): + """Test OpenEyeToolkitWrapper generate_conformers() for generating multiple conformers""" toolkit_wrapper = OpenEyeToolkitWrapper() - smiles = '[H]C([H])([H])C([H])([H])[H]' + smiles = 'CCCCCCC' molecule = toolkit_wrapper.from_smiles(smiles) - # Ensure that an exception is raised if no conformers are provided - with pytest.raises(Exception) as excinfo: - molecule.compute_partial_charges(toolkit_registry=toolkit_wrapper) - molecule.generate_conformers(toolkit_registry=toolkit_wrapper) - # Ensure that an exception is raised if an invalid charge model is passed in - with pytest.raises(Exception) as excinfo: - charge_model = 'notARealChargeModel' - molecule.compute_partial_charges(toolkit_registry=toolkit_wrapper, charge_model=charge_model) - - # TODO: Test all supported charge models - # Note: "amber" and "amberff94" only work for a subset of residue types, so we'll need to find testing data for - # those - # charge_model = [,'amber','amberff94'] - # TODO: 'mmff' and 'mmff94' often assign charges of 0, presumably if the molecule is unrecognized. - # charge_model = ['mmff', 'mmff94'] - for charge_model in ['noop', 'am1bcc', 'am1bccnosymspt', 'am1bccelf10']: - with pytest.raises(NotImplementedError) as excinfo: - molecule.compute_partial_charges(toolkit_registry=toolkit_wrapper) # , charge_model=charge_model) - charge_sum = 0 * unit.elementary_charge - for pc in molecule._partial_charges: - charge_sum += pc - assert charge_sum < 0.001 * unit.elementary_charge - - # For now, just test AM1-BCC while the SMIRNOFF spec for other charge models gets worked out - molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_wrapper) # , charge_model=charge_model) + molecule.generate_conformers(rms_cutoff=1*unit.angstrom, + n_conformers=100, + toolkit_registry=toolkit_wrapper) + assert molecule.n_conformers > 1 + assert not(molecule.conformers[0] == (0.*unit.angstrom)).all() + + # Ensure rms_cutoff kwarg is working + molecule2 = toolkit_wrapper.from_smiles(smiles) + molecule2.generate_conformers(rms_cutoff=0.1*unit.angstrom, + n_conformers=100, + toolkit_registry=toolkit_wrapper) + assert molecule2.n_conformers > molecule.n_conformers + + # Ensure n_conformers kwarg is working + molecule2 = toolkit_wrapper.from_smiles(smiles) + molecule2.generate_conformers(rms_cutoff=0.1*unit.angstrom, + n_conformers=10, + toolkit_registry=toolkit_wrapper) + assert molecule2.n_conformers == 10 + + @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') + def test_compute_partial_charges_am1bcc(self): + """Test OpenEyeToolkitWrapper compute_partial_charges_am1bcc()""" + toolkit_registry = ToolkitRegistry(toolkit_precedence=[OpenEyeToolkitWrapper]) + molecule = create_ethanol() + molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry) # , charge_model=charge_model) + charge_sum = 0 * unit.elementary_charge + abs_charge_sum = 0 * unit.elementary_charge + for pc in molecule._partial_charges: + charge_sum += pc + abs_charge_sum += abs(pc) + assert abs(charge_sum) < 0.005 * unit.elementary_charge + assert abs_charge_sum > 0.25 * unit.elementary_charge + + @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') + def test_compute_partial_charges_am1bcc_net_charge(self): + """Test OpenEyeToolkitWrapper assign_partial_charges() on a molecule with a net +1 charge""" + toolkit_registry = ToolkitRegistry(toolkit_precedence=[OpenEyeToolkitWrapper]) + molecule = create_acetate() + molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry) charge_sum = 0 * unit.elementary_charge for pc in molecule._partial_charges: charge_sum += pc - assert charge_sum < 0.001 * unit.elementary_charge + assert -0.999 * unit.elementary_charge > charge_sum > -1.001 * unit.elementary_charge @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') - def test_compute_partial_charges_net_charge(self): - """Test OpenEyeToolkitWrapper compute_partial_charges() on a molecule with a net +1 charge""" + def test_compute_partial_charges_am1bcc_wrong_n_confs(self): + """ + Test OpenEyeToolkitWrapper compute_partial_charges_am1bcc() when requesting to use an incorrect number of + conformers. This test is a bit shorter than that for AmberToolsToolkitWrapper because OETK uses the + ELF10 multiconformer method of AM1BCC, which doesn't have a maximum number of conformers. + """ + from openforcefield.tests.test_forcefield import create_ethanol + toolkit_registry = ToolkitRegistry(toolkit_precedence=[OpenEyeToolkitWrapper]) + molecule = create_ethanol() + molecule.generate_conformers(n_conformers=2, + rms_cutoff=0.1*unit.angstrom, + toolkit_registry=toolkit_registry) - toolkit_wrapper = OpenEyeToolkitWrapper() - smiles = '[H]C([H])([H])[N+]([H])([H])[H]' - molecule = toolkit_wrapper.from_smiles(smiles) - molecule.generate_conformers(toolkit_registry=toolkit_wrapper) + # Try again, with strict_n_confs as true, but not including use_confs, so the + # recommended number of confs will be generated + molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry, + strict_n_conformers=True) - with pytest.raises(NotImplementedError) as excinfo: - charge_model = 'notARealChargeModel' - molecule.compute_partial_charges(toolkit_registry=toolkit_wrapper)#, charge_model=charge_model) - - # TODO: Test all supported charge models - # TODO: "amber" and "amberff94" only work for a subset of residue types, so we'll need to find testing data for - # those - # charge_model = [,'amber','amberff94'] - # The 'noop' charge model doesn't add up to the formal charge, so we shouldn't test it - # charge_model = ['noop'] - for charge_model in ['mmff', 'mmff94', 'am1bcc', 'am1bccnosymspt', 'am1bccelf10']: - with pytest.raises(NotImplementedError) as excinfo: - molecule.compute_partial_charges(toolkit_registry=toolkit_wrapper) #, charge_model=charge_model) - charge_sum = 0 * unit.elementary_charge - for pc in molecule._partial_charges: - charge_sum += pc - assert 0.999 * unit.elementary_charge < charge_sum < 1.001 * unit.elementary_charge - # For now, I'm just testing AM1-BCC (will test more when the SMIRNOFF spec for other charges is finalized) - molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_wrapper) - charge_sum = 0 * unit.elementary_charge - for pc in molecule._partial_charges: + + @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') + @pytest.mark.parametrize("partial_charge_method", ['am1bcc', 'am1-mulliken', 'gasteiger']) + def test_assign_partial_charges_neutral(self, partial_charge_method): + """Test OpenEyeToolkitWrapper assign_partial_charges()""" + from openforcefield.tests.test_forcefield import create_ethanol + toolkit_registry = ToolkitRegistry(toolkit_precedence=[OpenEyeToolkitWrapper]) + molecule = create_ethanol() + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method) + charge_sum = 0. * unit.elementary_charge + for pc in molecule.partial_charges: + charge_sum += pc + assert -1.e-5 < charge_sum.value_in_unit(unit.elementary_charge) < 1.e-5 + + @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') + @pytest.mark.parametrize("partial_charge_method", ['am1bcc', 'am1-mulliken']) + def test_assign_partial_charges_conformer_dependence(self, partial_charge_method): + """Test OpenEyeToolkitWrapper assign_partial_charges()'s use_conformers kwarg + to ensure charges are really conformer dependent. Skip Gasteiger because it isn't + conformer dependent.""" + from openforcefield.tests.test_forcefield import create_ethanol + import copy + toolkit_registry = ToolkitRegistry(toolkit_precedence=[OpenEyeToolkitWrapper]) + molecule = create_ethanol() + molecule.generate_conformers(n_conformers=1) + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method, + use_conformers=molecule.conformers) + pcs1 = copy.deepcopy(molecule.partial_charges) + molecule._conformers[0][0][0] += 0.2 * unit.angstrom + molecule._conformers[0][1][1] -= 0.2 * unit.angstrom + molecule._conformers[0][2][1] += 0.2 * unit.angstrom + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method, + use_conformers=molecule.conformers) + for pc1, pc2 in zip(pcs1, molecule.partial_charges): + assert abs(pc1 - pc2) > 1.e-5 * unit.elementary_charge + + @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') + @pytest.mark.parametrize("partial_charge_method", ['am1bcc', 'am1-mulliken', 'gasteiger']) + def test_assign_partial_charges_net_charge(self, partial_charge_method): + """ + Test OpenEyeToolkitWrapper assign_partial_charges() on a molecule with net charge. + """ + from openforcefield.tests.test_forcefield import create_acetate + toolkit_registry = ToolkitRegistry(toolkit_precedence=[OpenEyeToolkitWrapper]) + molecule = create_acetate() + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method) + charge_sum = 0. * unit.elementary_charge + for pc in molecule.partial_charges: charge_sum += pc - assert 0.999 * unit.elementary_charge < charge_sum < 1.001 * unit.elementary_charge + assert -1.e-5 < charge_sum.value_in_unit(unit.elementary_charge) + 1. < 1.e-5 + + @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') + def test_assign_partial_charges_bad_charge_method(self): + """Test OpenEyeToolkitWrapper assign_partial_charges() for a nonexistent charge method""" + from openforcefield.tests.test_forcefield import create_ethanol + toolkit_registry = ToolkitRegistry(toolkit_precedence=[OpenEyeToolkitWrapper]) + molecule = create_ethanol() + + # Molecule.assign_partial_charges calls the ToolkitRegistry with raise_exception_types = [], + # which means it will only ever return ValueError + with pytest.raises(ValueError, match="is not available from OpenEyeToolkitWrapper") as excinfo: + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method="NotARealChargeMethod") + + # ToolkitWrappers raise a specific exception class, so we test that here + with pytest.raises(ChargeMethodUnavailableError, match="is not available from OpenEyeToolkitWrapper") as excinfo: + OETKW = OpenEyeToolkitWrapper() + OETKW.assign_partial_charges(molecule=molecule, + partial_charge_method="NotARealChargeMethod") + + @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') + @pytest.mark.parametrize("partial_charge_method,expected_n_confs", [('am1bcc', 1), + ('am1-mulliken', 1), + ('gasteiger', 0)]) + def test_assign_partial_charges_wrong_n_confs(self, partial_charge_method, expected_n_confs): + """ + Test OpenEyeToolkitWrapper assign_partial_charges() when requesting to use an incorrect number of + conformers + """ + from openforcefield.tests.test_forcefield import create_ethanol + toolkit_registry = ToolkitRegistry(toolkit_precedence=[OpenEyeToolkitWrapper]) + molecule = create_ethanol() + molecule.generate_conformers(n_conformers=2, rms_cutoff=0.01*unit.angstrom) + + # Try passing in the incorrect number of confs, but without specifying strict_n_conformers, + # which should produce a warning + with pytest.warns(IncorrectNumConformersWarning, + match=f"has 2 conformers, but charge method '{partial_charge_method}' " + f"expects exactly {expected_n_confs}."): + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method, + use_conformers=molecule.conformers, + strict_n_conformers=False) + + # Try again, with strict_n_confs as true, but not including use_confs, so the + # recommended number of confs will be generated + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method, + strict_n_conformers=True) + + # Test calling the ToolkitWrapper _indirectly_, though a ToolkitRegistry, + # which should aggregate any exceptions and bundle all of the messages + # in a failed task together in a single ValueError. + with pytest.raises(ValueError, + match=f"has 2 conformers, but charge method '{partial_charge_method}' " + f"expects exactly {expected_n_confs}."): + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method, + use_conformers=molecule.conformers, + strict_n_conformers=True) + + # Test calling the ToolkitWrapper _directly_, passing in the incorrect number of + # confs, and specify strict_n_conformers, which should produce an IncorrectNumConformersError + with pytest.raises(IncorrectNumConformersError, + match=f"has 2 conformers, but charge method '{partial_charge_method}' " + f"expects exactly {expected_n_confs}."): + OETKW = OpenEyeToolkitWrapper() + OETKW.assign_partial_charges(molecule=molecule, + partial_charge_method=partial_charge_method, + use_conformers=molecule.conformers, + strict_n_conformers=True) @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') def test_compute_partial_charges_failure(self): @@ -892,7 +1020,7 @@ def test_substructure_search_on_large_molecule(self): molecule = tk.from_smiles(smiles) query = "[C:1]~[C:2]" ret = molecule.chemical_environment_matches(query, toolkit_registry=tk) - assert len(ret) == 1198 + assert len(ret) == 1198 assert len(ret[0]) == 2 def test_find_rotatable_bonds(self): @@ -1475,6 +1603,32 @@ def test_generate_conformers(self): molecule.generate_conformers() # TODO: Make this test more robust + @pytest.mark.skipif(not RDKitToolkitWrapper.is_available(), reason='RDKit Toolkit not available') + def test_generate_multiple_conformers(self): + """Test RDKitToolkitWrapper generate_conformers() for generating multiple conformers""" + toolkit_wrapper = RDKitToolkitWrapper() + smiles = 'CCCCCCC' + molecule = toolkit_wrapper.from_smiles(smiles) + molecule.generate_conformers(rms_cutoff=1 * unit.angstrom, + n_conformers=100, + toolkit_registry=toolkit_wrapper) + assert molecule.n_conformers > 1 + assert not (molecule.conformers[0] == (0. * unit.angstrom)).all() + + # Ensure rms_cutoff kwarg is working + molecule2 = toolkit_wrapper.from_smiles(smiles) + molecule2.generate_conformers(rms_cutoff=0.1 * unit.angstrom, + n_conformers=100, + toolkit_registry=toolkit_wrapper) + assert molecule2.n_conformers > molecule.n_conformers + + # Ensure n_conformers kwarg is working + molecule2 = toolkit_wrapper.from_smiles(smiles) + molecule2.generate_conformers(rms_cutoff=0.1 * unit.angstrom, + n_conformers=10, + toolkit_registry=toolkit_wrapper) + assert molecule2.n_conformers == 10 + @pytest.mark.skipif(not RDKitToolkitWrapper.is_available(), reason='RDKit Toolkit not available') def test_find_rotatable_bonds(self): """Test finding rotatable bonds while ignoring some groups""" @@ -1557,7 +1711,6 @@ def test_to_rdkit_losing_aromaticity_(self): assert offatom.is_aromatic is rdatom.GetIsAromatic() - # TODO: Add test for higher bonds orders # TODO: Add test for aromaticity # TODO: Add test and molecule functionality for isotopes @@ -1569,72 +1722,226 @@ def test_to_rdkit_losing_aromaticity_(self): - + class TestAmberToolsToolkitWrapper: """Test the AmberToolsToolkitWrapper""" @pytest.mark.skipif(not RDKitToolkitWrapper.is_available() or not AmberToolsToolkitWrapper.is_available(), reason='RDKitToolkit and AmberToolsToolkit not available') - def test_compute_partial_charges(self): - """Test OpenEyeToolkitWrapper compute_partial_charges()""" + def test_compute_partial_charges_am1bcc(self): + """Test AmberToolsToolkitWrapper compute_partial_charges_am1bcc()""" toolkit_registry = ToolkitRegistry(toolkit_precedence=[AmberToolsToolkitWrapper, RDKitToolkitWrapper]) + molecule = create_ethanol() + molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry) + charge_sum = 0 * unit.elementary_charge + abs_charge_sum = 0 * unit.elementary_charge + for pc in molecule._partial_charges: + charge_sum += pc + abs_charge_sum += abs(pc) + assert abs(charge_sum) < 0.001 * unit.elementary_charge + assert abs_charge_sum > 0.25 * unit.elementary_charge - smiles = '[H]C([H])([H])C([H])([H])[H]' - molecule = Molecule.from_smiles(smiles, toolkit_registry=toolkit_registry) - molecule.generate_conformers(toolkit_registry=toolkit_registry) - - # TODO: Implementation of these tests is pending a decision on the API for our charge model - with pytest.raises(NotImplementedError) as excinfo: - charge_model = 'notARealChargeModel' - molecule.compute_partial_charges(toolkit_registry=toolkit_registry)#, charge_model=charge_model) - - # ['cm1', 'cm2'] - for charge_model in ['gas', 'mul', 'bcc']: - with pytest.raises(NotImplementedError) as excinfo: - molecule.compute_partial_charges(toolkit_registry=toolkit_registry)#, charge_model=charge_model) - charge_sum = 0 * unit.elementary_charge - for pc in molecule._partial_charges: - charge_sum += pc - assert charge_sum < 0.01 * unit.elementary_charge - - # For now, just test AM1-BCC while the SMIRNOFF spec for other charge models gets worked out - molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry) # , charge_model=charge_model) + @pytest.mark.skipif(not RDKitToolkitWrapper.is_available() or not AmberToolsToolkitWrapper.is_available(), + reason='RDKitToolkit and AmberToolsToolkit not available') + def test_compute_partial_charges_am1bcc_net_charge(self): + """Test AmberToolsToolkitWrapper assign_partial_charges() on a molecule with a net -1 charge""" + toolkit_registry = ToolkitRegistry(toolkit_precedence=[AmberToolsToolkitWrapper, RDKitToolkitWrapper]) + molecule = create_acetate() + molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry) charge_sum = 0 * unit.elementary_charge for pc in molecule._partial_charges: charge_sum += pc - assert charge_sum < 0.002 * unit.elementary_charge + assert -0.99 * unit.elementary_charge > charge_sum > -1.01 * unit.elementary_charge + @pytest.mark.skipif(not RDKitToolkitWrapper.is_available() or not AmberToolsToolkitWrapper.is_available(), + reason='RDKitToolkit and AmberToolsToolkit not available') + def test_compute_partial_charges_am1bcc_wrong_n_confs(self): + """ + Test AmberToolsToolkitWrapper compute_partial_charges_am1bcc() when requesting to use an incorrect number of + conformers + """ + from openforcefield.tests.test_forcefield import create_ethanol + toolkit_registry = ToolkitRegistry(toolkit_precedence=[AmberToolsToolkitWrapper, RDKitToolkitWrapper]) + molecule = create_ethanol() + molecule.generate_conformers(n_conformers=2, rms_cutoff=0.1*unit.angstrom) + + # Try passing in the incorrect number of confs, but without specifying strict_n_conformers, + # which should produce a warning + with pytest.warns(IncorrectNumConformersWarning, + match="has 2 conformers, but charge method 'am1bcc' expects exactly 1."): + molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry, + use_conformers=molecule.conformers, + strict_n_conformers=False) + + # Try again, with strict_n_confs as true, but not including use_confs, so the + # recommended number of confs will be generated + molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry, + strict_n_conformers=True) + + # Test calling the ToolkitWrapper _indirectly_, though the Molecule API, + # which should raise the first error encountered + with pytest.raises(ValueError, + match=f"has 2 conformers, but charge method 'am1bcc' " + f"expects exactly 1."): + molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry, + use_conformers=molecule.conformers, + strict_n_conformers=True) + + # Test calling the ToolkitWrapper _indirectly_, though a ToolkitRegistry, + # specifying raise_exception_types=[] + # which should aggregate any exceptions and bundle all of the messages + # in a failed task together in a single ValueError. + with pytest.raises(ValueError, + match=f"has 2 conformers, but charge method 'am1bcc' " + f"expects exactly 1."): + toolkit_registry.call('compute_partial_charges_am1bcc', + molecule=molecule, + use_conformers=molecule.conformers, + strict_n_conformers=True, + raise_exception_types=[]) + + # Test calling the ToolkitWrapper _directly_, passing in the incorrect number of + # confs, and specify strict_n_conformers, which should produce an IncorrectNumConformersError + with pytest.raises(IncorrectNumConformersError, + match=f"has 2 conformers, but charge method 'am1bcc' " + f"expects exactly 1."): + ATTKW = AmberToolsToolkitWrapper() + ATTKW.compute_partial_charges_am1bcc(molecule=molecule, + use_conformers=molecule.conformers, + strict_n_conformers=True) @pytest.mark.skipif(not RDKitToolkitWrapper.is_available() or not AmberToolsToolkitWrapper.is_available(), reason='RDKitToolkit and AmberToolsToolkit not available') - def test_compute_partial_charges_net_charge(self): - """Test OpenEyeToolkitWrapper compute_partial_charges() on a molecule with a net +1 charge""" + @pytest.mark.parametrize("partial_charge_method", ['am1bcc', 'am1-mulliken', 'gasteiger']) + def test_assign_partial_charges_neutral(self, partial_charge_method): + """Test AmberToolsToolkitWrapper assign_partial_charges()""" + from openforcefield.tests.test_forcefield import create_ethanol toolkit_registry = ToolkitRegistry(toolkit_precedence=[AmberToolsToolkitWrapper, RDKitToolkitWrapper]) - smiles = '[H]C([H])([H])[N+]([H])([H])[H]' - molecule = Molecule.from_smiles(smiles, toolkit_registry=toolkit_registry) - molecule.generate_conformers(toolkit_registry=toolkit_registry) + molecule = create_ethanol() + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method) + charge_sum = 0. * unit.elementary_charge + for pc in molecule.partial_charges: + charge_sum += pc + assert -1.e-5 < charge_sum.value_in_unit(unit.elementary_charge) < 1.e-5 + + @pytest.mark.skipif(not RDKitToolkitWrapper.is_available() or not AmberToolsToolkitWrapper.is_available(), + reason='RDKitToolkit and AmberToolsToolkit not available') + @pytest.mark.parametrize("partial_charge_method", ['am1bcc', 'am1-mulliken']) + def test_assign_partial_charges_conformer_dependence(self, partial_charge_method): + """Test AmberToolsToolkitWrapper assign_partial_charges()'s use_conformers kwarg + to ensure charges are really conformer dependent. Skip Gasteiger because it isn't + conformer dependent.""" + from openforcefield.tests.test_forcefield import create_ethanol + import copy + toolkit_registry = ToolkitRegistry(toolkit_precedence=[AmberToolsToolkitWrapper, RDKitToolkitWrapper]) + molecule = create_ethanol() + molecule.generate_conformers(n_conformers=1) + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method, + use_conformers=molecule.conformers) + pcs1 = copy.deepcopy(molecule.partial_charges) + # This test case needs a pretty extreme coordinate change since ambertools only + # stores partial charges to 1e-3 + molecule._conformers[0][0][0] += 3. * unit.angstrom + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method, + use_conformers=molecule.conformers) + for pc1, pc2 in zip(pcs1, molecule.partial_charges): + assert abs(pc1 - pc2) > 1.e-3 * unit.elementary_charge + + @pytest.mark.skipif(not RDKitToolkitWrapper.is_available() or not AmberToolsToolkitWrapper.is_available(), + reason='RDKitToolkit and AmberToolsToolkit not available') + @pytest.mark.parametrize("partial_charge_method", ['am1bcc', 'am1-mulliken', 'gasteiger']) + def test_assign_partial_charges_net_charge(self, partial_charge_method): + """ + Test AmberToolsToolkitWrapper assign_partial_charges(). + """ + from openforcefield.tests.test_forcefield import create_acetate + toolkit_registry = ToolkitRegistry(toolkit_precedence=[AmberToolsToolkitWrapper, RDKitToolkitWrapper]) + molecule = create_acetate() + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method) + charge_sum = 0. * unit.elementary_charge + for pc in molecule.partial_charges: + charge_sum += pc + assert -1.01 < charge_sum.value_in_unit(unit.elementary_charge) < -0.99 - with pytest.raises(NotImplementedError) as excinfo: - charge_model = 'notARealChargeModel' - molecule.compute_partial_charges(toolkit_registry=toolkit_registry)#, charge_model=charge_model) + @pytest.mark.skipif(not RDKitToolkitWrapper.is_available() or not AmberToolsToolkitWrapper.is_available(), + reason='RDKitToolkit and AmberToolsToolkit not available') + def test_assign_partial_charges_bad_charge_method(self): + """Test AmberToolsToolkitWrapper assign_partial_charges() for a nonexistent charge method""" + from openforcefield.tests.test_forcefield import create_ethanol + toolkit_registry = ToolkitRegistry(toolkit_precedence=[AmberToolsToolkitWrapper, RDKitToolkitWrapper]) + molecule = create_ethanol() - # TODO: Figure out why ['cm1', 'cm2'] fail - for charge_model in ['gas', 'mul', 'bcc']: - with pytest.raises(NotImplementedError) as excinfo: - molecule.compute_partial_charges(toolkit_registry=toolkit_registry)#, charge_model=charge_model) - charge_sum = 0 * unit.elementary_charge - for pc in molecule._partial_charges: - charge_sum += pc - assert 0.99 * unit.elementary_charge < charge_sum < 1.01 * unit.elementary_charge + # For now, ToolkitRegistries lose track of what exception type + # was thrown inside them, so we just check for a ValueError here + with pytest.raises(ValueError, match="is not available from AmberToolsToolkitWrapper") as excinfo: + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method="NotARealChargeMethod") - # For now, I'm just testing AM1-BCC (will test more when the SMIRNOFF spec for other charges is finalized) - molecule.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry) - charge_sum = 0 * unit.elementary_charge - for pc in molecule._partial_charges: - charge_sum += pc - assert 0.999 * unit.elementary_charge < charge_sum < 1.001 * unit.elementary_charge + # ToolkitWrappers raise a specific exception class, so we test that here + with pytest.raises(ChargeMethodUnavailableError, match="is not available from AmberToolsToolkitWrapper") as excinfo: + ATTKW = AmberToolsToolkitWrapper() + ATTKW.assign_partial_charges(molecule=molecule, + partial_charge_method="NotARealChargeMethod") + + + @pytest.mark.skipif(not RDKitToolkitWrapper.is_available() or not AmberToolsToolkitWrapper.is_available(), + reason='RDKitToolkit and AmberToolsToolkit not available') + @pytest.mark.parametrize("partial_charge_method,expected_n_confs", [('am1bcc', 1), + ('am1-mulliken', 1), + ('gasteiger', 0)]) + def test_assign_partial_charges_wrong_n_confs(self, partial_charge_method, expected_n_confs): + """ + Test AmberToolsToolkitWrapper assign_partial_charges() when requesting to use an incorrect number of + conformers + """ + from openforcefield.tests.test_forcefield import create_ethanol + toolkit_registry = ToolkitRegistry(toolkit_precedence=[AmberToolsToolkitWrapper, RDKitToolkitWrapper]) + molecule = create_ethanol() + molecule.generate_conformers(n_conformers=2, rms_cutoff=0.01*unit.angstrom) + + # Try passing in the incorrect number of confs, but without specifying strict_n_conformers, + # which should produce a warning + with pytest.warns(IncorrectNumConformersWarning, + match=f"has 2 conformers, but charge method '{partial_charge_method}' " + f"expects exactly {expected_n_confs}."): + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method, + use_conformers=molecule.conformers, + strict_n_conformers=False) + + # Try again, with strict_n_confs as true, but not including use_confs, so the + # recommended number of confs will be generated + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method, + strict_n_conformers=True) + + # Test calling the ToolkitWrapper _indirectly_, though the Molecule API + # which should aggregate any exceptions and bundle all of the messages + # in a failed task together in a single ValueError. + with pytest.raises(ValueError, + match=f"has 2 conformers, but charge method '{partial_charge_method}' " + f"expects exactly {expected_n_confs}."): + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method, + use_conformers=molecule.conformers, + strict_n_conformers=True) + + # Test calling the ToolkitWrapper _directly_, passing in the incorrect number of + # confs, and specify strict_n_conformers, which should produce an IncorrectNumConformersError + with pytest.raises(IncorrectNumConformersError, + match=f"has 2 conformers, but charge method '{partial_charge_method}' " + f"expects exactly {expected_n_confs}."): + ATTKW = AmberToolsToolkitWrapper() + ATTKW.assign_partial_charges(molecule=molecule, + partial_charge_method=partial_charge_method, + use_conformers=molecule.conformers, + strict_n_conformers=True) @pytest.mark.skipif(not RDKitToolkitWrapper.is_available() or not AmberToolsToolkitWrapper.is_available(), reason='RDKitToolkit and AmberToolsToolkit not available') @@ -1765,8 +2072,188 @@ def test_assign_fractional_bond_orders_double_bond(self): assert double_bond_has_wbo_near_2 +class TestBuiltInToolkitWrapper: + """Test the BuiltInToolkitWrapper""" + @pytest.mark.parametrize("partial_charge_method", ['zeros', 'formal_charge']) + def test_assign_partial_charges_neutral(self, partial_charge_method): + """Test BuiltInToolkitWrapper assign_partial_charges()""" + from openforcefield.tests.test_forcefield import create_ethanol + toolkit_registry = ToolkitRegistry(toolkit_precedence=[BuiltInToolkitWrapper]) + molecule = create_ethanol() + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method) + charge_sum = 0. * unit.elementary_charge + for pc in molecule.partial_charges: + charge_sum += pc + assert -1.e-6 < charge_sum.value_in_unit(unit.elementary_charge) < 1.e-6 + + @pytest.mark.parametrize("partial_charge_method", ['formal_charge']) + def test_assign_partial_charges_net_charge(self, partial_charge_method): + """ + Test BuiltInToolkitWrapper assign_partial_charges(). Only formal_charge is tested, since zeros will not + sum up to the proper number + """ + from openforcefield.tests.test_forcefield import create_acetate + toolkit_registry = ToolkitRegistry(toolkit_precedence=[BuiltInToolkitWrapper]) + molecule = create_acetate() + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method=partial_charge_method) + charge_sum = 0. * unit.elementary_charge + for pc in molecule.partial_charges: + charge_sum += pc + assert -1.e-6 < charge_sum.value_in_unit(unit.elementary_charge) + 1. < 1.e-6 + + + def test_assign_partial_charges_bad_charge_method(self): + """Test BuiltInToolkitWrapper assign_partial_charges() for a nonexistent charge method""" + from openforcefield.tests.test_forcefield import create_ethanol + toolkit_registry = ToolkitRegistry(toolkit_precedence=[BuiltInToolkitWrapper]) + molecule = create_ethanol() + + # For now, the Molecule API passes raise_exception_types=[] to ToolkitRegistry.call, + # which loses track of what exception type + # was thrown inside them, so we just check for a ValueError here + with pytest.raises(ValueError, match="is not supported by the Built-in toolkit") as excinfo: + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method="NotARealChargeMethod") + + # ToolkitWrappers raise a specific exception class, so we test that here + with pytest.raises(ChargeMethodUnavailableError, match="is not supported by the Built-in toolkit") as excinfo: + BITKW = BuiltInToolkitWrapper() + BITKW.assign_partial_charges(molecule=molecule, + partial_charge_method="NotARealChargeMethod") + + def test_assign_partial_charges_wrong_n_confs(self): + """ + Test BuiltInToolkitWrapper assign_partial_charges() when requesting to use an incorrect number of + conformers + """ + from openforcefield.tests.test_forcefield import create_ethanol + toolkit_registry = ToolkitRegistry(toolkit_precedence=[BuiltInToolkitWrapper]) + molecule = create_ethanol() + molecule.generate_conformers(n_conformers=1) + with pytest.warns(IncorrectNumConformersWarning, + match="has 1 conformers, but charge method 'zeros' expects exactly 0."): + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method="zeros", + use_conformers=molecule.conformers, + strict_n_conformers=False) + + # Specify strict_n_conformers=True, but not use_conformers, so a recommended number of + # conformers will be generated internally + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method="zeros", + strict_n_conformers=True) + + # For now, the Molecule API passes raise_exception_types=[] to ToolkitRegistry.call, + # which loses track of what exception type + # was thrown inside them, so we just check for a ValueError here + with pytest.raises(ValueError, + match=f"has 1 conformers, but charge method 'zeros' " + f"expects exactly 0."): + molecule.assign_partial_charges(toolkit_registry=toolkit_registry, + partial_charge_method="zeros", + use_conformers=molecule.conformers, + strict_n_conformers=True) + + # Test calling the ToolkitWrapper _directly_, passing in the incorrect number of + # confs, and specify strict_n_conformers, which should produce an IncorrectNumConformersError + with pytest.raises(IncorrectNumConformersError, + match=f"has 1 conformers, but charge method 'zeros' " + f"expects exactly 0."): + BITKW = BuiltInToolkitWrapper() + BITKW.assign_partial_charges(molecule=molecule, + partial_charge_method="zeros", + use_conformers=molecule.conformers, + strict_n_conformers=True) + +class TestToolkitWrapper: + """Test the ToolkitWrapper class""" + def test_check_n_conformers(self): + """Ensure that _check_n_conformers is working properly""" + tkw = ToolkitWrapper() + mol = create_ethanol() + + ## Test molecule with no conformers + # Check with no min or max should pass + tkw._check_n_conformers(mol, 'nocharge') + # Check with min=1 should warn + with pytest.warns(IncorrectNumConformersWarning, + match="has 0 conformers, but charge method 'nocharge' expects at least 1"): + tkw._check_n_conformers(mol, 'nocharge', min_confs=1) + # Check with min=1 and strict_n_conformers should raise an error + with pytest.raises(IncorrectNumConformersError, + match="has 0 conformers, but charge method 'nocharge' expects at least 1"): + tkw._check_n_conformers(mol, 'nocharge', min_confs=1, strict_n_conformers=True) + # Check with min=1, max=1 and strict_n_conformers should raise an error + with pytest.raises(IncorrectNumConformersError, + match="has 0 conformers, but charge method 'nocharge' expects exactly 1"): + tkw._check_n_conformers(mol, 'nocharge', min_confs=1, max_confs=1, strict_n_conformers=True) + # Check with min=1, max=2 and strict_n_conformers should raise an error + with pytest.raises(IncorrectNumConformersError, + match="has 0 conformers, but charge method 'nocharge' expects between 1 and 2"): + tkw._check_n_conformers(mol, 'nocharge', min_confs=1, max_confs=2, strict_n_conformers=True) + # Check with max=1 should pass + tkw._check_n_conformers(mol, 'nocharge', max_confs=1, strict_n_conformers=True) + + ## Test molecule with conformers + # Add some conformers + mol.generate_conformers(n_conformers=1) + for _ in range(9): + mol.add_conformer(mol.conformers[0]) + + # Check with no min or max should pass + tkw._check_n_conformers(mol, 'nocharge') + + ## min_confs checks + # Check with min=1 should be fine + tkw._check_n_conformers(mol, 'nocharge', min_confs=1) + # Check with min=10 should be fine + tkw._check_n_conformers(mol, 'nocharge', min_confs=10) + # Check with min=11 should warn + with pytest.warns(IncorrectNumConformersWarning, + match="has 10 conformers, but charge method 'nocharge' expects at least 11"): + tkw._check_n_conformers(mol, 'nocharge', min_confs=11) + # Check with min=11 and strict_n_conformers should raise an error + with pytest.raises(IncorrectNumConformersError, + match="has 10 conformers, but charge method 'nocharge' expects at least 11"): + tkw._check_n_conformers(mol, 'nocharge', min_confs=11, strict_n_conformers=True) + + ## max_confs checks + # Check with max=1 and strict_n_conformers should raise an error + with pytest.raises(IncorrectNumConformersError, + match="has 10 conformers, but charge method 'nocharge' expects at most 1"): + tkw._check_n_conformers(mol, 'nocharge', max_confs=1, strict_n_conformers=True) + # Check with max=10 and strict_n_conformers should be OK + tkw._check_n_conformers(mol, 'nocharge', max_confs=10, strict_n_conformers=True) + # Check with max=11 and strict_n_conformers should be OK + tkw._check_n_conformers(mol, 'nocharge', max_confs=11, strict_n_conformers=True) + + ## min_confs and max_confs checks + # Check with max=10 and min=10 and strict_n_conformers should be OK + tkw._check_n_conformers(mol, 'nocharge', min_confs=10, max_confs=10, strict_n_conformers=True) + # Check with max=10 and min=9 and strict_n_conformers should be OK + tkw._check_n_conformers(mol, 'nocharge', min_confs=9, max_confs=10, strict_n_conformers=True) + # Check with max=11 and min=10 and strict_n_conformers should be OK + tkw._check_n_conformers(mol, 'nocharge', min_confs=10, max_confs=11, strict_n_conformers=True) + # Check with max=11 and min=9 and strict_n_conformers should be OK + tkw._check_n_conformers(mol, 'nocharge', min_confs=9, max_confs=11, strict_n_conformers=True) + # Check with min=9 and max=9 and strict_n_conformers should raise an error + with pytest.raises(IncorrectNumConformersError, + match="has 10 conformers, but charge method 'nocharge' expects exactly 9"): + tkw._check_n_conformers(mol, 'nocharge', min_confs=9, max_confs=9, strict_n_conformers=True) + # Check with min=1 and max=9 and strict_n_conformers should raise an error + with pytest.raises(IncorrectNumConformersError, + match="has 10 conformers, but charge method 'nocharge' expects between 1 and 9"): + tkw._check_n_conformers(mol, 'nocharge', min_confs=1, max_confs=9, strict_n_conformers=True) + # Check with min=11 and max=12 and strict_n_conformers should raise an error + with pytest.raises(IncorrectNumConformersError, + match="has 10 conformers, but charge method 'nocharge' expects between 11 and 12"): + tkw._check_n_conformers(mol, 'nocharge', min_confs=11, max_confs=12, strict_n_conformers=True) + + class TestToolkitRegistry: - """Test the ToolkitRegistry""" + """Test the ToolkitRegistry class""" @pytest.mark.skipif(not OpenEyeToolkitWrapper.is_available(), reason='OpenEye Toolkit not available') def test_register_openeye(self): @@ -1819,8 +2306,8 @@ def test_register_ambertools(self): assert set([ type(c) for c in registry.registered_toolkits]) == set([AmberToolsToolkitWrapper,RDKitToolkitWrapper]) # Test ToolkitRegistry.resolve() - registry.resolve('compute_partial_charges') - assert registry.resolve('compute_partial_charges') == registry.registered_toolkits[0].compute_partial_charges + registry.resolve('assign_partial_charges') + assert registry.resolve('assign_partial_charges') == registry.registered_toolkits[0].assign_partial_charges # Test ToolkitRegistry.call() registry.register_toolkit(RDKitToolkitWrapper) @@ -1828,6 +2315,39 @@ def test_register_ambertools(self): molecule = registry.call('from_smiles', smiles) #partial_charges = registry.call('compute_partial_charges', molecule) + def test_register_builtintoolkit(self): + """Test creation of toolkit registry with Built-in toolkit""" + # Test registration of BuiltInToolkitWrapper + toolkit_precedence = [BuiltInToolkitWrapper] + registry = ToolkitRegistry(toolkit_precedence=toolkit_precedence, + register_imported_toolkit_wrappers=False) + #registry.register_toolkit(BuiltInToolkitWrapper) + assert set([ type(c) for c in registry.registered_toolkits]) == set([BuiltInToolkitWrapper]) + + # Test ToolkitRegistry.resolve() + assert registry.resolve('assign_partial_charges') == registry.registered_toolkits[0].assign_partial_charges + + @pytest.mark.skipif(not AmberToolsToolkitWrapper.is_available(), reason='AmberTools Toolkit not available') + def test_call_raise_first_error(self): + """Test to ensure proper behavior of raise_first_error kwarg to ToolkitRegistry.call""" + toolkit_precedence = [BuiltInToolkitWrapper, RDKitToolkitWrapper, AmberToolsToolkitWrapper] + registry = ToolkitRegistry(toolkit_precedence=toolkit_precedence, + register_imported_toolkit_wrappers=False) + mol = registry.call('from_smiles', 'C') + # Specify that the ToolkitRegistry should raise the first ChargeMethodUnavailableError it encounters + with pytest.raises(ChargeMethodUnavailableError, match='"notarealchargemethod"" is not supported by the Built-in toolkit.'): + registry.call('assign_partial_charges', + molecule=mol, + partial_charge_method="NotARealChargeMethod", + raise_exception_types=[ChargeMethodUnavailableError]) + # Specify that the ToolkitRegistry should collect all the errors it encounters and + # ensure it raises a single ValueError when no ToolkitWrappers succeed + with pytest.raises(ValueError, match="partial_charge_method \'notarealchargemethod\' is not available from AmberToolsToolkitWrapper"): + registry.call('assign_partial_charges', + molecule=mol, + partial_charge_method="NotARealChargeMethod", + raise_exception_types=[]) + @pytest.mark.skipif(not RDKitToolkitWrapper.is_available(), reason='RDKit Toolkit not available') def test_substructure_search_on_large_molecule(self): """Test RDKitToolkitWrapper substructure search when a large number hits are found""" @@ -1837,7 +2357,6 @@ def test_substructure_search_on_large_molecule(self): molecule = tk.from_smiles(smiles) query = "[C:1]~[C:2]" ret = molecule.chemical_environment_matches(query, toolkit_registry=tk) - assert len(ret) == 5998 + assert len(ret) == 5998 assert len(ret[0]) == 2 - diff --git a/openforcefield/tests/test_utils.py b/openforcefield/tests/test_utils.py index fce047041..24b5d9112 100644 --- a/openforcefield/tests/test_utils.py +++ b/openforcefield/tests/test_utils.py @@ -40,6 +40,14 @@ class FooSubSubclass(FooSubclass1): subclass_names = [ cls.__name__ for cls in utils.all_subclasses(Foo) ] assert set(subclass_names) == set(['FooSubclass1', 'FooSubclass2', 'FooSubSubclass']) +def test_temporary_cd(): + """Test temporary_cd() context manager""" + initial_dir = os.getcwd() + temporary_dir = '/' + from openforcefield.utils import temporary_cd + with temporary_cd(temporary_dir): + assert os.getcwd() == temporary_dir + assert os.getcwd() == initial_dir def test_get_data_file_path(): """Test get_data_file_path()""" diff --git a/openforcefield/topology/__init__.py b/openforcefield/topology/__init__.py index 8e9856250..2bf4fba40 100644 --- a/openforcefield/topology/__init__.py +++ b/openforcefield/topology/__init__.py @@ -6,6 +6,6 @@ from openforcefield.topology.topology import ( DuplicateUniqueMoleculeError, NotBondedError, - ValenceDict, ImproperDict, + ValenceDict, ImproperDict, SortedDict, TopologyAtom, TopologyBond, TopologyVirtualSite, TopologyMolecule, Topology ) diff --git a/openforcefield/topology/molecule.py b/openforcefield/topology/molecule.py index 26bada5bf..d7065cbe5 100644 --- a/openforcefield/topology/molecule.py +++ b/openforcefield/topology/molecule.py @@ -48,13 +48,15 @@ from networkx.algorithms.isomorphism import GraphMatcher import openforcefield -from openforcefield.utils import serialize_numpy, deserialize_numpy, quantity_to_string, string_to_quantity -from openforcefield.utils.toolkits import ToolkitRegistry, ToolkitWrapper, RDKitToolkitWrapper, OpenEyeToolkitWrapper, \ - InvalidToolkitError, GLOBAL_TOOLKIT_REGISTRY +from openforcefield.utils import serialize_numpy, deserialize_numpy, quantity_to_string, \ + string_to_quantity, check_units_are_compatible +from openforcefield.utils.toolkits import ToolkitRegistry, ToolkitWrapper, RDKitToolkitWrapper, OpenEyeToolkitWrapper,\ + InvalidToolkitError, UndefinedStereochemistryError, GLOBAL_TOOLKIT_REGISTRY from openforcefield.utils.toolkits import DEFAULT_AROMATICITY_MODEL from openforcefield.utils.serialization import Serializable + #============================================================================================= # GLOBAL PARAMETERS #============================================================================================= @@ -175,7 +177,7 @@ def __init__(self, ---------- atomic_number : int Atomic number of the atom - formal_charge : int + formal_charge : int or simtk.unit.Quantity-wrapped int with dimension "charge" Formal charge of the atom is_aromatic : bool If True, atom is aromatic; if False, not aromatic @@ -197,7 +199,8 @@ def __init__(self, """ self._atomic_number = atomic_number - self._formal_charge = formal_charge + # Use the setter here, since it will handle either ints or Quantities + self.formal_charge = formal_charge self._is_aromatic = is_aromatic self._stereochemistry = stereochemistry if name is None: @@ -245,7 +248,7 @@ def to_dict(self): # TODO atom_dict = OrderedDict() atom_dict['atomic_number'] = self._atomic_number - atom_dict['formal_charge'] = self._formal_charge + atom_dict['formal_charge'] = self._formal_charge / unit.elementary_charge atom_dict['is_aromatic'] = self._is_aromatic atom_dict['stereochemistry'] = self._stereochemistry # TODO: Should we let atoms have names? @@ -270,6 +273,17 @@ def formal_charge(self): """ return self._formal_charge + @formal_charge.setter + def formal_charge(self, other): + """ + Set the atom's formal charge. Accepts either ints or simtk.unit.Quantity-wrapped ints with units of charge. + """ + if isinstance(other, int): + self._formal_charge = other * unit.elementary_charge + else: + check_units_are_compatible("formal charge", other, unit.elementary_charge) + self._formal_charge = other + @property def partial_charge(self): """ @@ -1506,6 +1520,10 @@ def __init__(self, self._initialize() else: loaded = False + # Start a list of the ValueErrors the following logic encounters, so we can print it out + # if there turned out to be no way to load this input + value_errors = list() + if isinstance( other, openforcefield.topology.FrozenMolecule) and not (loaded): @@ -1519,36 +1537,60 @@ def __init__(self, if isinstance(other, OrderedDict) and not (loaded): self.__setstate__(other) loaded = True + + # Check through the toolkit registry to find a compatible wrapper for loading if not loaded: try: - result = toolkit_registry.call('from_object', other, allow_undefined_stereo=allow_undefined_stereo) - except NotImplementedError: - pass + # Each ToolkitWrapper may provide a from_object method, which turns some particular type(s) + # of object into OFFMols. For example, RDKitToolkitWrapper's from_object method will + # return an OFFMol if provided with an RDMol, or raise a ValueError if it is provided + # an OEMol (or anything else). This makes the assumption that any non-ValueError errors raised + # by the toolkit _really are_ bad and should be raised immediately, which may be a bad assumption. + result = toolkit_registry.call('from_object', + other, + allow_undefined_stereo=allow_undefined_stereo, + raise_exception_types=[UndefinedStereochemistryError]) + # NotImplementedError should never be raised... Only from_file and from_file_obj are provided + # in the base ToolkitWrapper class and require overwriting, so from_object should be excluded + # except NotImplementedError as e: + # raise e + # The toolkit registry will aggregate all errors except UndefinedStereochemistryErrors into a single + # ValueError, which we should catch and and store that here. + except ValueError as e: + value_errors.append(e) else: self._copy_initializer(result) loaded = True # TODO: Make this compatible with file-like objects (I couldn't figure out how to make an oemolistream # from a fileIO object) - if (isinstance(other, str) - or hasattr(other, 'read')) and not (loaded): - mol = Molecule.from_file( - other, - file_format=file_format, - toolkit_registry=toolkit_registry, - allow_undefined_stereo=allow_undefined_stereo - ) # returns a list only if multiple molecules are found - if type(mol) == list: - raise ValueError( - 'Specified file or file-like object must contain exactly one molecule' - ) - - self._copy_initializer(mol) - loaded = True + if (isinstance(other, str) or hasattr(other, 'read')) and not (loaded): + try: + mol = Molecule.from_file( + other, + file_format=file_format, + toolkit_registry=toolkit_registry, + allow_undefined_stereo=allow_undefined_stereo + ) # returns a list only if multiple molecules are found + if type(mol) == list: + raise ValueError( + 'Specified file or file-like object must contain exactly one molecule' + ) + except ValueError as e: + value_errors.append(e) + else: + self._copy_initializer(mol) + loaded = True + + # If none of the above methods worked, raise a ValueError summarizing the + # errors from the different loading attempts + if not (loaded): msg = 'Cannot construct openforcefield.topology.Molecule from {}\n'.format( other) - raise Exception(msg) + for value_error in value_errors: + msg += str(value_error) + raise ValueError(msg) @property def has_unique_atom_names(self): @@ -2253,9 +2295,11 @@ def is_isomorphic_with(self, other, **kwargs): def generate_conformers(self, toolkit_registry=GLOBAL_TOOLKIT_REGISTRY, n_conformers=10, + rms_cutoff=None, clear_existing=True): """ - Generate conformers for this molecule using an underlying toolkit + Generate conformers for this molecule using an underlying toolkit. If n_conformers=0, no toolkit wrapper + will be called. If n_conformers=0 and clear_existing=True, molecule.conformers will be set to None. Parameters ---------- @@ -2263,6 +2307,10 @@ def generate_conformers(self, :class:`ToolkitRegistry` or :class:`ToolkitWrapper` to use for SMILES-to-molecule conversion n_conformers : int, default=1 The maximum number of conformers to produce + rms_cutoff : simtk.Quantity-wrapped float, in units of distance, optional, default=None + The minimum RMS value at which two conformers are considered redundant and one is deleted. Precise + implementation of this cutoff may be toolkit-dependent. If None, the cutoff is set to be the default value + for each ToolkitWrapper (generally 1 Angstrom). clear_existing : bool, default=True Whether to overwrite existing conformers for the molecule @@ -2278,23 +2326,45 @@ def generate_conformers(self, If an invalid object is passed as the toolkit_registry parameter """ + # If no conformers are requested, do not call to a ToolkitWrapper at all + if n_conformers == 0: + if clear_existing: + self._conformers = None + return + if isinstance(toolkit_registry, ToolkitRegistry): - return toolkit_registry.call('generate_conformers', self, n_conformers=n_conformers, + return toolkit_registry.call('generate_conformers', + self, + n_conformers=n_conformers, + rms_cutoff=rms_cutoff, clear_existing=clear_existing) elif isinstance(toolkit_registry, ToolkitWrapper): toolkit = toolkit_registry - return toolkit.generate_conformers(self, n_conformers=n_conformers, clear_existing=clear_existing) + return toolkit.generate_conformers(self, + n_conformers=n_conformers, + rms_cutoff = rms_cutoff, + clear_existing=clear_existing) else: raise InvalidToolkitError( 'Invalid toolkit_registry passed to generate_conformers. Expected ToolkitRegistry or ToolkitWrapper. Got {}' .format(type(toolkit_registry))) - def compute_partial_charges_am1bcc(self, toolkit_registry=GLOBAL_TOOLKIT_REGISTRY): + def compute_partial_charges_am1bcc(self, + use_conformers=None, + strict_n_conformers=False, + toolkit_registry=GLOBAL_TOOLKIT_REGISTRY): """ Calculate partial atomic charges for this molecule using AM1-BCC run by an underlying toolkit + and assign them to this molecule's partial_charges attribute. Parameters ---------- + strict_n_conformers : bool, default=False + Whether to raise an exception if an invalid number of conformers is provided for the given charge method. + If this is False and an invalid number of conformers is found, a warning will be raised. + use_conformers : iterable of simtk.unit.Quantity-wrapped numpy arrays, each with shape (n_atoms, 3) and dimension of distance. Optional, default=None + Coordinates to use for partial charge calculation. + If None, an appropriate number of conformers for the given charge method will be generated. toolkit_registry : openforcefield.utils.toolkits.ToolkitRegistry or openforcefield.utils.toolkits.ToolkitWrapper, optional, default=None :class:`ToolkitRegistry` or :class:`ToolkitWrapper` to use for the calculation @@ -2311,42 +2381,38 @@ def compute_partial_charges_am1bcc(self, toolkit_registry=GLOBAL_TOOLKIT_REGISTR If an invalid object is passed as the toolkit_registry parameter """ - if isinstance(toolkit_registry, ToolkitRegistry): - charges = toolkit_registry.call( - 'compute_partial_charges_am1bcc', - self - ) - elif isinstance(toolkit_registry, ToolkitWrapper): - toolkit = toolkit_registry - charges = toolkit.compute_partial_charges_am1bcc(self) - else: - raise InvalidToolkitError( - 'Invalid toolkit_registry passed to compute_partial_charges_am1bcc. Expected ToolkitRegistry or ToolkitWrapper. Got {}' - .format(type(toolkit_registry))) - self.partial_charges = charges + self.assign_partial_charges(partial_charge_method='am1bcc', + use_conformers=use_conformers, + strict_n_conformers=strict_n_conformers, + toolkit_registry=toolkit_registry) - def compute_partial_charges(self, - #quantum_chemical_method='AM1-BCC', - #partial_charge_method='None', - toolkit_registry=GLOBAL_TOOLKIT_REGISTRY): + + def assign_partial_charges(self, + partial_charge_method, + strict_n_conformers=False, + use_conformers=None, + toolkit_registry=GLOBAL_TOOLKIT_REGISTRY): """ - **Warning! Not Implemented!** - Calculate partial atomic charges for this molecule using an underlying toolkit + Calculate partial atomic charges for this molecule using an underlying toolkit, and assign + the new values to the partial_charges attribute. Parameters ---------- - quantum_chemical_method : string, default='AM1-BCC' - The quantum chemical method to use for partial charge calculation. - partial_charge_method : string, default='None' + partial_charge_method : string The partial charge calculation method to use for partial charge calculation. + strict_n_conformers : bool, default=False + Whether to raise an exception if an invalid number of conformers is provided for the given charge method. + If this is False and an invalid number of conformers is found, a warning will be raised. + use_conformers : iterable of simtk.unit.Quantity-wrapped numpy arrays, each with shape (n_atoms, 3) and dimension of distance. Optional, default=None + Coordinates to use for partial charge calculation. If None, an appropriate number of conformers will be generated. toolkit_registry : openforcefield.utils.toolkits.ToolkitRegistry or openforcefield.utils.toolkits.ToolkitWrapper, optional, default=None - :class:`ToolkitRegistry` or :class:`ToolkitWrapper` to use for SMILES-to-molecule conversion + :class:`ToolkitRegistry` or :class:`ToolkitWrapper` to use for the calculation. Examples -------- >>> molecule = Molecule.from_smiles('CCCCCC') - >>> molecule.generate_conformers() + >>> molecule.assign_partial_charges('am1-mulliken') Raises ------ @@ -2354,24 +2420,30 @@ def compute_partial_charges(self, If an invalid object is passed as the toolkit_registry parameter """ - raise NotImplementedError - # TODO: Implement this in a way that's compliant with SMIRNOFF's tag when the spec gets finalized - # if isinstance(toolkit_registry, ToolkitRegistry): - # charges = toolkit_registry.call( - # 'compute_partial_charges_am1bcc', - # self, - # ) - # elif isinstance(toolkit_registry, ToolkitWrapper): - # toolkit = toolkit_registry - # charges = toolkit.compute_partial_charges_am1bcc( - # self, - # #quantum_chemical_method=quantum_chemical_method, - # #partial_charge_method=partial_charge_method - # ) - # else: - # raise InvalidToolkitError( - # 'Invalid toolkit_registry passed to compute_partial_charges_am1bcc. Expected ToolkitRegistry or ToolkitWrapper. Got {}' - # .format(type(toolkit_registry))) + if isinstance(toolkit_registry, ToolkitRegistry): + # We may need to try several toolkitwrappers to find one + # that supports the desired partial charge method, so we + # tell the ToolkitRegistry to continue trying ToolkitWrappers + # if one raises an error (raise_exception_types=[]) + toolkit_registry.call('assign_partial_charges', + self, + partial_charge_method=partial_charge_method, + use_conformers=use_conformers, + strict_n_conformers=strict_n_conformers, + raise_exception_types=[] + ) + elif isinstance(toolkit_registry, ToolkitWrapper): + toolkit = toolkit_registry + toolkit.assign_partial_charges(self, + partial_charge_method=partial_charge_method, + use_conformers=use_conformers, + strict_n_conformers=strict_n_conformers + ) + else: + raise InvalidToolkitError( + f'Invalid toolkit_registry passed to assign_partial_charges.' + f'Expected ToolkitRegistry or ToolkitWrapper. Got {type(toolkit_registry)}') + def assign_fractional_bond_orders(self, bond_order_model=None, @@ -3070,7 +3142,10 @@ def total_charge(self): """ Return the total charge on the molecule """ - return sum([atom.formal_charge for atom in self.atoms]) + charge_sum = 0. * unit.elementary_charge + for atom in self.atoms: + charge_sum += atom.formal_charge + return charge_sum @property def name(self): @@ -3851,8 +3926,8 @@ def to_qcschema(self, multiplicity=1, conformer=0, extras=None): raise InvalidConformerError('The molecule must have a conformation to produce a valid qcschema; ' f'no conformer was found at index {conformer}.') - # Gather the required qschema data - charge = sum([atom.formal_charge for atom in self.atoms]) + # Gather the required qcschema data + charge = self.total_charge / unit.elementary_charge connectivity = [(bond.atom1_index, bond.atom2_index, bond.bond_order) for bond in self.bonds] symbols = [Element.getByAtomicNumber(atom.atomic_number).symbol for atom in self.atoms] diff --git a/openforcefield/topology/topology.py b/openforcefield/topology/topology.py index 86c8e9acf..71c1063df 100644 --- a/openforcefield/topology/topology.py +++ b/openforcefield/topology/topology.py @@ -109,6 +109,15 @@ def __keytransform__(self, key): key = tuple(reversed(key)) return key +class SortedDict(_TransformedDict): + """Enforce uniqueness of atom index tuples, without any restrictions on atom reordering.""" + + def __keytransform__(self, key): + """Sort tuple from lowest to highest.""" + # Ensure key is a tuple. + key = tuple(sorted(key)) + # Reverse the key if the first element is bigger than the last. + return key class ImproperDict(_TransformedDict): """Symmetrize improper torsions.""" diff --git a/openforcefield/typing/engines/smirnoff/forcefield.py b/openforcefield/typing/engines/smirnoff/forcefield.py index b106fe902..6bd598b5a 100644 --- a/openforcefield/typing/engines/smirnoff/forcefield.py +++ b/openforcefield/typing/engines/smirnoff/forcefield.py @@ -885,8 +885,12 @@ def _load_smirnoff_data(self, smirnoff_data, allow_cosmetic_attributes=False): parameter_list_dict = {} else: parameter_list_dict = section_dict.pop(parameter_list_tagname, {}) - # Must be wrapped into its own tag. - parameter_list_dict = {parameter_list_tagname: parameter_list_dict} + + # If the parameter list isn't empty, it must be transferred into its own tag. + # This is necessary for deserializing SMIRNOFF force field sections which may or may + # not have any smirks-based elements (like an empty ChargeIncrementModel section) + if parameter_list_dict != {}: + parameter_list_dict = {parameter_list_tagname: parameter_list_dict} # Retrieve or create parameter handler, passing in section_dict to check for # compatibility if a handler for this parameter name already exists diff --git a/openforcefield/typing/engines/smirnoff/parameters.py b/openforcefield/typing/engines/smirnoff/parameters.py index 6a1170ab1..5284e4949 100644 --- a/openforcefield/typing/engines/smirnoff/parameters.py +++ b/openforcefield/typing/engines/smirnoff/parameters.py @@ -50,9 +50,9 @@ from simtk import openmm, unit from openforcefield.utils import attach_units, \ - extract_serialized_units_from_dict, ToolkitUnavailableException, MessageException, \ - object_to_quantity -from openforcefield.topology import ValenceDict, ImproperDict + extract_serialized_units_from_dict, MessageException, \ + object_to_quantity, GLOBAL_TOOLKIT_REGISTRY +from openforcefield.topology import ValenceDict, ImproperDict, SortedDict from openforcefield.topology.molecule import Molecule from openforcefield.typing.chemistry import ChemicalEnvironment from openforcefield.utils import IncompatibleUnitError @@ -3229,7 +3229,7 @@ def postprocess_system(self, system, topology, **kwargs): for top_particle in top_mol.particles: q, _, _ = force.getParticleParameters(top_particle.topology_particle_index) partial_charge_sum += q - if abs(float(formal_charge_sum) - (partial_charge_sum/unit.elementary_charge)) > 0.01: + if abs(formal_charge_sum - partial_charge_sum) > 0.01 * unit.elementary_charge: msg = f"Partial charge sum ({partial_charge_sum}) " \ f"for molecule '{top_mol.reference_molecule.name}' (SMILES " \ f"{top_mol.reference_molecule.to_smiles()} does not equal formal charge sum " \ @@ -3255,6 +3255,13 @@ class LibraryChargeType(ParameterType): name = ParameterAttribute(default=None) charge = IndexedParameterAttribute(unit=unit.elementary_charge) + def __init__(self, **kwargs): + super().__init__(**kwargs) + unique_tags, connectivity = GLOBAL_TOOLKIT_REGISTRY.call('get_tagged_smarts_connectivity', self.smirks) + if len(self.charge) != len(unique_tags): + raise SMIRNOFFSpecError(f"LibraryCharge {self} was initialized with unequal number of " + f"tagged atoms and charges") + _TAGNAME = 'LibraryCharges' # SMIRNOFF tag name to process _INFOTYPE = LibraryChargeType # info type to store @@ -3275,12 +3282,10 @@ def find_matches(self, entity): matching the tuple of particle indices in ``entity``. """ - # TODO: Right now, this method is only ever called with an entity that is a Topoogy. - # Should we reduce its scope and have a check here to make sure entity is a Topology? return self._find_matches(entity, transformed_dict_cls=dict) def create_force(self, system, topology, **kwargs): - from openforcefield.topology import FrozenMolecule, TopologyAtom, TopologyVirtualSite + from openforcefield.topology import FrozenMolecule force = super().create_force(system, topology, **kwargs) @@ -3369,7 +3374,7 @@ def check_handler_compatibility(self, pass def create_force(self, system, topology, **kwargs): - + import warnings from openforcefield.utils.toolkits import GLOBAL_TOOLKIT_REGISTRY from openforcefield.topology import FrozenMolecule, TopologyAtom, TopologyVirtualSite @@ -3387,8 +3392,13 @@ def create_force(self, system, topology, **kwargs): # If the molecule wasn't already assigned charge values, calculate them here toolkit_registry = kwargs.get('toolkit_registry', GLOBAL_TOOLKIT_REGISTRY) - temp_mol.generate_conformers(n_conformers=10, toolkit_registry=toolkit_registry) - temp_mol.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry) + try: + # We don't need to generate conformers here, since that will be done by default in + # compute_partial_charges_am1bcc if the use_conformers kwarg isn't defined + temp_mol.compute_partial_charges_am1bcc(toolkit_registry=toolkit_registry) + except Exception as e: + warnings.warn(str(e), Warning) + continue # Assign charges to relevant atoms for topology_molecule in topology._reference_molecule_to_topology_molecules[ref_mol]: @@ -3458,30 +3468,29 @@ class ChargeIncrementType(ParameterType): .. warning :: This API is experimental and subject to change. """ - _VALENCE_TYPE = 'Bond' # ChemicalEnvironment valence type expected for SMARTS + _VALENCE_TYPE = None # This disables the connectivity check when parsing LibraryChargeType objects _ELEMENT_NAME = 'ChargeIncrement' - chargeincrement = IndexedParameterAttribute(unit=unit.elementary_charge) + charge_increment = IndexedParameterAttribute(unit=unit.elementary_charge) + + def __init__(self, **kwargs): + super().__init__(**kwargs) + unique_tags, connectivity = GLOBAL_TOOLKIT_REGISTRY.call('get_tagged_smarts_connectivity', self.smirks) + if len(self.charge_increment) != len(unique_tags): + raise SMIRNOFFSpecError(f"ChargeIncrement {self} was initialized with unequal number of " + f"tagged atoms and charge increments") _TAGNAME = 'ChargeIncrementModel' # SMIRNOFF tag name to process _INFOTYPE = ChargeIncrementType # info type to store - # TODO: The structure of this is still undecided + _DEPENDENCIES = [vdWHandler, ElectrostaticsHandler, LibraryChargeHandler, ToolkitAM1BCCHandler] + + number_of_conformers = ParameterAttribute(default=1, converter=int) - number_of_conformers = ParameterAttribute(default=10, converter=int) - quantum_chemical_method = ParameterAttribute( - default='AM1', - converter=_allow_only(['AM1']) - ) partial_charge_method = ParameterAttribute( - default='CM2', - converter=_allow_only(['CM2']) + default='AM1-Mulliken', converter=str ) - def __init__(self, **kwargs): - raise NotImplementedError("ChangeIncrementHandler is not yet implemented, pending finalization of the " - "SMIRNOFF spec") - def check_handler_compatibility(self, other_handler, assume_missing_is_default=True): @@ -3500,17 +3509,36 @@ def check_handler_compatibility(self, """ int_attrs_to_compare = ['number_of_conformers'] - string_attrs_to_compare = ['quantum_chemical_method', 'partial_charge_method'] + string_attrs_to_compare = ['partial_charge_method'] self._check_attributes_are_equal(other_handler, identical_attrs=string_attrs_to_compare+int_attrs_to_compare) + def find_matches(self, entity): + """Find the elements of the topology/molecule matched by a parameter type. - def create_force(self, system, topology, **kwargs): + Parameters + ---------- + entity : openforcefield.topology.Topology + Topology to search. + Returns + --------- + matches : ValenceDict[Tuple[int], ParameterHandler._Match] + ``matches[particle_indices]`` is the ``ParameterType`` object + matching the tuple of particle indices in ``entity``. + """ + # Using SortedDict here leads to the desired deduplication behavior, BUT it mangles the order + # of the atom indices in the keys. Thankfully, the Match objects that are values in `matches` contain + # `match.environment_match.topology_atom_indices`, which has the tuple in the correct order + matches = self._find_matches(entity, transformed_dict_cls=SortedDict) + return matches + def create_force(self, system, topology, **kwargs): + import warnings from openforcefield.topology import FrozenMolecule, TopologyAtom, TopologyVirtualSite + # We only want one instance of this force type existing = [system.getForce(i) for i in range(system.getNumForces())] existing = [f for f in existing if type(f) == self._OPENMMTYPE] if len(existing) == 0: @@ -3528,12 +3556,19 @@ def create_force(self, system, topology, **kwargs): # Make a temporary copy of ref_mol to assign charges from charge_mol temp_mol = FrozenMolecule(ref_mol) - # If the molecule wasn't assigned parameters from a manually-input charge_mol, calculate them here - temp_mol.generate_conformers(n_conformers=10) - temp_mol.compute_partial_charges(quantum_chemical_method=self._quantum_chemical_method, - partial_charge_method=self._partial_charge_method) + toolkit_registry = kwargs.get('toolkit_registry', GLOBAL_TOOLKIT_REGISTRY) + try: + # If the molecule wasn't assigned parameters from a manually-input charge_mol, calculate them here + temp_mol.generate_conformers(n_conformers=self.number_of_conformers) + temp_mol.assign_partial_charges(partial_charge_method=self.partial_charge_method, + toolkit_registry=toolkit_registry) + except Exception as e: + warnings.warn(str(e), Warning) + continue - # Assign charges to relevant atoms + charges_to_assign = {} + + # Assign initial, un-incremented charges to relevant atoms for topology_molecule in topology._reference_molecule_to_topology_molecules[ref_mol]: for topology_particle in topology_molecule.particles: topology_particle_index = topology_particle.topology_particle_index @@ -3542,52 +3577,37 @@ def create_force(self, system, topology, **kwargs): if type(topology_particle) is TopologyVirtualSite: ref_mol_particle_index = topology_particle.virtual_site.molecule_particle_index particle_charge = temp_mol._partial_charges[ref_mol_particle_index] - - # Retrieve nonbonded parameters for reference atom (charge not set yet) - _, sigma, epsilon = force.getParticleParameters(topology_particle_index) - # Set the nonbonded force with the partial charge - force.setParticleParameters(topology_particle_index, - particle_charge, sigma, - epsilon) + charges_to_assign[topology_particle_index] = particle_charge + + # Find SMARTS-based matches for charge increments + charge_increment_matches = self.find_matches(topology) + + # We ignore the atom index order in the keys here, since they have been + # sorted in order to deduplicate matches and let us identify when one parameter overwrites another + # in the SMIRNOFF parameter hierarchy. Since they are sorted, the position of the atom index + # in the key tuple DOES NOT correspond to the appropriate chargeincrementX value. + # Instead, the correct ordering of the match indices is found in + # charge_increment_match.environment_match.topology_atom_indices + for (_, charge_increment_match) in charge_increment_matches.items(): + # Adjust the values in the charges_to_assign dict by adding any + # charge increments onto the existing values + atom_indices = charge_increment_match.environment_match.topology_atom_indices + charge_increment = charge_increment_match.parameter_type + for top_particle_idx, charge_increment in zip(atom_indices, charge_increment.charge_increment): + if top_particle_idx in charges_to_assign: + charges_to_assign[top_particle_idx] += charge_increment + + # Set the incremented charges on the System particles + for topology_particle_index, charge_to_assign in charges_to_assign.items(): + _, sigma, epsilon = force.getParticleParameters(topology_particle_index) + force.setParticleParameters(topology_particle_index, + charge_to_assign, sigma, + epsilon) # Finally, mark that charges were assigned for this reference molecule self.mark_charges_assigned(ref_mol, topology) - - - # TODO: Move chargeModel and library residue charges to SMIRNOFF spec - def postprocess_system(self, system, topology, **kwargs): - bond_matches = self.find_matches(topology) - - # Apply bond charge increments to all appropriate force groups - # QUESTION: Should we instead apply this to the Topology in a preprocessing step, prior to spreading out charge onto virtual sites? - for force in system.getForces(): - if force.__class__.__name__ in [ - 'NonbondedForce' - ]: # TODO: We need to apply this to all Force types that involve charges, such as (Custom)GBSA forces and CustomNonbondedForce - for (atoms, bond_match) in bond_matches.items(): - bond = bond_match.parameter_type - - # Get corresponding particle indices in Topology - particle_indices = tuple( - [atom.particle_index for atom in atoms]) - # Retrieve parameters - [charge0, sigma0, epsilon0] = force.getParticleParameters( - particle_indices[0]) - [charge1, sigma1, epsilon1] = force.getParticleParameters( - particle_indices[1]) - # Apply bond charge increment - charge0 -= bond.increment - charge1 += bond.increment - # Update charges - force.setParticleParameters(particle_indices[0], charge0, - sigma0, epsilon0) - force.setParticleParameters(particle_indices[1], charge1, - sigma1, epsilon1) - # TODO: Calculate exceptions - - class GBSAHandler(ParameterHandler): """Handle SMIRNOFF ```` tags @@ -3605,12 +3625,11 @@ class GBSAType(ParameterType): radius = ParameterAttribute(unit=unit.angstrom) scale = ParameterAttribute(converter=float) - def __init__(self, **kwargs): - super().__init__(**kwargs) - _TAGNAME = 'GBSA' _INFOTYPE = GBSAType _OPENMMTYPE = openmm.GBSAOBCForce + # It's important that this runs AFTER partial charges are assigned to all particles, since this will need to + # collect and assign them to the GBSA particles _DEPENDENCIES = [vdWHandler, ElectrostaticsHandler, ToolkitAM1BCCHandler, ChargeIncrementModelHandler, LibraryChargeHandler] diff --git a/openforcefield/utils/toolkits.py b/openforcefield/utils/toolkits.py index acc4c5bf1..7d07b3435 100644 --- a/openforcefield/utils/toolkits.py +++ b/openforcefield/utils/toolkits.py @@ -61,7 +61,7 @@ from simtk import unit import numpy as np -from openforcefield.utils import all_subclasses, MessageException, inherit_docstrings +from openforcefield.utils import all_subclasses, MessageException, inherit_docstrings, temporary_cd #============================================================================================= @@ -121,6 +121,26 @@ class GAFFAtomTypeWarning(RuntimeWarning): pass +class ChargeMethodUnavailableError(MessageException): + """A toolkit does not support the requested partial_charge_method combination""" + pass + + +class IncorrectNumConformersError(MessageException): + """The requested partial_charge_method expects a different number of conformers than was provided""" + pass + + +class IncorrectNumConformersWarning(Warning): + """The requested partial_charge_method expects a different number of conformers than was provided""" + pass + + +class ChargeCalculationError(MessageException): + """An unhandled error occured in an external toolkit during charge calculation""" + pass + + #============================================================================================= # TOOLKIT UTILITY DECORATORS #============================================================================================= @@ -142,8 +162,6 @@ class ToolkitWrapper: _is_available = None # True if toolkit is available _toolkit_name = None # Name of the toolkit _toolkit_installation_instructions = None # Installation instructions for the toolkit - _toolkit_file_read_formats = None # The file types that this toolkit can read - _toolkit_file_write_formats = None # The file types that this toolkit can write #@staticmethod # TODO: Right now, to access the class definition, I have to make this a classmethod @@ -261,6 +279,168 @@ def from_file_obj(self, return NotImplementedError + def _check_n_conformers(self, + molecule, + partial_charge_method, + min_confs=None, + max_confs=None, + strict_n_conformers=False): + """ + Private method for validating the number of conformers on a molecule prior to partial + charge calculation + + Parameters + ---------- + molecule : Molecule + Molecule for which partial charges are to be computed + partial_charge_method : str, optional, default=None + The name of the charge method being used + min_confs : int, optional, default=None + The minimum number of conformers required to use this charge method + max_confs : int, optional, default=None + The maximum number of conformers required to use this charge method + strict_n_conformers : bool, default=False + Whether to raise an exception if an invalid number of conformers is provided. + If this is False and an invalid number of conformers is found, a warning will be raised. + + Raises + ------ + IncorrectNumConformersError + If the wrong number of conformers is attached to the input molecule, and strict_n_conformers is True. + """ + import warnings + n_confs = molecule.n_conformers + wrong_confs_msg = f"Molecule '{molecule}' has {n_confs} conformers, " \ + f"but charge method '{partial_charge_method}' expects" + exception_suffix = "You can disable this error by setting `strict_n_conformers=False' " \ + "when calling 'molecule.assign_partial_charges'." + # If there's no n_confs filter, then this molecule automatically passes + if min_confs is None and max_confs is None: + return + # If there's constraints on both ends, check both limits + elif min_confs is not None and max_confs is not None: + if not(min_confs <= n_confs <= max_confs): + if min_confs == max_confs: + wrong_confs_msg += f" exactly {min_confs}." + else: + wrong_confs_msg += f" between {min_confs} and {max_confs}." + + else: + return + # If there's only a max constraint, check that + elif min_confs is not None and max_confs is None: + if not(min_confs <= n_confs): + wrong_confs_msg += f" at least {min_confs}." + else: + return + # If there's only a maximum constraint, check that + elif min_confs is None and max_confs is not None: + if not(n_confs <= max_confs): + wrong_confs_msg += f" at most {max_confs}." + else: + return + # If we've made it this far, the molecule has the wrong number of conformers + if strict_n_conformers: + wrong_confs_msg += exception_suffix + raise IncorrectNumConformersError(wrong_confs_msg) + else: + warnings.warn(wrong_confs_msg, IncorrectNumConformersWarning) + + +@inherit_docstrings +class BuiltInToolkitWrapper(ToolkitWrapper): + """ + Built-in ToolkitWrapper for very basic functionality. This is intended for use in testing and not much more. + + .. warning :: This API is experimental and subject to change. + """ + _toolkit_name = 'Built-in Toolkit' + _toolkit_installation_instructions = 'This toolkit is installed with the Open Force Field Toolkit and does ' \ + 'not require additional dependencies.' + + def __init__(self): + super().__init__() + + self._toolkit_file_read_formats = [] + self._toolkit_file_write_formats = [] + + + def assign_partial_charges(self, + molecule, + partial_charge_method=None, + use_conformers=None, + strict_n_conformers=False): + """ + Compute partial charges with the built-in toolkit using simple arithmetic operations, and assign + the new values to the partial_charges attribute. + + .. warning :: This API is experimental and subject to change. + + Parameters + ---------- + molecule : openforcefield.topology.Molecule + Molecule for which partial charges are to be computed + partial_charge_method: str, optional, default=None + The charge model to use. One of ['zeros', 'formal_charge']. If None, 'formal_charge' will be used. + use_conformers : iterable of simtk.unit.Quantity-wrapped numpy arrays, each with shape (n_atoms, 3) and dimension of distance. Optional, default = None + Coordinates to use for partial charge calculation. If None, an appropriate number of conformers + will be generated. + strict_n_conformers : bool, default=False + Whether to raise an exception if an invalid number of conformers is provided for the given charge method. + If this is False and an invalid number of conformers is found, a warning will be raised + instead of an Exception. + + Raises + ------ + ChargeMethodUnavailableError if the requested charge method can not be handled by this toolkit + + IncorrectNumConformersError if strict_n_conformers is True and use_conformers is provided and specifies an + invalid number of conformers for the requested method + + ChargeCalculationError if the charge calculation is supported by this toolkit, but fails + """ + from openforcefield.topology import Molecule + PARTIAL_CHARGE_METHODS = {'zeros': {'rec_confs':0, + 'min_confs':0, + 'max_confs':0}, + 'formal_charge': {'rec_confs':0, + 'min_confs':0, + 'max_confs':0} + } + + if partial_charge_method is None: + partial_charge_method = 'formal_charge' + + # Make a temporary copy of the molecule, since we'll be messing with its conformers + mol_copy = Molecule(molecule) + + partial_charge_method = partial_charge_method.lower() + if partial_charge_method not in PARTIAL_CHARGE_METHODS: + raise ChargeMethodUnavailableError(f'Partial charge method "{partial_charge_method}"" is not supported by ' + f'the Built-in toolkit. Available charge methods are ' + f'{list(PARTIAL_CHARGE_METHODS.keys())}') + + if use_conformers is None: + # Note that this refers back to the GLOBAL_TOOLKIT_REGISTRY by default, since + # BuiltInToolkitWrapper can't generate conformers + mol_copy.generate_conformers(n_conformers=PARTIAL_CHARGE_METHODS[partial_charge_method]['rec_confs']) + else: + mol_copy._conformers = None + for conformer in use_conformers: + mol_copy.add_conformer(conformer) + self._check_n_conformers(mol_copy, partial_charge_method=partial_charge_method, min_confs=0, + max_confs=0,strict_n_conformers=strict_n_conformers) + + partial_charges = unit.Quantity(np.zeros((molecule.n_particles)), unit.elementary_charge) + if partial_charge_method == 'zeroes': + pass + elif partial_charge_method == 'formal_charge': + for part_idx, particle in enumerate(molecule.particles): + partial_charges[part_idx] = particle.formal_charge + + molecule.partial_charges = partial_charges + + @inherit_docstrings class OpenEyeToolkitWrapper(ToolkitWrapper): """ @@ -350,7 +530,6 @@ def is_available( def from_object(self, object, allow_undefined_stereo=False): """ If given an OEMol (or OEMol-derived object), this function will load it into an openforcefield.topology.molecule - Otherwise, it will return False. Parameters ---------- @@ -1008,7 +1187,7 @@ def describe_oeatom(oeatom): oe_idx = oeatom.GetIdx() map_id = oeatom.GetMapIdx() atomic_number = oeatom.GetAtomicNum() - formal_charge = oeatom.GetFormalCharge() + formal_charge = oeatom.GetFormalCharge() * unit.elementary_charge is_aromatic = oeatom.IsAromatic() stereochemistry = OpenEyeToolkitWrapper._openeye_cip_atom_stereochemistry( oemol, oeatom) @@ -1150,7 +1329,7 @@ def to_openeye(molecule, aromaticity_model=DEFAULT_AROMATICITY_MODEL): oemol_atoms = list() # list of corresponding oemol atoms for atom in molecule.atoms: oeatom = oemol.NewAtom(atom.atomic_number) - oeatom.SetFormalCharge(atom.formal_charge) + oeatom.SetFormalCharge(atom.formal_charge.value_in_unit(unit.elementary_charge)) # simtk.unit.Quantity(1, unit.elementary_charge) oeatom.SetAromatic(atom.is_aromatic) oeatom.SetData('name', atom.name) oeatom.SetPartialCharge(float('nan')) @@ -1548,7 +1727,7 @@ def from_inchi(self, inchi, allow_undefined_stereo=False): return molecule - def generate_conformers(self, molecule, n_conformers=1, clear_existing=True): + def generate_conformers(self, molecule, n_conformers=1, rms_cutoff=None, clear_existing=True): """ Generate molecule conformers using OpenEye Omega. @@ -1566,6 +1745,9 @@ def generate_conformers(self, molecule, n_conformers=1, clear_existing=True): The molecule to generate conformers for. n_conformers : int, default=1 The maximum number of conformers to generate. + rms_cutoff : simtk.Quantity-wrapped float, in units of distance, optional, default=None + The minimum RMS value at which two conformers are considered redundant and one is deleted. + If None, the cutoff is set to 1 Angstrom clear_existing : bool, default=True Whether to overwrite existing conformers for the molecule @@ -1577,8 +1759,11 @@ def generate_conformers(self, molecule, n_conformers=1, clear_existing=True): omega.SetCanonOrder(False) omega.SetSampleHydrogens(True) omega.SetEnergyWindow(15.0) #unit? - omega.SetRMSThreshold(1.0) - #Don't generate random stereoisomer if not specified + if rms_cutoff is None: + omega.SetRMSThreshold(1.0) + else: + omega.SetRMSThreshold(rms_cutoff.value_in_unit(unit.angstrom)) + # Don't generate random stereoisomer if not specified omega.SetStrictStereo(True) status = omega(oemol) @@ -1597,9 +1782,14 @@ def generate_conformers(self, molecule, n_conformers=1, clear_existing=True): for conformer in molecule2._conformers: molecule._add_conformer(conformer) - def compute_partial_charges(self, molecule, quantum_chemical_method="AM1-BCC", partial_charge_method='None'): + def assign_partial_charges(self, + molecule, + partial_charge_method=None, + use_conformers=None, + strict_n_conformers=False): """ - Compute partial charges with OpenEye quacpac + Compute partial charges with OpenEye quacpac, and assign + the new values to the partial_charges attribute. .. warning :: This API is experimental and subject to change. @@ -1611,130 +1801,113 @@ def compute_partial_charges(self, molecule, quantum_chemical_method="AM1-BCC", p Parameters ---------- - molecule : Molecule + molecule : openforcefield.topology.Molecule Molecule for which partial charges are to be computed - charge_model : str, optional, default=None - The charge model to use. One of ['noop', 'mmff', 'mmff94', 'am1bcc', 'am1bccnosymspt', 'amber', - 'amberff94', 'am1bccelf10'] - If None, 'am1bcc' will be used. - - Returns - ------- - charges : numpy.array of shape (natoms) of type float - The partial charges - - """ - raise NotImplementedError - # TODO: Implement this in a way that's compliant with SMIRNOFF's tag when the spec gets finalized - - # from openeye import oequacpac - # import numpy as np - # - # if molecule.n_conformers == 0: - # raise Exception( - # "No conformers present in molecule submitted for partial charge calculation. Consider " - # "loading the molecule from a file with geometry already present or running " - # "molecule.generate_conformers() before calling molecule.compute_partial_charges" - # ) - # oemol = molecule.to_openeye() - # - # ## This seems like a big decision. Implemented a simple solution here. Not to be considered final. - # ## Some discussion at https://github.com/openforcefield/openforcefield/pull/86#issuecomment-350111236 - # - # if charge_model is None: - # charge_model = "am1bcc" - # - # if charge_model == "noop": - # result = oequacpac.OEAssignCharges(oemol, - # oequacpac.OEChargeEngineNoOp()) - # elif charge_model == "mmff" or charge_model == "mmff94": - # result = oequacpac.OEAssignCharges(oemol, - # oequacpac.OEMMFF94Charges()) - # elif charge_model == "am1bcc": - # result = oequacpac.OEAssignCharges(oemol, - # oequacpac.OEAM1BCCCharges()) - # elif charge_model == "am1bccnosymspt": - # optimize = True - # symmetrize = True - # result = oequacpac.OEAssignCharges( - # oemol, oequacpac.OEAM1BCCCharges(not optimize, not symmetrize)) - # elif charge_model == "amber" or charge_model == "amberff94": - # result = oequacpac.OEAssignCharges(oemol, - # oequacpac.OEAmberFF94Charges()) - # elif charge_model == "am1bccelf10": - # result = oequacpac.OEAssignCharges( - # oemol, oequacpac.OEAM1BCCELF10Charges()) - # else: - # raise ValueError('charge_model {} unknown'.format(charge_model)) - # - # if result is False: - # raise Exception('Unable to assign charges') - # - # # Extract and return charges - # ## TODO: Behavior when given multiple conformations? - # ## TODO: Make sure atom mapping remains constant - # - # charges = unit.Quantity( - # np.zeros([oemol.NumAtoms()], np.float64), unit.elementary_charge) - # for index, atom in enumerate(oemol.GetAtoms()): - # charge = atom.GetPartialCharge() - # charge = charge * unit.elementary_charge - # charges[index] = charge - # - # if ((charges / unit.elementary_charge) == 0. - # ).all() and not (charge_model == 'noop'): - # # TODO: These will be 0 if the charging failed. What behavior do we want in that case? - # raise Exception( - # "Partial charge calculation failed. Charges from compute_partial_charges() are all 0." - # ) - # molecule.set_partial_charges(charges) - + partial_charge_method : str, optional, default=None + The charge model to use. One of ['amberff94', 'mmff', 'mmff94', `am1-mulliken`, 'am1bcc', + 'am1bccnosymspt', 'am1bccelf10'] + If None, 'am1-mulliken' will be used. + use_conformers : iterable of simtk.unit.Quantity-wrapped numpy arrays, each with shape (n_atoms, 3) and dimension of distance. Optional, default = None + Coordinates to use for partial charge calculation. If None, an appropriate number of conformers will be generated. + strict_n_conformers : bool, default=False + Whether to raise an exception if an invalid number of conformers is provided for the given charge method. + If this is False and an invalid number of conformers is found, a warning will be raised. + Raises + ------ + ChargeMethodUnavailableError if the requested charge method can not be handled by this toolkit - def compute_partial_charges_am1bcc(self, molecule): + ChargeCalculationError if the charge method is supported by this toolkit, but fails """ - Compute AM1BCC partial charges with OpenEye quacpac. This function will attempt to use - the OEAM1BCCELF10 charge generation method, but may print a warning and fall back to - normal OEAM1BCC if an error is encountered. This error is known to occur with some - carboxylic acids, and is under investigation by OpenEye. - - - .. warning :: This API is experimental and subject to change. - .. todo :: + from openeye import oequacpac, oechem + import numpy as np + from openforcefield.topology import Molecule - * Should the default be ELF? - * Can we expose more charge models? + SUPPORTED_CHARGE_METHODS = {'am1bcc': {'oe_charge_method': oequacpac.OEAM1BCCCharges, + 'min_confs': 1, + 'max_confs': 1, + 'rec_confs': 1, + }, + 'am1-mulliken': {'oe_charge_method': oequacpac.OEAM1Charges, + 'min_confs': 1, + 'max_confs': 1, + 'rec_confs': 1, + }, + 'gasteiger': {'oe_charge_method': oequacpac.OEGasteigerCharges, + 'min_confs': 0, + 'max_confs': 0, + 'rec_confs': 0, + }, + 'mmff94': {'oe_charge_method': oequacpac.OEMMFF94Charges, + 'min_confs': 0, + 'max_confs': 0, + 'rec_confs': 0, + }, + 'am1bccnosymspt': {'oe_charge_method': oequacpac.OEAM1BCCCharges, + 'min_confs': 1, + 'max_confs': 1, + 'rec_confs': 1, + }, + 'am1bccelf10': {'oe_charge_method': oequacpac.OEAM1BCCELF10Charges, + 'min_confs': 1, + 'max_confs': None, + 'rec_confs': 500, + }, + } + + if partial_charge_method is None: + partial_charge_method = "am1-mulliken" + + partial_charge_method = partial_charge_method.lower() + + if partial_charge_method not in SUPPORTED_CHARGE_METHODS: + raise ChargeMethodUnavailableError( + f"partial_charge_method '{partial_charge_method}' is not available from OpenEyeToolkitWrapper. " + f"Available charge methods are {list(SUPPORTED_CHARGE_METHODS.keys())} " + ) + charge_method = SUPPORTED_CHARGE_METHODS[partial_charge_method] - Parameters - ---------- - molecule : Molecule - Molecule for which partial charges are to be computed + # Make a temporary copy of the molecule, since we'll be messing with its conformers + mol_copy = Molecule(molecule) - Returns - ------- - charges : numpy.array of shape (natoms) of type float - The partial charges - - """ - from openeye import oequacpac, oechem - import numpy as np + if use_conformers is None: + if charge_method['rec_confs'] == 0: + mol_copy._conformers = None + else: + self.generate_conformers(mol_copy, + n_conformers=charge_method['rec_confs'], + rms_cutoff=0.25*unit.angstrom) + # TODO: What's a "best practice" RMS cutoff to use here? + else: + mol_copy._conformers = None + for conformer in use_conformers: + mol_copy.add_conformer(conformer) + self._check_n_conformers(mol_copy, + partial_charge_method=partial_charge_method, + min_confs=charge_method['min_confs'], + max_confs=charge_method['max_confs'], + strict_n_conformers=strict_n_conformers) - if molecule.n_conformers == 0: - raise Exception( - "No conformers present in molecule submitted for partial charge calculation. Consider " - "loading the molecule from a file with geometry already present or running " - "molecule.generate_conformers() before calling molecule.compute_partial_charges" - ) - oemol = self.to_openeye(molecule) + oemol = mol_copy.to_openeye() errfs = oechem.oeosstream() oechem.OEThrow.SetOutputStream(errfs) oechem.OEThrow.Clear() - quacpac_status = oequacpac.OEAssignCharges(oemol, oequacpac.OEAM1BCCELF10Charges()) + + # The OpenFF toolkit has always supported a version of AM1BCC with no geometry optimization + # or symmetry correction. So we include this keyword to provide a special configuration of quacpac + # if requested. + if partial_charge_method == "am1bccnosymspt": + optimize = False + symmetrize = False + quacpac_status = oequacpac.OEAssignCharges(oemol, charge_method['oe_charge_method'](optimize, symmetrize)) + else: + quacpac_status = oequacpac.OEAssignCharges(oemol, charge_method['oe_charge_method']()) + oechem.OEThrow.SetOutputStream(oechem.oeerr) # restoring to original state - # This logic handles errors encountered in #34 + # This logic handles errors encountered in #34, which can occur when using ELF10 conformer selection if not quacpac_status: if "SelectElfPop: issue with removing trans COOH conformers" in (errfs.str().decode("UTF-8")): logger.warning("Warning: OEAM1BCCELF10 charge assignment failed due to a known bug (toolkit issue " @@ -1743,7 +1916,7 @@ def compute_partial_charges_am1bcc(self, molecule): quacpac_status = oequacpac.OEAssignCharges(oemol, oequacpac.OEAM1BCCCharges()) if quacpac_status is False: - raise Exception('Unable to assign charges; OE error: "{}"'.format(errfs.str().decode("UTF-8"))) + raise ChargeCalculationError(f'Unable to assign charges: {errfs.str().decode("UTF-8")}') # Extract and return charges ## TODO: Make sure atom mapping remains constant @@ -1756,7 +1929,46 @@ def compute_partial_charges_am1bcc(self, molecule): charge = charge * unit.elementary_charge charges[index] = charge - return charges + molecule.partial_charges = charges + + + def compute_partial_charges_am1bcc(self, molecule, use_conformers=None, strict_n_conformers=False): + """ + Compute AM1BCC partial charges with OpenEye quacpac. This function will attempt to use + the OEAM1BCCELF10 charge generation method, but may print a warning and fall back to + normal OEAM1BCC if an error is encountered. This error is known to occur with some + carboxylic acids, and is under investigation by OpenEye. + + + .. warning :: This API is experimental and subject to change. + + Parameters + ---------- + molecule : Molecule + Molecule for which partial charges are to be computed + use_conformers : iterable of simtk.unit.Quantity-wrapped numpy arrays, each with shape (n_atoms, 3) and dimension of distance. Optional, default = None + Coordinates to use for partial charge calculation. If None, an appropriate number of conformers + will be generated. + strict_n_conformers : bool, default=False + Whether to raise an exception if an invalid number of conformers is provided. + If this is False and an invalid number of conformers is found, a warning will be raised + instead of an Exception. + + Returns + ------- + charges : numpy.array of shape (natoms) of type float + The partial charges + """ + + import warnings + warnings.warn("compute_partial_charges_am1bcc will be deprecated in an upcoming release. " + "Use assign_partial_charges(partial_charge_method='am1bccelf10') instead.", + DeprecationWarning) + self.assign_partial_charges(molecule, + partial_charge_method='am1bccelf10', + use_conformers=use_conformers, + strict_n_conformers=strict_n_conformers) + return molecule.partial_charges def assign_fractional_bond_orders(self, molecule, bond_order_model=None, use_conformers=None): """ @@ -2598,7 +2810,7 @@ def from_inchi(self, inchi, allow_undefined_stereo=False): return molecule - def generate_conformers(self, molecule, n_conformers=1, clear_existing=True): + def generate_conformers(self, molecule, n_conformers=1, rms_cutoff=None, clear_existing=True): """ Generate molecule conformers using RDKit. @@ -2615,18 +2827,24 @@ def generate_conformers(self, molecule, n_conformers=1, clear_existing=True): The molecule to generate conformers for. n_conformers : int, default=1 Maximum number of conformers to generate. + rms_cutoff : simtk.Quantity-wrapped float, in units of distance, optional, default=None + The minimum RMS value at which two conformers are considered redundant and one is deleted. + If None, the cutoff is set to 1 Angstrom + clear_existing : bool, default=True Whether to overwrite existing conformers for the molecule. """ from rdkit.Chem import AllChem + if rms_cutoff is None: + rms_cutoff = 1. * unit.angstrom rdmol = self.to_rdkit(molecule) # TODO: This generates way more conformations than omega, given the same nConfs and RMS threshold. Is there some way to set an energy cutoff as well? AllChem.EmbedMultipleConfs( rdmol, numConfs=n_conformers, - pruneRmsThresh=1.0, + pruneRmsThresh=rms_cutoff / unit.angstrom, randomSeed=1, #params=AllChem.ETKDG() ) @@ -2731,7 +2949,7 @@ def from_rdkit(self, rdmol, allow_undefined_stereo=False): # create a new atom #atomic_number = oemol.NewAtom(rda.GetAtomicNum()) atomic_number = rda.GetAtomicNum() - formal_charge = rda.GetFormalCharge() + formal_charge = rda.GetFormalCharge() * unit.elementary_charge is_aromatic = rda.GetIsAromatic() if rda.HasProp('_Name'): name = rda.GetProp('_Name') @@ -2838,7 +3056,7 @@ def from_rdkit(self, rdmol, allow_undefined_stereo=False): else: # If some other atoms had partial charges but this one doesn't, raise an Exception if any_atom_has_partial_charge: - raise Exception( + raise ValueError( "Some atoms in rdmol have partial charges, but others do not." ) if any_atom_has_partial_charge: @@ -2913,7 +3131,7 @@ def to_rdkit(cls, molecule, aromaticity_model=DEFAULT_AROMATICITY_MODEL): for index, atom in enumerate(molecule.atoms): rdatom = Chem.Atom(atom.atomic_number) - rdatom.SetFormalCharge(atom.formal_charge) + rdatom.SetFormalCharge(atom.formal_charge.value_in_unit(unit.elementary_charge)) rdatom.SetIsAromatic(atom.is_aromatic) rdatom.SetProp('_Name', atom.name) @@ -3550,12 +3768,19 @@ def is_available(): ANTECHAMBER_PATH = find_executable("antechamber") if ANTECHAMBER_PATH is None: return False - else: - return True + # AmberToolsToolkitWrapper needs RDKit to do basically anything, since its interface requires SDF I/O + if not(RDKitToolkitWrapper.is_available()): + return False + return True - def compute_partial_charges(self, molecule, charge_model=None): + def assign_partial_charges(self, + molecule, + partial_charge_method=None, + use_conformers=None, + strict_n_conformers=False): """ - Compute partial charges with AmberTools using antechamber/sqm + Compute partial charges with AmberTools using antechamber/sqm, and assign + the new values to the partial_charges attribute. .. warning :: This API experimental and subject to change. @@ -3565,180 +3790,160 @@ def compute_partial_charges(self, molecule, charge_model=None): Parameters ---------- - molecule : Molecule + molecule : openforcefield.topology.Molecule Molecule for which partial charges are to be computed - charge_model : str, optional, default=None - The charge model to use. One of ['gas', 'mul', 'bcc']. If None, 'bcc' will be used. - + partial_charge_method : str, optional, default=None + The charge model to use. One of ['gasteiger', 'am1bcc', 'am1-mulliken']. If None, 'am1-mulliken' will be used. + use_conformers : iterable of simtk.unit.Quantity-wrapped numpy arrays, each with shape (n_atoms, 3) and dimension of distance. Optional, default = None + List of (n_atoms x 3) simtk.unit.Quantities to use for partial charge calculation. + If None, an appropriate number of conformers will be generated. + strict_n_conformers : bool, default=False + Whether to raise an exception if an invalid number of conformers is provided for the given charge method. + If this is False and an invalid number of conformers is found, a warning will be raised. Raises ------ - ValueError if the requested charge method could not be handled - - Notes - ----- - Currently only sdf file supported as input and mol2 as output - https://github.com/choderalab/openmoltools/blob/master/openmoltools/packmol.py + ChargeMethodUnavailableError if the requested charge method can not be handled by this toolkit + ChargeCalculationError if the charge method is supported by this toolkit, but fails """ - raise NotImplementedError - # TODO: Implement this in a way that's compliant with SMIRNOFF's tag when the spec gets finalized - # import os - # from simtk import unit - # - # if charge_model is None: - # charge_model = 'bcc' - # - # # Check that the requested charge method is supported - # # Needs to be fixed: 'cm1', 'cm2', - # SUPPORTED_ANTECHAMBER_CHARGE_MODELS = ['gas', 'mul', 'bcc'] - # if charge_model not in SUPPORTED_ANTECHAMBER_CHARGE_MODELS: - # raise ValueError( - # 'Requested charge method {} not among supported charge ' - # 'methods {}'.format(charge_model, - # SUPPORTED_ANTECHAMBER_CHARGE_MODELS)) - # - # # Find the path to antechamber - # # TODO: How should we implement find_executable? - # ANTECHAMBER_PATH = find_executable("antechamber") - # if ANTECHAMBER_PATH is None: - # raise (IOError("Antechamber not found, cannot run charge_mol()")) - # - # if len(molecule._conformers) == 0: - # raise Exception( - # "No conformers present in molecule submitted for partial charge calculation. Consider " - # "loading the molecule from a file with geometry already present or running " - # "molecule.generate_conformers() before calling molecule.compute_partial_charges" - # ) - # - # - # # Compute charges - # import tempfile - # with tempfile.TemporaryDirectory() as tmpdir: - # net_charge = molecule.total_charge - # # Write out molecule in SDF format - # ## TODO: How should we handle multiple conformers? - # self._rdkit_toolkit_wrapper.to_file( - # molecule, 'molecule.sdf', file_format='sdf') - # #os.system('ls') - # #os.system('cat molecule.sdf') - # # Compute desired charges - # # TODO: Add error handling if antechamber chokes - # # TODO: Add something cleaner than os.system - # os.system( - # "antechamber -i molecule.sdf -fi sdf -o charged.mol2 -fo mol2 -pf " - # "yes -c {} -nc {}".format(charge_model, net_charge)) - # #os.system('cat charged.mol2') - # - # # Write out just charges - # os.system( - # "antechamber -i charged.mol2 -fi mol2 -o charges2.mol2 -fo mol2 -c wc " - # "-cf charges.txt -pf yes") - # #os.system('cat charges.txt') - # # Check to ensure charges were actually produced - # if not os.path.exists('charges.txt'): - # # TODO: copy files into local directory to aid debugging? - # raise Exception( - # "Antechamber/sqm partial charge calculation failed on " - # "molecule {} (SMILES {})".format( - # molecule.name, molecule.to_smiles())) - # # Read the charges - # with open('charges.txt', 'r') as infile: - # contents = infile.read() - # text_charges = contents.split() - # charges = np.zeros([molecule.n_atoms], np.float64) - # for index, token in enumerate(text_charges): - # charges[index] = float(token) - # # TODO: Ensure that the atoms in charged.mol2 are in the same order as in molecule.sdf - # - # charges = unit.Quantity(charges, unit.elementary_charge) - # - # molecule.set_partial_charges(charges) + import os + import subprocess + from openforcefield.topology import Molecule + + if partial_charge_method is None: + partial_charge_method = 'am1-mulliken' + else: + # Standardize method name for string comparisons + partial_charge_method = partial_charge_method.lower() + + SUPPORTED_CHARGE_METHODS = {'am1bcc': {'antechamber_keyword':'bcc', + 'min_confs': 1, + 'max_confs': 1, + 'rec_confs': 1, + }, + 'am1-mulliken': {'antechamber_keyword': 'mul', + 'min_confs': 1, + 'max_confs': 1, + 'rec_confs': 1, + }, + 'gasteiger': {'antechamber_keyword': 'gas', + 'min_confs': 0, + 'max_confs': 0, + 'rec_confs': 0, + } + } + + + if partial_charge_method not in SUPPORTED_CHARGE_METHODS: + raise ChargeMethodUnavailableError( + f"partial_charge_method '{partial_charge_method}' is not available from AmberToolsToolkitWrapper. " + f"Available charge methods are {list(SUPPORTED_CHARGE_METHODS.keys())} " + ) + + charge_method = SUPPORTED_CHARGE_METHODS[partial_charge_method] + + # Make a temporary copy of the molecule, since we'll be messing with its conformers + mol_copy = Molecule(molecule) + + if use_conformers is None: + if charge_method['rec_confs'] == 0: + mol_copy._conformers = None + else: + mol_copy.generate_conformers(n_conformers=charge_method['rec_confs'], + rms_cutoff=0.25*unit.angstrom, + toolkit_registry=RDKitToolkitWrapper()) + # TODO: What's a "best practice" RMS cutoff to use here? + else: + mol_copy._conformers = None + for conformer in use_conformers: + mol_copy.add_conformer(conformer) + self._check_n_conformers(mol_copy, + partial_charge_method=partial_charge_method, + min_confs=charge_method['min_confs'], + max_confs=charge_method['max_confs'], + strict_n_conformers=strict_n_conformers) - def compute_partial_charges_am1bcc(self, molecule): + # Find the path to antechamber + # TODO: How should we implement find_executable? + ANTECHAMBER_PATH = find_executable("antechamber") + if ANTECHAMBER_PATH is None: + raise IOError("Antechamber not found, cannot run charge_mol()") + + # Compute charges + with tempfile.TemporaryDirectory() as tmpdir: + with temporary_cd(tmpdir): + net_charge = mol_copy.total_charge / unit.elementary_charge + # Write out molecule in SDF format + ## TODO: How should we handle multiple conformers? + self._rdkit_toolkit_wrapper.to_file( + mol_copy, 'molecule.sdf', file_format='sdf') + # Compute desired charges + # TODO: Add error handling if antechamber chokes + # TODO: Add something cleaner than os.system + short_charge_method = charge_method['antechamber_keyword'] + subprocess.check_output([ + "antechamber", "-i", "molecule.sdf", "-fi", "sdf", "-o", "charged.mol2", "-fo", + "mol2", "-pf", "yes", "-dr", "n", "-c", short_charge_method, "-nc", str(net_charge)] + ) + # Write out just charges + subprocess.check_output([ + "antechamber", "-dr", "n", "-i", "charged.mol2", "-fi", "mol2", "-o", "charges2.mol2", "-fo", + "mol2", "-c", "wc", "-cf", "charges.txt", "-pf", "yes"]) + # Check to ensure charges were actually produced + if not os.path.exists('charges.txt'): + # TODO: copy files into local directory to aid debugging? + raise ChargeCalculationError( + "Antechamber/sqm partial charge calculation failed on " + "molecule {} (SMILES {})".format( + molecule.name, molecule.to_smiles())) + # Read the charges + with open('charges.txt', 'r') as infile: + contents = infile.read() + text_charges = contents.split() + charges = np.zeros([molecule.n_atoms], np.float64) + for index, token in enumerate(text_charges): + charges[index] = float(token) + # TODO: Ensure that the atoms in charged.mol2 are in the same order as in molecule.sdf + charges = unit.Quantity(charges, unit.elementary_charge) + molecule.partial_charges = charges + + def compute_partial_charges_am1bcc(self, molecule, use_conformers=None, strict_n_conformers=False): """ Compute partial charges with AmberTools using antechamber/sqm. This will calculate AM1-BCC charges on the first conformer only. - .. warning :: This API experimental and subject to change. + .. warning :: This API is experimental and subject to change. Parameters ---------- molecule : Molecule Molecule for which partial charges are to be computed + use_conformers : iterable of simtk.unit.Quantity-wrapped numpy arrays, each with shape (n_atoms, 3) and dimension of distance. Optional, default = None + Coordinates to use for partial charge calculation. If None, an appropriate number of conformers + will be generated. + strict_n_conformers : bool, default=False + Whether to raise an exception if an invalid number of conformers is provided. + If this is False and an invalid number of conformers is found, a warning will be raised + instead of an Exception. - - Raises - ------ - ValueError if the requested charge method could not be handled - + Returns + ------- + charges : numpy.array of shape (natoms) of type float + The partial charges """ - import os - from simtk import unit + import warnings + warnings.warn("compute_partial_charges_am1bcc will be deprecated in an upcoming release. " + "Use assign_partial_charges(partial_charge_method='am1bcc') instead.", + DeprecationWarning) - # Find the path to antechamber - # TODO: How should we implement find_executable? - ANTECHAMBER_PATH = find_executable("antechamber") - if ANTECHAMBER_PATH is None: - raise (IOError("Antechamber not found, cannot run " - "AmberToolsToolkitWrapper.compute_partial_charges_am1bcc()")) - - if len(molecule._conformers) == 0: - raise ValueError( - "No conformers present in molecule submitted for partial charge calculation. Consider " - "loading the molecule from a file with geometry already present or running " - "molecule.generate_conformers() before calling molecule.compute_partial_charges" - ) - if len(molecule._conformers) > 1: - logger.warning("Warning: In AmberToolsToolkitwrapper.compute_partial_charges_am1bcc: " - "Molecule '{}' has more than one conformer, but this function " - "will only generate charges for the first one.".format(molecule.name)) - - - # Compute charges - with tempfile.TemporaryDirectory() as tmpdir: - net_charge = molecule.total_charge - # Write out molecule in SDF format - ## TODO: How should we handle multiple conformers? - self._rdkit_toolkit_wrapper.to_file( - molecule, 'molecule.sdf', file_format='sdf') - #os.system('ls') - #os.system('cat molecule.sdf') - # Compute desired charges - # TODO: Add error handling if antechamber chokes - # TODO: Add something cleaner than os.system - subprocess.check_output([ - "antechamber", "-i", "molecule.sdf", "-fi", "sdf", "-o", "charged.mol2", "-fo", "mol2", "-pf", "yes", "-c", "bcc", - "-nc", str(net_charge) - ]) - #os.system('cat charged.mol2') - - # Write out just charges - subprocess.check_output([ - "antechamber", "-i", "charged.mol2", "-fi", "mol2", "-o", "charges2.mol2", "-fo", "mol2", "-c", "wc", "-cf", - "charges.txt", "-pf", "yes" - ]) - #os.system('cat charges.txt') - # Check to ensure charges were actually produced - if not os.path.exists('charges.txt'): - # TODO: copy files into local directory to aid debugging? - raise Exception( - "Antechamber/sqm partial charge calculation failed on " - "molecule {} (SMILES {})".format( - molecule.name, molecule.to_smiles())) - # Read the charges - with open('charges.txt', 'r') as infile: - contents = infile.read() - text_charges = contents.split() - charges = np.zeros([molecule.n_atoms], np.float64) - for index, token in enumerate(text_charges): - charges[index] = float(token) - # TODO: Ensure that the atoms in charged.mol2 are in the same order as in molecule.sdf - - charges = unit.Quantity(charges, unit.elementary_charge) - return charges + self.assign_partial_charges(molecule, + partial_charge_method="AM1BCC", + use_conformers=use_conformers, + strict_n_conformers=strict_n_conformers) + return molecule.partial_charges def _modify_sqm_in_to_request_bond_orders(self, file_path): """ @@ -3908,25 +4113,26 @@ def assign_fractional_bond_orders(self, molecule, bond_order_model=None, use_con ac_charge_keyword = bond_order_model_to_antechamber_keyword[bond_order_model] with tempfile.TemporaryDirectory() as tmpdir: - net_charge = temp_mol.total_charge - # Write out molecule in SDF format - self._rdkit_toolkit_wrapper.to_file( - temp_mol, 'molecule.sdf', file_format='sdf') - # Prepare sqm.in file as if we were going to run charge calc - # TODO: Add error handling if antechamber chokes - subprocess.check_output([ - "antechamber", "-i", "molecule.sdf", "-fi", "sdf", "-o", - "sqm.in", "-fo", "sqmcrt", "-pf", "yes", "-c", ac_charge_keyword, - "-nc", str(net_charge) - ]) - # Modify sqm.in to request bond order calculation - self._modify_sqm_in_to_request_bond_orders("sqm.in") - # Run sqm to get bond orders - subprocess.check_output(["sqm", "-i", "sqm.in", "-o", "sqm.out", "-O"]) - # Ensure that antechamber/sqm did not change the indexing by checking against - # an ordered list of element symbols for this molecule - expected_elements = [at.element.symbol for at in molecule.atoms] - bond_orders = self._get_fractional_bond_orders_from_sqm_out('sqm.out', + with temporary_cd(tmpdir): + net_charge = temp_mol.total_charge + # Write out molecule in SDF format + self._rdkit_toolkit_wrapper.to_file( + temp_mol, 'molecule.sdf', file_format='sdf') + # Prepare sqm.in file as if we were going to run charge calc + # TODO: Add error handling if antechamber chokes + subprocess.check_output([ + "antechamber", "-i", "molecule.sdf", "-fi", "sdf", "-o", + "sqm.in", "-fo", "sqmcrt", "-pf", "yes", "-c", ac_charge_keyword, + "-nc", str(net_charge) + ]) + # Modify sqm.in to request bond order calculation + self._modify_sqm_in_to_request_bond_orders("sqm.in") + # Run sqm to get bond orders + subprocess.check_output(["sqm", "-i", "sqm.in", "-o", "sqm.out", "-O"]) + # Ensure that antechamber/sqm did not change the indexing by checking against + # an ordered list of element symbols for this molecule + expected_elements = [at.element.symbol for at in molecule.atoms] + bond_orders = self._get_fractional_bond_orders_from_sqm_out('sqm.out', validate_elements=expected_elements) # Note that sqm calculate WBOs for ALL PAIRS of atoms, not just those that have @@ -4011,7 +4217,7 @@ def __init__(self, if toolkit_precedence is None: toolkit_precedence = [ OpenEyeToolkitWrapper, RDKitToolkitWrapper, - AmberToolsToolkitWrapper + AmberToolsToolkitWrapper, BuiltInToolkitWrapper ] if register_imported_toolkit_wrappers: @@ -4154,7 +4360,7 @@ def resolve(self, method_name): raise NotImplementedError(msg) # TODO: Can we instead register available methods directly with `ToolkitRegistry`, so we can just use `ToolkitRegistry.method()`? - def call(self, method_name, *args, **kwargs): + def call(self, method_name, *args, raise_exception_types=None, **kwargs): """ Execute the requested method by attempting to use all registered toolkits in order of precedence. @@ -4166,11 +4372,23 @@ def call(self, method_name, *args, **kwargs): ---------- method_name : str The name of the method to execute + raise_exception_types : list of Exception subclasses, default=None + A list of exception-derived types to catch and raise immediately. If None, this will be set to [Exception], + which will raise an error immediately if the first ToolkitWrapper in the registry fails. To try each + ToolkitWrapper that provides a suitably-named method, set this to the empty list ([]). If all + ToolkitWrappers run without raising any exceptions in this list, a single ValueError will be raised + containing the each ToolkitWrapper that was tried and the exception it raised. Raises ------ NotImplementedError if the requested method cannot be found among the registered toolkits + ValueError if no exceptions in the raise_exception_types list were raised by ToolkitWrappers, and + all ToolkitWrappers in the ToolkitRegistry were tried. + + Other forms of exceptions are possible if raise_exception_types is specified. + These are defined by the ToolkitWrapper method being called. + Examples -------- @@ -4182,28 +4400,31 @@ def call(self, method_name, *args, **kwargs): >>> smiles = toolkit_registry.call('to_smiles', molecule) """ - # TODO: catch ValueError and compile list of methods that exist but rejected the specific parameters because they did not implement the requested methods + if raise_exception_types is None: + raise_exception_types = [Exception] - value_errors = list() + errors = list() for toolkit in self._toolkits: if hasattr(toolkit, method_name): method = getattr(toolkit, method_name) try: return method(*args, **kwargs) - except NotImplementedError: - pass - except ValueError as value_error: - value_errors.append((toolkit, value_error)) + except Exception as e: + for exception_type in raise_exception_types: + if isinstance(e, exception_type): + raise e + errors.append((toolkit, e)) # No toolkit was found to provide the requested capability # TODO: Can we help developers by providing a check for typos in expected method names? - msg = 'No registered toolkits can provide the capability "{}".\n'.format( - method_name) + msg = f'No registered toolkits can provide the capability "{method_name}" ' \ + f'for args "{args}" and kwargs "{kwargs}"\n' + msg += 'Available toolkits are: {}\n'.format(self.registered_toolkits) # Append information about toolkits that implemented the method, but could not handle the provided parameters - for toolkit, value_error in value_errors: - msg += ' {} : {}\n'.format(toolkit, value_error) - raise NotImplementedError(msg) + for toolkit, error in errors: + msg += ' {} {} : {}\n'.format(toolkit, type(error), error) + raise ValueError(msg) #============================================================================================= diff --git a/openforcefield/utils/utils.py b/openforcefield/utils/utils.py index 715c8d760..0cea73c16 100644 --- a/openforcefield/utils/utils.py +++ b/openforcefield/utils/utils.py @@ -9,6 +9,7 @@ 'IncompatibleUnitError', 'inherit_docstrings', 'all_subclasses', + 'temporary_cd', 'get_data_file_path', 'unit_to_string', 'quantity_to_string', @@ -35,6 +36,7 @@ from simtk import unit import logging +import contextlib #============================================================================================= @@ -86,6 +88,26 @@ def all_subclasses(cls): g for s in cls.__subclasses__() for g in all_subclasses(s) ] +@contextlib.contextmanager +def temporary_cd(dir_path): + """Context to temporary change the working directory. + Parameters + ---------- + dir_path : str + The directory path to enter within the context + Examples + -------- + >>> dir_path = '/tmp' + >>> with temporary_cd(dir_path): + ... pass # do something in dir_path + """ + import os + prev_dir = os.getcwd() + os.chdir(os.path.abspath(dir_path)) + try: + yield + finally: + os.chdir(prev_dir) def get_data_file_path(relative_path): """Get the full path to one of the reference files in testsystems. @@ -383,17 +405,25 @@ def check_units_are_compatible(object_name, object, unit_to_check, context=None) IncompatibleUnitError """ from simtk import unit + + # If context is not provided, explicitly make it a blank string + if context is None: + context = '' + # Otherwise add a space after the end of it to correct message printing + else: + context += ' ' + if isinstance(object, list): for sub_object in object: check_units_are_compatible(object_name, sub_object, unit_to_check, context=context) elif isinstance(object, unit.Quantity): if not object.unit.is_compatible(unit_to_check): - msg = f"{context} {object_name} with " \ - f"value {object}, is incompatible with expected unit {unit_to_check}" + msg = f"{context}{object_name} with " \ + f"value {object} is incompatible with expected unit {unit_to_check}" raise IncompatibleUnitError(msg) else: - msg = f"{context} {object_name} with " \ - f"value {object}, is incompatible with expected unit {unit_to_check}" + msg = f"{context}{object_name} with " \ + f"value {object} is incompatible with expected unit {unit_to_check}" raise IncompatibleUnitError(msg)