From b685637bbe17486bbb3c0759251c3d46fda76a1e Mon Sep 17 00:00:00 2001 From: Matt Thompson Date: Fri, 3 May 2024 10:07:59 -0500 Subject: [PATCH] Temporarily remove GAFF because of parmchk2 issues (#329) * Temporarily remove GAFF because of parmchk2 issues * Drop OpenEye from tests * Drop automerge action * Update Python and AmberTools versions, run tests in parallel * Remove Espaloma from test environment * Remove Espaloma CI * Refactor some tests to enable skipping GAFF * Skip SPC/E test when iterating over small molecule force fields * Update default force field in `SystemGenerator` * Streamline testing * Update testing setup * Mask CodeCov upload failures * Add release notes --- .github/workflows/ci.yaml | 37 +-- .github/workflows/dependabot_automerge.yml | 19 -- .github/workflows/espaloma_ci.yaml | 66 ----- README.md | 11 + devtools/conda-envs/test_env.yaml | 7 +- .../generators/system_generators.py | 4 +- .../generators/template_generators.py | 6 + openmmforcefields/tests/conftest.py | 23 +- .../tests/test_system_generator.py | 23 +- .../tests/test_template_generators.py | 280 +++++++++--------- 10 files changed, 213 insertions(+), 263 deletions(-) delete mode 100644 .github/workflows/dependabot_automerge.yml delete mode 100644 .github/workflows/espaloma_ci.yaml diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 5425d11f..96aa84e8 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -28,10 +28,9 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macos-latest] - python-version: ["3.9", "3.10"] # Add 3.11 in with AmberTools 23 - exclude: - - python-version: "3.10" - os: macos-latest + python-version: ["3.10", "3.11", "3.12"] + amber-conversion: [false] + charmm-conversion: [false] steps: - uses: actions/checkout@v4 @@ -43,14 +42,6 @@ jobs: create-args: >- python=${{ matrix.python-version }} - - name: License OpenEye - shell: bash -l {0} - run: | - echo "${SECRET_OE_LICENSE}" > ${OE_LICENSE} - python -c "from openeye import oechem; assert oechem.OEChemIsLicensed()" - env: - SECRET_OE_LICENSE: ${{ secrets.OE_LICENSE }} - - name: Install Package run: | pip list @@ -66,29 +57,31 @@ jobs: - name: Test Installed Package run: | - pytest -v --log-cli-level $LOGLEVEL $COV_ARGS --durations=20 \ - openmmforcefields/tests/test_amber_import.py \ - openmmforcefields/tests/test_template_generators.py \ - openmmforcefields/tests/test_system_generator.py \ - -k "not TestEspalomaTemplateGenerator" + pytest -v -n logical --log-cli-level $LOGLEVEL $COV_ARGS --durations=20 \ + openmmforcefields/tests/ + + # add --rungaff / --runespaloma to run those tests + env: COV_ARGS: --cov=openmmforcefields --cov-config=setup.cfg --cov-append --cov-report=xml LOGLEVEL: "INFO" KMP_DUPLICATE_LIB_OK: "True" - name: Test AMBER conversion + if: ${{ matrix.amber-conversion == 'true' }} run: | python convert_amber.py --input gaff.yaml --log gaff-tests.csv --verbose working-directory: ./amber - name: Test CHARMM conversion + if: ${{ matrix.charmm-conversion == 'true' }} run: | - # TODO: Uncomment these tests when new ParmEd is released + # TODO: Turn on these tests (in job matrix) when new ParmEd is released # TODO: Find a way to avoid timing out when running full charmm36.yaml conversion below - # python convert_charmm.py --verbose --in files/waters.yaml && python convert_charmm.py --verbose --in files/charmm36.yaml - # python convert_charmm.py --verbose --in files/waters.yaml - # python convert_charmm.py --verbose + python convert_charmm.py --verbose --in files/waters.yaml && python convert_charmm.py --verbose --in files/charmm36.yaml + python convert_charmm.py --verbose --in files/waters.yaml + python convert_charmm.py --verbose # TODO: Uncomment this when tests are expected to work # python test_charmm.py --verbose @@ -107,4 +100,4 @@ jobs: with: token: ${{ secrets.CODECOV_TOKEN }} file: ./coverage.xml - fail_ci_if_error: true + fail_ci_if_error: false # I don't want this to cause CI failures when a developer pushes to a fork diff --git a/.github/workflows/dependabot_automerge.yml b/.github/workflows/dependabot_automerge.yml deleted file mode 100644 index 6a41391d..00000000 --- a/.github/workflows/dependabot_automerge.yml +++ /dev/null @@ -1,19 +0,0 @@ -### .github/workflows/dependabot_automerge.yml -### As of 2300 UTC on 11 March, this workflow has secrets and a read-write token -### https://github.com/dependabot/dependabot-core/issues/3253#issuecomment-852541544 -name: Dependabot Workflow -on: - pull_request_target - -permissions: - # down scope as necessary via https://docs.github.com/en/actions/reference/authentication-in-a-workflow#modifying-the-permissions-for-the-github_token - -jobs: - do-stuff: - runs-on: ubuntu-latest - if: ${{ github.actor == 'dependabot[bot]' }} - steps: - - uses: actions/checkout@v4 - with: - ref: ${{ github.event.pull_request.head.sha }} - github-token: ${{ secrets.GITHUB_TOKEN }} diff --git a/.github/workflows/espaloma_ci.yaml b/.github/workflows/espaloma_ci.yaml deleted file mode 100644 index 8c9a83b2..00000000 --- a/.github/workflows/espaloma_ci.yaml +++ /dev/null @@ -1,66 +0,0 @@ -name: EspalomaCI - -on: - push: - branches: - - "main" - pull_request: - branches: - - "main" - schedule: - - cron: "0 0 * * *" - -defaults: - run: - shell: bash -l {0} - -concurrency: - group: "${{ github.workflow }}-${{ github.ref }}" - cancel-in-progress: true - -jobs: - test: - name: Test on ${{ matrix.os }}, Python ${{ matrix.python-version }}, Latest openff-toolkit ${{ matrix.latest-openff-toolkit }} - runs-on: ${{ matrix.os }} - env: - OE_LICENSE: ${{ github.workspace }}/oe_license.txt - strategy: - fail-fast: false - matrix: - os: [ubuntu-latest, macos-latest] - python-version: ["3.9", "3.10"] # Add 3.11 in with AmberTools 23 - exclude: - - python-version: "3.10" - os: macos-latest - - steps: - - uses: actions/checkout@v4 - - - name: Setup Conda Environment - uses: mamba-org/setup-micromamba@v1 - with: - environment-file: devtools/conda-envs/test_env.yaml - create-args: >- - python=${{ matrix.python-version }} - - - name: Install Package - run: | - pip list - micromamba list - micromamba remove --force openmmforcefields - python -m pip install . - - - name: Conda Environment Information - run: | - micromamba info - micromamba list - python -c "from openmmforcefields import __version__, __file__; print(__version__, __file__)" - - - name: Test Installed Package - run: | - pytest -v --log-cli-level $LOGLEVEL $COV_ARGS --durations=20 \ - -m "espaloma" openmmforcefields/tests --runespaloma - env: - COV_ARGS: --cov=openmmforcefields --cov-config=setup.cfg --cov-append --cov-report=xml - LOGLEVEL: "INFO" - KMP_DUPLICATE_LIB_OK: "True" diff --git a/README.md b/README.md index 651356b8..4b5879e1 100644 --- a/README.md +++ b/README.md @@ -337,6 +337,17 @@ See the corresponding directories for information on how to use the provided con # [Changelog](https://github.com/openmm/openmmforcefields/releases) +## 0.13.0 Temporarily remove GAFFTemplateGenerator + +This release temporarily removes GAFFTemplateGenerator because of packaging incompatibilities with +AmberTools 23. This functionality is planned to be re-introduced in 0.14.0. + +This release is expected to work with Python 3.10-3.12. + +Other changes include +* The default force field of `SystemGenerator` was updated from `openff-1.0.0` (code name Parsley) to + `openff-2.0.0` (code name Sage). + ## 0.12.0 Updates for espaloma and support a offxml string in SystemGenerator See our [0.12.0 release page](https://github.com/openmm/openmmforcefields/releases/tag/0.12.0) for more details. diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index 36af94ee..62feed31 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -2,7 +2,6 @@ name: openmmforcefields-test channels: - conda-forge - - openeye dependencies: # Base depends @@ -22,14 +21,12 @@ dependencies: # Requirements for conversion tools - pyyaml - - ambertools <23 + - ambertools >=22,<24 - lxml - networkx # JSON caching of residue templates - tinydb - # OpenEye toolkits are only used to speed up testing; they are not required for use - - openeye-toolkits # Validating URLs - validators @@ -37,4 +34,4 @@ dependencies: # # Espaloma requirements below # - - espaloma >=0.3.1 + # espaloma >=0.3.1 diff --git a/openmmforcefields/generators/system_generators.py b/openmmforcefields/generators/system_generators.py index 971ffb7b..41dd25cf 100644 --- a/openmmforcefields/generators/system_generators.py +++ b/openmmforcefields/generators/system_generators.py @@ -65,7 +65,7 @@ class SystemGenerator: postprocess_system : method If not None, this method will be called as ``system = postprocess_system(system)`` to post-process the System object for create_system(topology) before it is returned. """ - def __init__(self, forcefields=None, small_molecule_forcefield='openff-1.0.0', forcefield_kwargs=None, nonperiodic_forcefield_kwargs=None, periodic_forcefield_kwargs=None, template_generator_kwargs=None, barostat=None, molecules=None, cache=None, postprocess_system=None): + def __init__(self, forcefields=None, small_molecule_forcefield='openff-2.2.0', forcefield_kwargs=None, nonperiodic_forcefield_kwargs=None, periodic_forcefield_kwargs=None, template_generator_kwargs=None, barostat=None, molecules=None, cache=None, postprocess_system=None): """ This is a utility class to generate OpenMM Systems from Open Force Field Topology objects using AMBER protein force fields and GAFF small molecule force fields. @@ -206,7 +206,7 @@ def __init__(self, forcefields=None, small_molecule_forcefield='openff-1.0.0', f _logger.debug(f'Trying {template_generator_cls.__name__} to load {small_molecule_forcefield}') self.template_generator = template_generator_cls(forcefield=small_molecule_forcefield, cache=cache, template_generator_kwargs=self.template_generator_kwargs) break - except ValueError as e: + except (ValueError, NotImplementedError) as e: _logger.debug(f' {template_generator_cls.__name__} cannot load {small_molecule_forcefield}') _logger.debug(e) if self.template_generator is None: diff --git a/openmmforcefields/generators/template_generators.py b/openmmforcefields/generators/template_generators.py index 4465838c..6bd54835 100644 --- a/openmmforcefields/generators/template_generators.py +++ b/openmmforcefields/generators/template_generators.py @@ -465,6 +465,12 @@ def __init__(self, molecules=None, forcefield=None, cache=None, **kwargs): Newly parameterized molecules will be written to the cache, saving time next time! """ + raise NotImplementedError( + "This release (0.13.x) of openmmforcefields temporarily drops GAFF support and " + "thereby the GAFFTemplateGenerator class. Support will be re-introduced in " + "future releases (0.14.x). To use this class, install version 0.12.0 or older." + ) + # Initialize molecules and cache super().__init__(molecules=molecules, cache=cache) diff --git a/openmmforcefields/tests/conftest.py b/openmmforcefields/tests/conftest.py index cbbb0949..e1ae5493 100644 --- a/openmmforcefields/tests/conftest.py +++ b/openmmforcefields/tests/conftest.py @@ -4,7 +4,17 @@ def pytest_addoption(parser): parser.addoption( - "--runespaloma", action="store_true", default=False, help="run espaloma tests" + "--runespaloma", + action="store_true", + default=False, + help="run espaloma tests", + ) + + parser.addoption( + "--rungaff", + action="store_true", + default=False, + help="run gaff tests", ) @@ -18,4 +28,13 @@ def pytest_collection_modifyitems(config, items): if not config.getoption("--runespaloma"): for item in items: if "espaloma" in item.keywords: - item.add_marker(skip_slow) \ No newline at end of file + item.add_marker(skip_slow) + + del skip_slow + + skip_slow = pytest.mark.skip(reason="need --rungaff option to run") + + if not config.getoption("--rungaff"): + for item in items: + if "gaff" in item.keywords: + item.add_marker(skip_slow) diff --git a/openmmforcefields/tests/test_system_generator.py b/openmmforcefields/tests/test_system_generator.py index 537db94e..57cafa8c 100644 --- a/openmmforcefields/tests/test_system_generator.py +++ b/openmmforcefields/tests/test_system_generator.py @@ -176,7 +176,7 @@ def test_barostat(self): assert force.getFrequency() == frequency @pytest.mark.parametrize("small_molecule_forcefield", [ - 'gaff-2.11', + pytest.param('gaff-2.11', marks=pytest.mark.gaff), 'openff-2.0.0', pytest.param('espaloma-0.3.2', marks=pytest.mark.espaloma)]) def test_create_with_template_generator(self, small_molecule_forcefield): @@ -198,7 +198,7 @@ def test_create_with_template_generator(self, small_molecule_forcefield): del generator @pytest.mark.parametrize("small_molecule_forcefield", [ - 'gaff-2.11', + pytest.param('gaff-2.11', marks=pytest.mark.gaff), 'openff-2.0.0', pytest.param('espaloma-0.3.2', marks=pytest.mark.espaloma)]) def test_forcefield_default_kwargs(self, small_molecule_forcefield, test_systems): @@ -234,7 +234,7 @@ def test_forcefield_default_kwargs(self, small_molecule_forcefield, test_systems assert forces['NonbondedForce'].getNonbondedMethod() == openmm.NonbondedForce.PME, "Expected LJPME, got {forces['NonbondedForce'].getNonbondedMethod()}" @pytest.mark.parametrize("small_molecule_forcefield", [ - 'gaff-2.11', + pytest.param('gaff-2.11', marks=pytest.mark.gaff), 'openff-2.0.0', pytest.param('espaloma-0.3.2', marks=pytest.mark.espaloma)]) def test_forcefield_kwargs(self, small_molecule_forcefield, test_systems): @@ -280,7 +280,7 @@ def test_forcefield_kwargs(self, small_molecule_forcefield, test_systems): 'NonbondedForce'].getNonbondedMethod() == openmm.NonbondedForce.LJPME, "Expected LJPME, got {forces['NonbondedForce'].getNonbondedMethod()}" @pytest.mark.parametrize("small_molecule_forcefield", [ - 'gaff-2.11', + pytest.param('gaff-2.11', marks=pytest.mark.gaff), 'openff-2.0.0', pytest.param('espaloma-0.3.2', marks=pytest.mark.espaloma)]) def test_parameterize_molecules_from_creation(self, test_systems, small_molecule_forcefield): @@ -307,7 +307,7 @@ def test_parameterize_molecules_from_creation(self, test_systems, small_molecule assert (t2.interval() < t1.interval()) @pytest.mark.parametrize("small_molecule_forcefield", [ - 'gaff-2.11', + pytest.param('gaff-2.11', marks=pytest.mark.gaff), 'openff-2.0.0', pytest.param('espaloma-0.3.2', marks=pytest.mark.espaloma)]) def test_parameterize_molecules_specified_during_create_system(self, test_systems, small_molecule_forcefield): @@ -325,10 +325,13 @@ def test_parameterize_molecules_specified_during_create_system(self, test_system # Specify molecules during system creation system = generator.create_system(openmm_topology, molecules=molecules) - @pytest.mark.parametrize("small_molecule_forcefield", [ - 'gaff-2.11', - 'openff-2.0.0', - pytest.param('espaloma-0.3.2', marks=pytest.mark.espaloma)]) + @pytest.mark.parametrize( + "small_molecule_forcefield", [ + pytest.param('gaff-2.11', marks=pytest.mark.gaff), + 'openff-2.0.0', + pytest.param('espaloma-0.3.2', marks=pytest.mark.espaloma), + ] + ) def test_add_molecules(self, test_systems, small_molecule_forcefield): """Test that Molecules can be added to SystemGenerator later""" # Create a SystemGenerator for this force field @@ -355,7 +358,7 @@ def test_add_molecules(self, test_systems, small_molecule_forcefield): assert (t2.interval() < t1.interval()) @pytest.mark.parametrize("small_molecule_forcefield", [ - 'gaff-2.11', + pytest.param('gaff-2.11', marks=pytest.mark.gaff), 'openff-2.0.0', pytest.param('espaloma-0.3.2', marks=pytest.mark.espaloma)]) def test_cache(self, test_systems, small_molecule_forcefield): diff --git a/openmmforcefields/tests/test_template_generators.py b/openmmforcefields/tests/test_template_generators.py index 6c43e3d3..d117fbe6 100644 --- a/openmmforcefields/tests/test_template_generators.py +++ b/openmmforcefields/tests/test_template_generators.py @@ -28,10 +28,7 @@ # Tests ################################################################################ -class TestGAFFTemplateGenerator(unittest.TestCase): - TEMPLATE_GENERATOR = GAFFTemplateGenerator - - amber_forcefields = ['amber/protein.ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml'] +class TemplateGeneratorBaseCase(unittest.TestCase): def filter_molecules(self, molecules): """ @@ -89,6 +86,144 @@ def setUp(self): for name in ['parmed', 'matplotlib']: logging.getLogger(name).setLevel(logging.WARNING) + @classmethod + def compare_energies(cls, molecule, template_generated_system, reference_system): + """Compare energies between OpenMM System generated by reference method and OpenMM System generated by ForceField template. + + The OpenMM System object created internally by the reference method is used to + avoid any issues with stochasticity of partial charges due to conformer generation. + + Parameters + ---------- + molecule : openff.toolkit.topology.Molecule + The Molecule object to compare energy components (including positions) + template_generated_system : openmm.System + System generated by OpenMM ForceField template + reference_system : openmm.System + System generated by reference parmaeterization engine + + """ + from openff.units.openmm import ensure_quantity + + # Compute energies + reference_energy, reference_forces = cls.compute_energy( + template_generated_system, + ensure_quantity(molecule.conformers[0], "openmm"), + ) + template_energy, template_forces = cls.compute_energy( + reference_system, + ensure_quantity(molecule.conformers[0], "openmm"), + ) + + from openmm import unit + + def write_xml(filename, system): + with open(filename, 'w') as outfile: + print(f'Writing {filename}...') + outfile.write(openmm.XmlSerializer.serialize(system)) + # DEBUG + print(openmm.XmlSerializer.serialize(system)) + + # Make sure both systems contain the same energy components + reference_components = set(reference_energy['components']) + template_components = set(template_energy['components']) + if len(reference_components.difference(template_components)) > 0: + raise Exception(f'Reference system contains components {reference_components.difference(template_components)} that do not appear in template-generated system.') + if len(template_components.difference(reference_components)) > 0: + raise Exception(f'Template-generated system contains components {template_components.difference(reference_components)} that do not appear in reference system.') + components = reference_components + + # Compare energies + ENERGY_DEVIATION_TOLERANCE = 1.0e-2 * unit.kilocalories_per_mole + delta = (template_energy['total'] - reference_energy['total']) + if abs(delta) > ENERGY_DEVIATION_TOLERANCE: + # Show breakdown by components + print('Energy components:') + print(f"{'component':24} {'Template (kcal/mol)':>20} {'Reference (kcal/mol)':>20}") + for key in components: + reference_component_energy = reference_energy['components'][key] + template_component_energy = template_energy['components'][key] + print(f'{key:24} {(template_component_energy/unit.kilocalories_per_mole):20.3f} {(reference_component_energy/unit.kilocalories_per_mole):20.3f} kcal/mol') + print(f'{"TOTAL":24} {(template_energy["total"]/unit.kilocalories_per_mole):20.3f} {(reference_energy["total"]/unit.kilocalories_per_mole):20.3f} kcal/mol') + write_xml('reference_system.xml', reference_system) + write_xml('template_system.xml', template_system) # What's this? This variable does not exist + raise Exception(f'Energy deviation for {molecule.to_smiles()} ({delta/unit.kilocalories_per_mole} kcal/mol) exceeds threshold ({ENERGY_DEVIATION_TOLERANCE})') + + # Compare forces + def norm(x): + N = x.shape[0] + return np.sqrt((1.0/N) * (x**2).sum()) + def relative_deviation(x, y): + FORCE_UNIT = unit.kilocalories_per_mole / unit.angstroms + if hasattr(x, 'value_in_unit'): + x = x / FORCE_UNIT + if hasattr(y, 'value_in_unit'): + y = y / FORCE_UNIT + + if norm(y) > 0: + return norm(x-y) / np.sqrt(norm(x)**2 + norm(y)**2) + else: + return 0 + + RELATIVE_FORCE_DEVIATION_TOLERANCE = 1.0e-5 + relative_force_deviation = relative_deviation(template_forces['total'], reference_forces['total']) + if relative_force_deviation > RELATIVE_FORCE_DEVIATION_TOLERANCE: + # Show breakdown by components + print('Force components:') + print(f"{'component':24} {'relative deviation':>24}") + for key in components: + print(f"{key:24} {relative_deviation(template_forces['components'][key], reference_forces['components'][key]):24.10f}") + print(f'{"TOTAL":24} {relative_force_deviation:24.10f}') + write_xml('system-smirnoff.xml', reference_system) + write_xml('openmm-smirnoff.xml', template_generated_system) + raise Exception(f'Relative force deviation for {molecule.to_smiles()} ({relative_force_deviation}) exceeds threshold ({RELATIVE_FORCE_DEVIATION_TOLERANCE})') + + @staticmethod + def compute_energy(system, positions): + """Compute potential energy and Force components for an OpenMM system. + + Parameters + ---------- + system : openmm.System + The System object + positions : openmm.unit.Quantity of shape (nparticles,3) with units compatible with nanometers + The positions for which energy is to be computed + + Returns + ------- + openmm_energy : dict of str : openmm.unit.Quantity + openmm_energy['total'] is the total potential energy + openmm_energy['components'][forcename] is the potential energy for the specified component force + openmm_forces : dict of str : openmm.unit.Quantity + openmm_forces['total'] is the total force + openmm_forces['components'][forcename] is the force for the specified component force + """ + system = copy.deepcopy(system) + for index, force in enumerate(system.getForces()): + force.setForceGroup(index) + platform = openmm.Platform.getPlatformByName('Reference') + integrator = openmm.VerletIntegrator(0.001) + context = openmm.Context(system, integrator, platform) + context.setPositions(positions) + openmm_energy = { + 'total' : context.getState(getEnergy=True).getPotentialEnergy(), + 'components' : { system.getForce(index).__class__.__name__ : context.getState(getEnergy=True, groups=(1 << index)).getPotentialEnergy() for index in range(system.getNumForces()) }, + } + + openmm_forces = { + 'total' : context.getState(getForces=True).getForces(asNumpy=True), + 'components' : { system.getForce(index).__class__.__name__ : context.getState(getForces=True, groups=(1 << index)).getForces(asNumpy=True) for index in range(system.getNumForces()) }, + } + + del context, integrator + return openmm_energy, openmm_forces + +@pytest.mark.skip(reason='Skip GAFF tests') +class TestGAFFTemplateGenerator(TemplateGeneratorBaseCase): + TEMPLATE_GENERATOR = GAFFTemplateGenerator + + amber_forcefields = ['amber/protein.ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml'] + def test_version(self): """Test version""" for forcefield in GAFFTemplateGenerator.INSTALLED_FORCEFIELDS: @@ -584,139 +719,8 @@ def test_multiple_registration(self): system = forcefield.createSystem(openmm_topology, nonbondedMethod=NoCutoff) assert system.getNumParticles() == molecule.n_atoms - @staticmethod - def compute_energy(system, positions): - """Compute potential energy and Force components for an OpenMM system. - - Parameters - ---------- - system : openmm.System - The System object - positions : openmm.unit.Quantity of shape (nparticles,3) with units compatible with nanometers - The positions for which energy is to be computed - - Returns - ------- - openmm_energy : dict of str : openmm.unit.Quantity - openmm_energy['total'] is the total potential energy - openmm_energy['components'][forcename] is the potential energy for the specified component force - openmm_forces : dict of str : openmm.unit.Quantity - openmm_forces['total'] is the total force - openmm_forces['components'][forcename] is the force for the specified component force - """ - system = copy.deepcopy(system) - for index, force in enumerate(system.getForces()): - force.setForceGroup(index) - platform = openmm.Platform.getPlatformByName('Reference') - integrator = openmm.VerletIntegrator(0.001) - context = openmm.Context(system, integrator, platform) - context.setPositions(positions) - openmm_energy = { - 'total' : context.getState(getEnergy=True).getPotentialEnergy(), - 'components' : { system.getForce(index).__class__.__name__ : context.getState(getEnergy=True, groups=(1 << index)).getPotentialEnergy() for index in range(system.getNumForces()) }, - } - - openmm_forces = { - 'total' : context.getState(getForces=True).getForces(asNumpy=True), - 'components' : { system.getForce(index).__class__.__name__ : context.getState(getForces=True, groups=(1 << index)).getForces(asNumpy=True) for index in range(system.getNumForces()) }, - } - - del context, integrator - return openmm_energy, openmm_forces - - @classmethod - def compare_energies(cls, molecule, template_generated_system, reference_system): - """Compare energies between OpenMM System generated by reference method and OpenMM System generated by ForceField template. - - The OpenMM System object created internally by the reference method is used to - avoid any issues with stochasticity of partial charges due to conformer generation. - - Parameters - ---------- - molecule : openff.toolkit.topology.Molecule - The Molecule object to compare energy components (including positions) - template_generated_system : openmm.System - System generated by OpenMM ForceField template - reference_system : openmm.System - System generated by reference parmaeterization engine - - """ - from openff.units.openmm import ensure_quantity - # Compute energies - reference_energy, reference_forces = cls.compute_energy( - template_generated_system, - ensure_quantity(molecule.conformers[0], "openmm"), - ) - template_energy, template_forces = cls.compute_energy( - reference_system, - ensure_quantity(molecule.conformers[0], "openmm"), - ) - - from openmm import unit - - def write_xml(filename, system): - with open(filename, 'w') as outfile: - print(f'Writing {filename}...') - outfile.write(openmm.XmlSerializer.serialize(system)) - # DEBUG - print(openmm.XmlSerializer.serialize(system)) - - # Make sure both systems contain the same energy components - reference_components = set(reference_energy['components']) - template_components = set(template_energy['components']) - if len(reference_components.difference(template_components)) > 0: - raise Exception(f'Reference system contains components {reference_components.difference(template_components)} that do not appear in template-generated system.') - if len(template_components.difference(reference_components)) > 0: - raise Exception(f'Template-generated system contains components {template_components.difference(reference_components)} that do not appear in reference system.') - components = reference_components - - # Compare energies - ENERGY_DEVIATION_TOLERANCE = 1.0e-2 * unit.kilocalories_per_mole - delta = (template_energy['total'] - reference_energy['total']) - if abs(delta) > ENERGY_DEVIATION_TOLERANCE: - # Show breakdown by components - print('Energy components:') - print(f"{'component':24} {'Template (kcal/mol)':>20} {'Reference (kcal/mol)':>20}") - for key in components: - reference_component_energy = reference_energy['components'][key] - template_component_energy = template_energy['components'][key] - print(f'{key:24} {(template_component_energy/unit.kilocalories_per_mole):20.3f} {(reference_component_energy/unit.kilocalories_per_mole):20.3f} kcal/mol') - print(f'{"TOTAL":24} {(template_energy["total"]/unit.kilocalories_per_mole):20.3f} {(reference_energy["total"]/unit.kilocalories_per_mole):20.3f} kcal/mol') - write_xml('reference_system.xml', reference_system) - write_xml('template_system.xml', template_system) # What's this? This variable does not exist - raise Exception(f'Energy deviation for {molecule.to_smiles()} ({delta/unit.kilocalories_per_mole} kcal/mol) exceeds threshold ({ENERGY_DEVIATION_TOLERANCE})') - - # Compare forces - def norm(x): - N = x.shape[0] - return np.sqrt((1.0/N) * (x**2).sum()) - def relative_deviation(x, y): - FORCE_UNIT = unit.kilocalories_per_mole / unit.angstroms - if hasattr(x, 'value_in_unit'): - x = x / FORCE_UNIT - if hasattr(y, 'value_in_unit'): - y = y / FORCE_UNIT - - if norm(y) > 0: - return norm(x-y) / np.sqrt(norm(x)**2 + norm(y)**2) - else: - return 0 - - RELATIVE_FORCE_DEVIATION_TOLERANCE = 1.0e-5 - relative_force_deviation = relative_deviation(template_forces['total'], reference_forces['total']) - if relative_force_deviation > RELATIVE_FORCE_DEVIATION_TOLERANCE: - # Show breakdown by components - print('Force components:') - print(f"{'component':24} {'relative deviation':>24}") - for key in components: - print(f"{key:24} {relative_deviation(template_forces['components'][key], reference_forces['components'][key]):24.10f}") - print(f'{"TOTAL":24} {relative_force_deviation:24.10f}') - write_xml('system-smirnoff.xml', reference_system) - write_xml('openmm-smirnoff.xml', template_generated_system) - raise Exception(f'Relative force deviation for {molecule.to_smiles()} ({relative_force_deviation}) exceeds threshold ({RELATIVE_FORCE_DEVIATION_TOLERANCE})') - -class TestSMIRNOFFTemplateGenerator(TestGAFFTemplateGenerator): +class TestSMIRNOFFTemplateGenerator(TemplateGeneratorBaseCase): TEMPLATE_GENERATOR = SMIRNOFFTemplateGenerator def propagate_dynamics(self, molecule, system): @@ -840,6 +844,8 @@ def test_partial_charges_are_none(self): continue if "opc" in small_molecule_forcefield: continue + if "spce" in small_molecule_forcefield: + continue print(f'Testing energies for {small_molecule_forcefield}...') # Create a generator that knows about a few molecules @@ -931,7 +937,7 @@ def test_14_scaling_from_offxml(self): assert exception[4]._value == 0. @pytest.mark.espaloma -class TestEspalomaTemplateGenerator(TestGAFFTemplateGenerator): +class TestEspalomaTemplateGenerator(TemplateGeneratorBaseCase): TEMPLATE_GENERATOR = EspalomaTemplateGenerator def propagate_dynamics(self, molecule, system):