From e9dec9c5c197f6e80f2d3c30388e17b12bd3aec1 Mon Sep 17 00:00:00 2001 From: Charlles Abreu Date: Thu, 7 Mar 2024 19:12:21 -0500 Subject: [PATCH 1/6] Reformatted doctests --- cvpack/angle.py | 8 ++++---- cvpack/distance.py | 8 ++++---- cvpack/radius_of_gyration.py | 6 +++--- cvpack/radius_of_gyration_sq.py | 6 +++--- cvpack/torsion.py | 8 ++++---- cvpack/torsion_similarity.py | 6 +++--- 6 files changed, 21 insertions(+), 21 deletions(-) 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/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/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 From 91b29cc3c770f03f17fea05edef0e165c0310e58 Mon Sep 17 00:00:00 2001 From: Charlles Abreu Date: Thu, 7 Mar 2024 19:35:56 -0500 Subject: [PATCH 2/6] Implemented generic collective variable --- cvpack/__init__.py | 2 ++ cvpack/generic.py | 77 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+) create mode 100644 cvpack/generic.py diff --git a/cvpack/__init__.py b/cvpack/__init__.py index 10d41396..7e507321 100644 --- a/cvpack/__init__.py +++ b/cvpack/__init__.py @@ -10,6 +10,7 @@ from .centroid_function import CentroidFunction # noqa: F401 from .composite_rmsd import CompositeRMSD # noqa: F401 from .distance import Distance # noqa: F401 +from .generic import Generic # noqa: F401 from .helix_angle_content import HelixAngleContent # noqa: F401 from .helix_hbond_content import HelixHBondContent # noqa: F401 from .helix_rmsd_content import HelixRMSDContent # noqa: F401 @@ -31,6 +32,7 @@ CentroidFunction, CompositeRMSD, Distance, + Generic, HelixAngleContent, HelixHBondContent, HelixRMSDContent, diff --git a/cvpack/generic.py b/cvpack/generic.py new file mode 100644 index 00000000..0b7a0ff1 --- /dev/null +++ b/cvpack/generic.py @@ -0,0 +1,77 @@ +""" +.. 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 Generic(BaseCollectiveVariable): + r""" + A generic collective variable. + + Parameters + ---------- + atom1 + The index of the first atom + atom2 + The index of the second atom + atom3 + The index of the third atom + pbc + Whether to use periodic boundary conditions + + Example: + >>> import cvpack + >>> import numpy as np + >>> import openmm + >>> from openmm import unit + >>> from openmmtools import testsystems + >>> model = testsystems.AlanineDipeptideVacuum() + >>> cv1 = cvpack.Angle(0, 1, 2) + >>> model.system.addForce(cv1) + 5 + >>> angle = openmm.CustomAngleForce("theta") + >>> angle.addAngle(0, 1, 2) + 0 + >>> cv2 = cvpack.Generic(angle, unit.radian, period=2*np.pi) + >>> model.system.addForce(cv2) + 6 + >>> integrator = openmm.VerletIntegrator(0) + >>> platform = openmm.Platform.getPlatformByName('Reference') + >>> context = openmm.Context(model.system, integrator, platform) + >>> context.setPositions(model.positions) + >>> print(cv1.getValue(context)) + 1.911... rad + >>> print(cv2.getValue(context)) + 1.911... rad + >>> assert isinstance(cv2, openmm.CustomAngleForce) + """ + + yaml_tag = "!cvpack.Generic" + + 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.__dict__.update(force_copy.__dict__) + self.__class__.__bases__ += (force_copy.__class__,) + self._registerCV(unit, openmmForce, unit, period) + if period is not None: + self._registerPeriod(period) From 1601bd28fbde13e78c6d9e4cb27c9421d0e933b9 Mon Sep 17 00:00:00 2001 From: Charlles Abreu Date: Thu, 7 Mar 2024 20:13:01 -0500 Subject: [PATCH 3/6] Changed name to OpenMMForceWrapper --- cvpack/__init__.py | 4 +- cvpack/generic.py | 77 ---------------------------- cvpack/openmm_force_wrapper.py | 92 ++++++++++++++++++++++++++++++++++ 3 files changed, 94 insertions(+), 79 deletions(-) delete mode 100644 cvpack/generic.py create mode 100644 cvpack/openmm_force_wrapper.py diff --git a/cvpack/__init__.py b/cvpack/__init__.py index 7e507321..7ac0ee68 100644 --- a/cvpack/__init__.py +++ b/cvpack/__init__.py @@ -10,12 +10,12 @@ from .centroid_function import CentroidFunction # noqa: F401 from .composite_rmsd import CompositeRMSD # noqa: F401 from .distance import Distance # noqa: F401 -from .generic import Generic # noqa: F401 from .helix_angle_content import HelixAngleContent # noqa: F401 from .helix_hbond_content import HelixHBondContent # noqa: F401 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 @@ -32,12 +32,12 @@ CentroidFunction, CompositeRMSD, Distance, - Generic, HelixAngleContent, HelixHBondContent, HelixRMSDContent, HelixTorsionContent, NumberOfContacts, + OpenMMForceWrapper, PathInCVSpace, RadiusOfGyration, RadiusOfGyrationSq, diff --git a/cvpack/generic.py b/cvpack/generic.py deleted file mode 100644 index 0b7a0ff1..00000000 --- a/cvpack/generic.py +++ /dev/null @@ -1,77 +0,0 @@ -""" -.. 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 Generic(BaseCollectiveVariable): - r""" - A generic collective variable. - - Parameters - ---------- - atom1 - The index of the first atom - atom2 - The index of the second atom - atom3 - The index of the third atom - pbc - Whether to use periodic boundary conditions - - Example: - >>> import cvpack - >>> import numpy as np - >>> import openmm - >>> from openmm import unit - >>> from openmmtools import testsystems - >>> model = testsystems.AlanineDipeptideVacuum() - >>> cv1 = cvpack.Angle(0, 1, 2) - >>> model.system.addForce(cv1) - 5 - >>> angle = openmm.CustomAngleForce("theta") - >>> angle.addAngle(0, 1, 2) - 0 - >>> cv2 = cvpack.Generic(angle, unit.radian, period=2*np.pi) - >>> model.system.addForce(cv2) - 6 - >>> integrator = openmm.VerletIntegrator(0) - >>> platform = openmm.Platform.getPlatformByName('Reference') - >>> context = openmm.Context(model.system, integrator, platform) - >>> context.setPositions(model.positions) - >>> print(cv1.getValue(context)) - 1.911... rad - >>> print(cv2.getValue(context)) - 1.911... rad - >>> assert isinstance(cv2, openmm.CustomAngleForce) - """ - - yaml_tag = "!cvpack.Generic" - - 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.__dict__.update(force_copy.__dict__) - self.__class__.__bases__ += (force_copy.__class__,) - self._registerCV(unit, openmmForce, unit, period) - if period is not None: - self._registerPeriod(period) diff --git a/cvpack/openmm_force_wrapper.py b/cvpack/openmm_force_wrapper.py new file mode 100644 index 00000000..cc0cd199 --- /dev/null +++ b/cvpack/openmm_force_wrapper.py @@ -0,0 +1,92 @@ +""" +.. 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 generic collective variable. + + Parameters + ---------- + function + The function to be evaluated. It must be a valid + :OpenMM:`CustomCentroidBondForce` expression + 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` + groups + The groups of atoms to be used in the function. Each group must be specified + as a list of atom indices with arbitrary length + collections + The indices of the groups in each collection, passed as a 2D array-like object + of shape `(m, n)`, where `m` is the number of collections and `n` is the number + groups per collection. If a 1D object is passed, it is assumed that `m` is 1 and + `n` is the length of the object. + period + The period of the collective variable if it is periodic, or `None` if it is not + pbc + Whether to use periodic boundary conditions + weighByMass + Whether to define the centroid as the center of mass of the group instead of + the geometric center + + 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.__dict__.update(force_copy.__dict__) + self.__class__.__bases__ += (force_copy.__class__,) + self._registerCV(unit, openmmForce, unit, period) + if period is not None: + self._registerPeriod(period) From 75148a46f7e29f0aee768c0317f599f61e7d7e77 Mon Sep 17 00:00:00 2001 From: Charlles Abreu Date: Fri, 15 Mar 2024 15:10:22 -0400 Subject: [PATCH 4/6] Improved base-class handling --- cvpack/openmm_force_wrapper.py | 30 +++++++++--------------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/cvpack/openmm_force_wrapper.py b/cvpack/openmm_force_wrapper.py index cc0cd199..44dd786c 100644 --- a/cvpack/openmm_force_wrapper.py +++ b/cvpack/openmm_force_wrapper.py @@ -22,30 +22,18 @@ class OpenMMForceWrapper(BaseCollectiveVariable): Parameters ---------- - function - The function to be evaluated. It must be a valid - :OpenMM:`CustomCentroidBondForce` expression + 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` - groups - The groups of atoms to be used in the function. Each group must be specified - as a list of atom indices with arbitrary length - collections - The indices of the groups in each collection, passed as a 2D array-like object - of shape `(m, n)`, where `m` is the number of collections and `n` is the number - groups per collection. If a 1D object is passed, it is assumed that `m` is 1 and - `n` is the length of the object. + `dimensionless`. period - The period of the collective variable if it is periodic, or `None` if it is not - pbc - Whether to use periodic boundary conditions - weighByMass - Whether to define the centroid as the center of mass of the group instead of - the geometric center + The period of the collective variable if it is periodic, or `None` if it is not. Example: >>> import cvpack @@ -64,7 +52,7 @@ class OpenMMForceWrapper(BaseCollectiveVariable): >>> model.system.addForce(cv) 5 >>> integrator = openmm.VerletIntegrator(0) - >>> platform = openmm.Platform.getPlatformByName('Reference') + >>> platform = openmm.Platform.getPlatformByName("Reference") >>> context = openmm.Context(model.system, integrator, platform) >>> context.setPositions(model.positions) >>> print(cv.getValue(context)) @@ -85,8 +73,8 @@ def __init__( # pylint: disable=too-many-arguments, super-init-not-called openmmForce = openmm.XmlSerializer.serialize(openmmForce) unit = mmunit.SerializableUnit(unit) force_copy = openmm.XmlSerializer.deserialize(openmmForce) - self.__dict__.update(force_copy.__dict__) - self.__class__.__bases__ += (force_copy.__class__,) + 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) From ce4dadd683214359c0aad86565296d49444fe5a0 Mon Sep 17 00:00:00 2001 From: Charlles Abreu Date: Fri, 15 Mar 2024 15:10:30 -0400 Subject: [PATCH 5/6] Added unit test --- cvpack/tests/test_cvpack.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) 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) From 3787161547077d434e185dc6fa35cf59a1f3dcab Mon Sep 17 00:00:00 2001 From: Charlles Abreu Date: Fri, 15 Mar 2024 15:58:42 -0400 Subject: [PATCH 6/6] Updated docstring and README.md --- README.md | 1 + cvpack/openmm_force_wrapper.py | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index af918e2c..8768bc1b 100644 --- a/README.md +++ b/README.md @@ -46,6 +46,7 @@ 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 | +| [OpenMM Force wrapper] | converts an OpenMM Force object into a CVPack CV | | [PathInCVSpace] | progress along (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 | diff --git a/cvpack/openmm_force_wrapper.py b/cvpack/openmm_force_wrapper.py index 44dd786c..c15967b6 100644 --- a/cvpack/openmm_force_wrapper.py +++ b/cvpack/openmm_force_wrapper.py @@ -18,7 +18,8 @@ class OpenMMForceWrapper(BaseCollectiveVariable): r""" - A generic collective variable. + A collective variable whose numerical value is computed from the potential energy, + in kJ/mol, of an OpenMM force object. Parameters ----------