diff --git a/README.md b/README.md index 5634f77b..039e3284 100644 --- a/README.md +++ b/README.md @@ -46,7 +46,8 @@ The CVs implemented in CVPack are listed in the table below. | [Helix RMSD content] | alpha-helix RMSD content of a sequence of residues | | [Helix torsion content] | alpha-helix Ramachandran content of a sequence of residues | | [Number of contacts] | number of contacts between two groups of atoms | -| [PathInCVSpace] | progress along (or deviation from) a predefined path in CV space | +| [OpenMM Force wrapper] | converts an OpenMM Force object into a CVPack CV | +| [Path in CV space] | progress along (or deviation from) a predefined path in CV space | | [Radius of gyration] | radius of gyration of a group of atoms | | [(Radius of gyration)^2]| square of the radius of gyration of a group of atoms | | [Residue coordination] | number of contacts between two disjoint groups of residues | @@ -109,7 +110,8 @@ Initial project based on the [CMS Cookiecutter] version 1.1. [Helix RMSD content]: https://redesignscience.github.io/cvpack/latest/api/HelixRMSDContent.html [Helix torsion content]: https://redesignscience.github.io/cvpack/latest/api/HelixTorsionContent.html [Number of contacts]: https://redesignscience.github.io/cvpack/latest/api/NumberOfContacts.html -[PathInCVSpace]: https://redesignscience.github.io/cvpack/latest/api/PathInCVSpace.html +[OpenMM Force wrapper]: https://redesignscience.github.io/cvpack/latest/api/OpenMMForceWrapper.html +[Path in CV space]: https://redesignscience.github.io/cvpack/latest/api/PathInCVSpace.html [Radius of gyration]: https://redesignscience.github.io/cvpack/latest/api/RadiusOfGyration.html [(Radius of gyration)^2]: https://redesignscience.github.io/cvpack/latest/api/RadiusOfGyrationSq.html [Residue coordination]: https://redesignscience.github.io/cvpack/latest/api/ResidueCoordination.html diff --git a/cvpack/__init__.py b/cvpack/__init__.py index 10d41396..7ac0ee68 100644 --- a/cvpack/__init__.py +++ b/cvpack/__init__.py @@ -15,6 +15,7 @@ from .helix_rmsd_content import HelixRMSDContent # noqa: F401 from .helix_torsion_content import HelixTorsionContent # noqa: F401 from .number_of_contacts import NumberOfContacts # noqa: F401 +from .openmm_force_wrapper import OpenMMForceWrapper # noqa: F401 from .path_in_cv_space import PathInCVSpace # noqa: F401 from .radius_of_gyration import RadiusOfGyration # noqa: F401 from .radius_of_gyration_sq import RadiusOfGyrationSq # noqa: F401 @@ -36,6 +37,7 @@ HelixRMSDContent, HelixTorsionContent, NumberOfContacts, + OpenMMForceWrapper, PathInCVSpace, RadiusOfGyration, RadiusOfGyrationSq, diff --git a/cvpack/angle.py b/cvpack/angle.py index 3929178e..45d81a91 100644 --- a/cvpack/angle.py +++ b/cvpack/angle.py @@ -44,15 +44,15 @@ class Angle(openmm.CustomAngleForce, BaseCollectiveVariable): Example: >>> import cvpack >>> import openmm - >>> system =openmm.System() + >>> system = openmm.System() >>> [system.addParticle(1) for i in range(3)] [0, 1, 2] >>> angle = cvpack.Angle(0, 1, 2) >>> system.addForce(angle) 0 - >>> integrator =openmm.VerletIntegrator(0) - >>> platform =openmm.Platform.getPlatformByName('Reference') - >>> context =openmm.Context(system, integrator, platform) + >>> integrator = openmm.VerletIntegrator(0) + >>> platform = openmm.Platform.getPlatformByName('Reference') + >>> context = openmm.Context(system, integrator, platform) >>> positions = [[0, 0, 0], [1, 0, 0], [1, 1, 0]] >>> context.setPositions([openmm.Vec3(*pos) for pos in positions]) >>> print(angle.getValue(context, digits=6)) diff --git a/cvpack/distance.py b/cvpack/distance.py index f8b8aba0..881d00bd 100644 --- a/cvpack/distance.py +++ b/cvpack/distance.py @@ -34,15 +34,15 @@ class Distance(openmm.CustomBondForce, BaseCollectiveVariable): Example: >>> import cvpack >>> import openmm - >>> system =openmm.System() + >>> system = openmm.System() >>> [system.addParticle(1) for i in range(2)] [0, 1] >>> distance = cvpack.Distance(0, 1) >>> system.addForce(distance) 0 - >>> integrator =openmm.VerletIntegrator(0) - >>> platform =openmm.Platform.getPlatformByName('Reference') - >>> context =openmm.Context(system, integrator, platform) + >>> integrator = openmm.VerletIntegrator(0) + >>> platform = openmm.Platform.getPlatformByName('Reference') + >>> context = openmm.Context(system, integrator, platform) >>> context.setPositions([openmm.Vec3(0, 0, 0),openmm.Vec3(1, 1, 1)]) >>> print(distance.getValue(context, digits=5)) 1.73205 nm diff --git a/cvpack/openmm_force_wrapper.py b/cvpack/openmm_force_wrapper.py new file mode 100644 index 00000000..c15967b6 --- /dev/null +++ b/cvpack/openmm_force_wrapper.py @@ -0,0 +1,81 @@ +""" +.. class:: Generic + :platform: Linux, MacOS, Windows + :synopsis: Generic collective variable + +.. classauthor:: Charlles Abreu + +""" + +import typing as t + +import openmm + +from cvpack import unit as mmunit + +from .cvpack import BaseCollectiveVariable + + +class OpenMMForceWrapper(BaseCollectiveVariable): + r""" + A collective variable whose numerical value is computed from the potential energy, + in kJ/mol, of an OpenMM force object. + + Parameters + ---------- + openmmForce + The OpenMM force whose potential energy will be used to define the collective + variable. It can be passed as an instance of :OpenMM:`Force` or as a string + containing the XML serialization of the force. + unit + The unit of measurement of the collective variable. It must be compatible + with the MD unit system (mass in `daltons`, distance in `nanometers`, time + in `picoseconds`, temperature in `kelvin`, energy in `kilojoules_per_mol`, + angle in `radians`). If the collective variables does not have a unit, use + `dimensionless`. + period + The period of the collective variable if it is periodic, or `None` if it is not. + + Example: + >>> import cvpack + >>> import numpy as np + >>> import openmm + >>> from openmm import unit + >>> from openmmtools import testsystems + >>> model = testsystems.AlanineDipeptideVacuum() + >>> angle = openmm.CustomAngleForce("theta") + >>> angle.addAngle(0, 1, 2) + 0 + >>> cv = cvpack.OpenMMForceWrapper(angle, unit.radian, period=2*np.pi) + >>> assert isinstance(cv, openmm.CustomAngleForce) + >>> cv.setUnusedForceGroup(0, model.system) + 1 + >>> model.system.addForce(cv) + 5 + >>> integrator = openmm.VerletIntegrator(0) + >>> platform = openmm.Platform.getPlatformByName("Reference") + >>> context = openmm.Context(model.system, integrator, platform) + >>> context.setPositions(model.positions) + >>> print(cv.getValue(context)) + 1.911... rad + >>> print(cv.getEffectiveMass(context)) + 0.00538... nm**2 Da/(rad**2) + """ + + yaml_tag = "!cvpack.OpenMMForceWrapper" + + def __init__( # pylint: disable=too-many-arguments, super-init-not-called + self, + openmmForce: t.Union[openmm.Force, str], + unit: mmunit.Unit, + period: t.Optional[mmunit.ScalarQuantity] = None, + ) -> None: + if isinstance(openmmForce, openmm.Force): + openmmForce = openmm.XmlSerializer.serialize(openmmForce) + unit = mmunit.SerializableUnit(unit) + force_copy = openmm.XmlSerializer.deserialize(openmmForce) + self.this = force_copy.this + self.__class__.__bases__ = (BaseCollectiveVariable, type(force_copy)) + self._registerCV(unit, openmmForce, unit, period) + if period is not None: + self._registerPeriod(period) diff --git a/cvpack/radius_of_gyration.py b/cvpack/radius_of_gyration.py index 3bc556c3..58e6ad92 100644 --- a/cvpack/radius_of_gyration.py +++ b/cvpack/radius_of_gyration.py @@ -66,9 +66,9 @@ class RadiusOfGyration(BaseRadiusOfGyration): 1 >>> model.system.addForce(radius_of_gyration) 5 - >>> platform =openmm.Platform.getPlatformByName('Reference') - >>> integrator =openmm.VerletIntegrator(0) - >>> context =openmm.Context(model.system, integrator, platform) + >>> platform = openmm.Platform.getPlatformByName('Reference') + >>> integrator = openmm.VerletIntegrator(0) + >>> context = openmm.Context(model.system, integrator, platform) >>> context.setPositions(model.positions) >>> print(radius_of_gyration.getValue(context, digits=6)) 0.2951431 nm diff --git a/cvpack/radius_of_gyration_sq.py b/cvpack/radius_of_gyration_sq.py index bab78d47..1ceabf02 100644 --- a/cvpack/radius_of_gyration_sq.py +++ b/cvpack/radius_of_gyration_sq.py @@ -68,9 +68,9 @@ class RadiusOfGyrationSq(BaseRadiusOfGyration): 1 >>> model.system.addForce(rgsq) 5 - >>> platform =openmm.Platform.getPlatformByName('Reference') - >>> integrator =openmm.VerletIntegrator(0) - >>> context =openmm.Context(model.system, integrator, platform) + >>> platform = openmm.Platform.getPlatformByName('Reference') + >>> integrator = openmm.VerletIntegrator(0) + >>> context = openmm.Context(model.system, integrator, platform) >>> context.setPositions(model.positions) >>> print(rgsq.getValue(context, digits=6)) # doctest: +ELLIPSIS 0.0871... nm**2 diff --git a/cvpack/tests/test_cvpack.py b/cvpack/tests/test_cvpack.py index d08d162a..1bf39b9c 100644 --- a/cvpack/tests/test_cvpack.py +++ b/cvpack/tests/test_cvpack.py @@ -840,3 +840,30 @@ def logsumexp(x, a=None): computed_value = -2 * sigma**2 * logsumexp(x) assert cv_value / cv_value.unit == pytest.approx(computed_value) perform_common_tests(var, context) + + +def test_openmm_force_wrapper(): + """ + Test whether an OpenMM force wrapper is computed correctly. + + """ + model = testsystems.AlanineDipeptideVacuum() + angle = openmm.CustomAngleForce("theta") + angle.addAngle(0, 1, 2) + cv = cvpack.OpenMMForceWrapper(angle, unit.radian, period=2 * np.pi) + assert isinstance(cv, openmm.CustomAngleForce) + group = cv.setUnusedForceGroup(0, model.system) + model.system.addForce(cv) + angle.setForceGroup(group + 1) + model.system.addForce(angle) + integrator = openmm.VerletIntegrator(0) + platform = openmm.Platform.getPlatformByName("Reference") + context = openmm.Context(model.system, integrator, platform) + context.setPositions(model.positions) + state = context.getState( # pylint: disable=unexpected-keyword-arg + getEnergy=True, + groups={group + 1}, + ) + angle_value = state.getPotentialEnergy().value_in_unit(unit.kilojoule_per_mole) + assert cv.getValue(context) / unit.radians == pytest.approx(angle_value) + perform_common_tests(cv, context) diff --git a/cvpack/torsion.py b/cvpack/torsion.py index 923d357b..541cbceb 100644 --- a/cvpack/torsion.py +++ b/cvpack/torsion.py @@ -51,15 +51,15 @@ class Torsion(openmm.CustomTorsionForce, BaseCollectiveVariable): Example: >>> import cvpack >>> import openmm - >>> system =openmm.System() + >>> system = openmm.System() >>> [system.addParticle(1) for i in range(4)] [0, 1, 2, 3] >>> torsion = cvpack.Torsion(0, 1, 2, 3, pbc=False) >>> system.addForce(torsion) 0 - >>> integrator =openmm.VerletIntegrator(0) - >>> platform =openmm.Platform.getPlatformByName('Reference') - >>> context =openmm.Context(system, integrator, platform) + >>> integrator = openmm.VerletIntegrator(0) + >>> platform = openmm.Platform.getPlatformByName('Reference') + >>> context = openmm.Context(system, integrator, platform) >>> positions = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1]] >>> context.setPositions([openmm.Vec3(*pos) for pos in positions]) >>> print(torsion.getValue(context, digits=6)) diff --git a/cvpack/torsion_similarity.py b/cvpack/torsion_similarity.py index 74172aaf..6a70f676 100644 --- a/cvpack/torsion_similarity.py +++ b/cvpack/torsion_similarity.py @@ -71,9 +71,9 @@ class TorsionSimilarity(openmm.CustomCompoundBondForce, BaseCollectiveVariable): 1 >>> model.system.addForce(torsion_similarity) 6 - >>> integrator =openmm.VerletIntegrator(0) - >>> platform =openmm.Platform.getPlatformByName("Reference") - >>> context =openmm.Context(model.system, integrator, platform) + >>> integrator = openmm.VerletIntegrator(0) + >>> platform = openmm.Platform.getPlatformByName("Reference") + >>> context = openmm.Context(model.system, integrator, platform) >>> context.setPositions(model.positions) >>> print(torsion_similarity.getValue(context, digits=6)) 18.65992 dimensionless