Unverified Commit 54fa65ed authored by Peter Eastman's avatar Peter Eastman Committed by GitHub
Browse files

Created OrientationRestraintForce (#5048)

* Created OrientationRestraintForce

* Python API

* Documentation

* Serialization

* GPU implementation of OrientationRestraintForce

* Bug fix

* Remove unused code
parent 44524bd1
...@@ -1175,6 +1175,42 @@ private: ...@@ -1175,6 +1175,42 @@ private:
std::vector<int> particles; std::vector<int> particles;
}; };
/**
* This kernel is invoked by OrientationRestraintForce to calculate the forces acting on the system and the energy of the system.
*/
class ReferenceCalcOrientationRestraintForceKernel : public CalcOrientationRestraintForceKernel {
public:
ReferenceCalcOrientationRestraintForceKernel(std::string name, const Platform& platform) : CalcOrientationRestraintForceKernel(name, platform) {
}
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the OrientationRestraintForce this kernel will be used for
*/
void initialize(const System& system, const OrientationRestraintForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the OrientationRestraintForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const OrientationRestraintForce& force);
private:
double k;
std::vector<Vec3> referencePos;
std::vector<int> particles;
};
/** /**
* This kernel is invoked by VerletIntegrator to take one time step. * This kernel is invoked by VerletIntegrator to take one time step.
*/ */
......
/* Portions copyright (c) 2025 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __ReferenceOrientationRestraintForce_H__
#define __ReferenceOrientationRestraintForce_H__
#include "openmm/OrientationRestraintForce.h"
#include <vector>
namespace OpenMM {
class ReferenceOrientationRestraintForce {
private:
double k;
std::vector<OpenMM::Vec3> referencePos;
std::vector<int> particles;
public:
/**
* Constructor
*/
ReferenceOrientationRestraintForce(double k, std::vector<OpenMM::Vec3>& referencePos, std::vector<int>& particles);
/**
* Destructor
*/
~ReferenceOrientationRestraintForce();
/**
* Calculate the interaction.
*
* @param atomCoordinates atom coordinates
* @param forces the forces are added to this
* @return the energy of the interaction
*/
double calculateIxn(std::vector<OpenMM::Vec3>& atomCoordinates, std::vector<OpenMM::Vec3>& forces) const;
};
} // namespace OpenMM
#endif // __ReferenceOrientationRestraintForce_H__
...@@ -84,6 +84,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla ...@@ -84,6 +84,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
return new ReferenceCalcATMForceKernel(name, platform); return new ReferenceCalcATMForceKernel(name, platform);
if (name == CalcCustomCPPForceKernel::Name()) if (name == CalcCustomCPPForceKernel::Name())
return new ReferenceCalcCustomCPPForceKernel(name, platform); return new ReferenceCalcCustomCPPForceKernel(name, platform);
if (name == CalcOrientationRestraintForceKernel::Name())
return new ReferenceCalcOrientationRestraintForceKernel(name, platform);
if (name == CalcRGForceKernel::Name()) if (name == CalcRGForceKernel::Name())
return new ReferenceCalcRGForceKernel(name, platform); return new ReferenceCalcRGForceKernel(name, platform);
if (name == CalcRMSDForceKernel::Name()) if (name == CalcRMSDForceKernel::Name())
......
...@@ -58,6 +58,7 @@ ...@@ -58,6 +58,7 @@
#include "ReferenceMonteCarloBarostat.h" #include "ReferenceMonteCarloBarostat.h"
#include "ReferenceNoseHooverChain.h" #include "ReferenceNoseHooverChain.h"
#include "ReferenceNoseHooverDynamics.h" #include "ReferenceNoseHooverDynamics.h"
#include "ReferenceOrientationRestraintForce.h"
#include "ReferencePointFunctions.h" #include "ReferencePointFunctions.h"
#include "ReferenceProperDihedralBond.h" #include "ReferenceProperDihedralBond.h"
#include "ReferenceQTBDynamics.h" #include "ReferenceQTBDynamics.h"
...@@ -2323,6 +2324,45 @@ double ReferenceCalcRGForceKernel::execute(ContextImpl& context, bool includeFor ...@@ -2323,6 +2324,45 @@ double ReferenceCalcRGForceKernel::execute(ContextImpl& context, bool includeFor
return rg.calculateIxn(posData, forceData); return rg.calculateIxn(posData, forceData);
} }
void ReferenceCalcOrientationRestraintForceKernel::initialize(const System& system, const OrientationRestraintForce& force) {
k = force.getK();
particles = force.getParticles();
if (particles.size() == 0)
for (int i = 0; i < system.getNumParticles(); i++)
particles.push_back(i);
referencePos = force.getReferencePositions();
Vec3 center;
for (int i : particles)
center += referencePos[i];
center /= particles.size();
for (Vec3& p : referencePos)
p -= center;
}
double ReferenceCalcOrientationRestraintForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
ReferenceOrientationRestraintForce f(k, referencePos, particles);
return f.calculateIxn(posData, forceData);
}
void ReferenceCalcOrientationRestraintForceKernel::copyParametersToContext(ContextImpl& context, const OrientationRestraintForce& force) {
if (referencePos.size() != force.getReferencePositions().size())
throw OpenMMException("updateParametersInContext: The number of reference positions has changed");
k = force.getK();
particles = force.getParticles();
if (particles.size() == 0)
for (int i = 0; i < referencePos.size(); i++)
particles.push_back(i);
referencePos = force.getReferencePositions();
Vec3 center;
for (int i : particles)
center += referencePos[i];
center /= particles.size();
for (Vec3& p : referencePos)
p -= center;
}
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() { ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
if (dynamics) if (dynamics)
delete dynamics; delete dynamics;
......
...@@ -65,6 +65,7 @@ ReferencePlatform::ReferencePlatform() { ...@@ -65,6 +65,7 @@ ReferencePlatform::ReferencePlatform() {
registerKernelFactory(CalcCustomCPPForceKernel::Name(), factory); registerKernelFactory(CalcCustomCPPForceKernel::Name(), factory);
registerKernelFactory(CalcCustomCVForceKernel::Name(), factory); registerKernelFactory(CalcCustomCVForceKernel::Name(), factory);
registerKernelFactory(CalcATMForceKernel::Name(), factory); registerKernelFactory(CalcATMForceKernel::Name(), factory);
registerKernelFactory(CalcOrientationRestraintForceKernel::Name(), factory);
registerKernelFactory(CalcRGForceKernel::Name(), factory); registerKernelFactory(CalcRGForceKernel::Name(), factory);
registerKernelFactory(CalcRMSDForceKernel::Name(), factory); registerKernelFactory(CalcRMSDForceKernel::Name(), factory);
registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory); registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
......
/* Portions copyright (c) 2025 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include "ReferenceOrientationRestraintForce.h"
#include "jama_eig.h"
using namespace OpenMM;
using namespace std;
ReferenceOrientationRestraintForce::ReferenceOrientationRestraintForce(double k, vector<Vec3>& referencePos, vector<int>& particles) :
k(k), referencePos(referencePos), particles(particles) {
}
ReferenceOrientationRestraintForce::~ReferenceOrientationRestraintForce() {
}
double ReferenceOrientationRestraintForce::calculateIxn(vector<Vec3>& atomCoordinates, vector<Vec3>& forces) const {
// Find the optimal transformation using the algorithm described in Coutsias et al,
// "Using quaternions to calculate RMSD" (doi: 10.1002/jcc.20110). First subtract
// the centroid from the atom positions. The reference positions have already been centered.
int numParticles = particles.size();
Vec3 center;
for (int i : particles)
center += atomCoordinates[i];
center /= numParticles;
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++)
positions[i] = atomCoordinates[particles[i]]-center;
// Compute the correlation matrix.
double R[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
for (int k = 0; k < numParticles; k++) {
int index = particles[k];
R[i][j] += positions[k][j]*referencePos[index][i];
}
// Compute the F matrix.
Array2D<double> F(4, 4);
F[0][0] = R[0][0] + R[1][1] + R[2][2];
F[1][0] = R[1][2] - R[2][1];
F[2][0] = R[2][0] - R[0][2];
F[3][0] = R[0][1] - R[1][0];
F[0][1] = R[1][2] - R[2][1];
F[1][1] = R[0][0] - R[1][1] - R[2][2];
F[2][1] = R[0][1] + R[1][0];
F[3][1] = R[0][2] + R[2][0];
F[0][2] = R[2][0] - R[0][2];
F[1][2] = R[0][1] + R[1][0];
F[2][2] = -R[0][0] + R[1][1] - R[2][2];
F[3][2] = R[1][2] + R[2][1];
F[0][3] = R[0][1] - R[1][0];
F[1][3] = R[0][2] + R[2][0];
F[2][3] = R[1][2] + R[2][1];
F[3][3] = -R[0][0] - R[1][1] + R[2][2];
// Find the maximum eigenvalue and eigenvector.
JAMA::Eigenvalue<double> eigen(F);
Array1D<double> values;
eigen.getRealEigenvalues(values);
Array2D<double> vectors;
eigen.getV(vectors);
// Construct the quaternion and use it to compute the energy.
double q[] = {vectors[0][3], vectors[1][3], vectors[2][3], vectors[3][3]};
double energy = 2*k*(1.0-q[0]*q[0]);
// Compute the forces. The algorithm is modelled after the calculation of the
// orientationAngle CV in Colvars. I can't find any published source for it.
if (q[0]*q[0] < 1.0) {
double theta = 2*asin(sqrt(1.0-q[0]*q[0]));
double dxdq = 4.0*k*sin(theta/2)*cos(theta/2)/sqrt(1.0-q[0]*q[0]);
if (vectors[0][3] > 0)
dxdq = -dxdq;
for (int index = 0; index < numParticles; index++) {
const Vec3& p = referencePos[particles[index]];
Vec3 ds[4][4] = {
{Vec3(p[0], p[1], p[2]), Vec3(0.0, -p[2], p[1]), Vec3(p[2], 0.0, -p[0]), Vec3(-p[1], p[0], 0.0)},
{Vec3(0.0, -p[2], p[1]), Vec3(p[0], -p[1], -p[2]), Vec3(p[1], p[0], 0.0), Vec3(p[2], 0.0, p[0])},
{Vec3(p[2], 0.0, -p[0]), Vec3(p[1], p[0], 0.0), Vec3(-p[0], p[1], -p[2]), Vec3(0.0, p[2], p[1])},
{Vec3(-p[1], p[0], 0.0), Vec3(p[2], 0.0, p[0]), Vec3(0.0, p[2], p[1]), Vec3(-p[0], -p[1], p[2])}
};
Vec3 dq;
for (int i = 0; i < 4; i++) {
for (int j = 0; j < 4; j++) {
dq += ((vectors[i][2]*vectors[0][2]) / (values[3]-values[2]) +
(vectors[i][1]*vectors[0][1]) / (values[3]-values[1]) +
(vectors[i][0]*vectors[0][0]) / (values[3]-values[0])) * vectors[j][3] * ds[i][j];
}
}
forces[particles[index]] -= dxdq*dq;
}
}
return energy;
}
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ReferenceTests.h"
#include "TestOrientationRestraintForce.h"
void runPlatformTests() {
}
#ifndef OPENMM_ORIENTATIONRESTRAINTFORCE_PROXY_H_
#define OPENMM_ORIENTATIONRESTRAINTFORCE_PROXY_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/windowsExport.h"
#include "openmm/serialization/SerializationProxy.h"
namespace OpenMM {
/**
* This is a proxy for serializing OrientationRestraintForce objects.
*/
class OPENMM_EXPORT OrientationRestraintForceProxy : public SerializationProxy {
public:
OrientationRestraintForceProxy();
void serialize(const void* object, SerializationNode& node) const;
void* deserialize(const SerializationNode& node) const;
};
} // namespace OpenMM
#endif /*OPENMM_ORIENTATIONRESTRAINTFORCE_PROXY_H_*/
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/serialization/OrientationRestraintForceProxy.h"
#include "openmm/serialization/SerializationNode.h"
#include "openmm/Force.h"
#include "openmm/OrientationRestraintForce.h"
#include <sstream>
using namespace OpenMM;
using namespace std;
OrientationRestraintForceProxy::OrientationRestraintForceProxy() : SerializationProxy("OrientationRestraintForce") {
}
void OrientationRestraintForceProxy::serialize(const void* object, SerializationNode& node) const {
node.setIntProperty("version", 0);
const OrientationRestraintForce& force = *reinterpret_cast<const OrientationRestraintForce*>(object);
node.setIntProperty("forceGroup", force.getForceGroup());
node.setStringProperty("name", force.getName());
node.setDoubleProperty("k", force.getK());
SerializationNode& positionsNode = node.createChildNode("ReferencePositions");
for (const Vec3& pos : force.getReferencePositions())
positionsNode.createChildNode("Position").setDoubleProperty("x", pos[0]).setDoubleProperty("y", pos[1]).setDoubleProperty("z", pos[2]);
SerializationNode& particlesNode = node.createChildNode("Particles");
for (int i : force.getParticles())
particlesNode.createChildNode("Particle").setIntProperty("index", i);
}
void* OrientationRestraintForceProxy::deserialize(const SerializationNode& node) const {
int version = node.getIntProperty("version");
if (version != 0)
throw OpenMMException("Unsupported version number");
OrientationRestraintForce* force = NULL;
try {
double k = node.getDoubleProperty("k");
vector<Vec3> positions;
for (auto& pos : node.getChildNode("ReferencePositions").getChildren())
positions.push_back(Vec3(pos.getDoubleProperty("x"), pos.getDoubleProperty("y"), pos.getDoubleProperty("z")));
vector<int> particles;
for (auto& particle : node.getChildNode("Particles").getChildren())
particles.push_back(particle.getIntProperty("index"));
force = new OrientationRestraintForce(k, positions, particles);
force->setForceGroup(node.getIntProperty("forceGroup", 0));
force->setName(node.getStringProperty("name", force->getName()));
return force;
}
catch (...) {
if (force != NULL)
delete force;
throw;
}
}
...@@ -61,6 +61,7 @@ ...@@ -61,6 +61,7 @@
#include "openmm/MonteCarloMembraneBarostat.h" #include "openmm/MonteCarloMembraneBarostat.h"
#include "openmm/NonbondedForce.h" #include "openmm/NonbondedForce.h"
#include "openmm/NoseHooverIntegrator.h" #include "openmm/NoseHooverIntegrator.h"
#include "openmm/OrientationRestraintForce.h"
#include "openmm/PeriodicTorsionForce.h" #include "openmm/PeriodicTorsionForce.h"
#include "openmm/QTBIntegrator.h" #include "openmm/QTBIntegrator.h"
#include "openmm/RBTorsionForce.h" #include "openmm/RBTorsionForce.h"
...@@ -105,6 +106,7 @@ ...@@ -105,6 +106,7 @@
#include "openmm/serialization/MonteCarloMembraneBarostatProxy.h" #include "openmm/serialization/MonteCarloMembraneBarostatProxy.h"
#include "openmm/serialization/NonbondedForceProxy.h" #include "openmm/serialization/NonbondedForceProxy.h"
#include "openmm/serialization/NoseHooverIntegratorProxy.h" #include "openmm/serialization/NoseHooverIntegratorProxy.h"
#include "openmm/serialization/OrientationRestraintForceProxy.h"
#include "openmm/serialization/PeriodicTorsionForceProxy.h" #include "openmm/serialization/PeriodicTorsionForceProxy.h"
#include "openmm/serialization/QTBIntegratorProxy.h" #include "openmm/serialization/QTBIntegratorProxy.h"
#include "openmm/serialization/RBTorsionForceProxy.h" #include "openmm/serialization/RBTorsionForceProxy.h"
...@@ -170,6 +172,7 @@ extern "C" void registerSerializationProxies() { ...@@ -170,6 +172,7 @@ extern "C" void registerSerializationProxies() {
SerializationProxy::registerProxy(typeid(MonteCarloMembraneBarostat), new MonteCarloMembraneBarostatProxy()); SerializationProxy::registerProxy(typeid(MonteCarloMembraneBarostat), new MonteCarloMembraneBarostatProxy());
SerializationProxy::registerProxy(typeid(NonbondedForce), new NonbondedForceProxy()); SerializationProxy::registerProxy(typeid(NonbondedForce), new NonbondedForceProxy());
SerializationProxy::registerProxy(typeid(NoseHooverIntegrator), new NoseHooverIntegratorProxy()); SerializationProxy::registerProxy(typeid(NoseHooverIntegrator), new NoseHooverIntegratorProxy());
SerializationProxy::registerProxy(typeid(OrientationRestraintForce), new OrientationRestraintForceProxy());
SerializationProxy::registerProxy(typeid(PeriodicTorsionForce), new PeriodicTorsionForceProxy()); SerializationProxy::registerProxy(typeid(PeriodicTorsionForce), new PeriodicTorsionForceProxy());
SerializationProxy::registerProxy(typeid(QTBIntegrator), new QTBIntegratorProxy()); SerializationProxy::registerProxy(typeid(QTBIntegrator), new QTBIntegratorProxy());
SerializationProxy::registerProxy(typeid(RBTorsionForce), new RBTorsionForceProxy()); SerializationProxy::registerProxy(typeid(RBTorsionForce), new RBTorsionForceProxy());
......
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/OrientationRestraintForce.h"
#include "openmm/serialization/XmlSerializer.h"
#include <iostream>
#include <sstream>
using namespace OpenMM;
using namespace std;
void testSerialization() {
// Create a Force.
vector<Vec3> refPos;
for (int i = 0; i < 10; i++)
refPos.push_back(Vec3(i/5.0, i*1.2, i*i/3.5));
vector<int> particles;
for (int i = 0; i < 5; i++)
particles.push_back(i*i);
OrientationRestraintForce force(2.3, refPos, particles);
force.setForceGroup(3);
force.setName("custom name");
// Serialize and then deserialize it.
stringstream buffer;
XmlSerializer::serialize<OrientationRestraintForce>(&force, "Force", buffer);
OrientationRestraintForce* copy = XmlSerializer::deserialize<OrientationRestraintForce>(buffer);
// Compare the two forces to see if they are identical.
OrientationRestraintForce& force2 = *copy;
ASSERT_EQUAL(force.getK(), force2.getK());
ASSERT_EQUAL(force.getForceGroup(), force2.getForceGroup());
ASSERT_EQUAL(force.getName(), force2.getName());
ASSERT_EQUAL(force.getReferencePositions().size(), force2.getReferencePositions().size());
for (int i = 0; i < force.getReferencePositions().size(); i++)
ASSERT_EQUAL_VEC(force.getReferencePositions()[i], force2.getReferencePositions()[i], 0.0);
ASSERT_EQUAL(force.getParticles().size(), force2.getParticles().size());
for (int i = 0; i < force.getParticles().size(); i++)
ASSERT_EQUAL(force.getParticles()[i], force2.getParticles()[i]);
}
int main() {
try {
testSerialization();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/NonbondedForce.h"
#include "openmm/OrientationRestraintForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <cmath>
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
void testOrientationRestraint() {
const int numParticles = 20;
const double k = 3.5;
System system;
vector<Vec3> referencePos(numParticles);
vector<Vec3> positions(numParticles);
vector<int> particles;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
Vec3 center;
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
referencePos[i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*10;
if (i%5 != 0) {
particles.push_back(i);
center += referencePos[i];
}
}
center /= particles.size();
OrientationRestraintForce* force = new OrientationRestraintForce(k, referencePos, particles);
system.addForce(force);
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
// Randomly transform the reference positions and see if the energy is correct.
for (int i = 0; i < 20; i++) {
// Select a random axis, angle, and translation, and transform the particles from the reference positions.
Vec3 translation(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
Vec3 axis(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
axis /= sqrt(axis.dot(axis));
double angle = 3*genrand_real2(sfmt);
for (int j : particles) {
Vec3 p = referencePos[j]-center;
Vec3 cross1 = axis.cross(p);
p += sin(angle)*cross1 + (1-cos(angle))*axis.cross(cross1);
p += translation;
positions[j] = p;
}
context.setPositions(positions);
State state = context.getState(State::Energy);
double s = sin(angle/2);
ASSERT_EQUAL_TOL(2*k*s*s, state.getPotentialEnergy(), 1e-6);
}
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
for (int i = 0; i < 20; i++) {
for (int j = 0; j < numParticles; ++j)
positions[j] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*10;
context.setPositions(positions);
vector<Vec3> forces = context.getState(State::Forces).getForces();
double norm = 0.0;
for (int i = 0; i < forces.size(); ++i)
norm += forces[i].dot(forces[i]);
norm = sqrt(norm);
const double stepSize = 0.01;
double step = 0.5*stepSize/norm;
vector<Vec3> positions2(numParticles), positions3(numParticles);
for (int i = 0; i < positions.size(); ++i) {
Vec3 p = positions[i];
Vec3 f = forces[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
}
context.setPositions(positions2);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-3);
}
// When the current positions equal the reference positions, all forces should be zero.
context.setPositions(referencePos);
vector<Vec3> forces = context.getState(State::Forces).getForces();
Vec3 zero;
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(zero, forces[i], 1e-4);
// Check that updateParametersInContext() works correctly.
context.setPositions(positions);
double e1 = context.getState(State::Energy).getPotentialEnergy();
force->setK(2*k);
force->updateParametersInContext(context);
double e2 = context.getState(State::Energy).getPotentialEnergy();
force->setReferencePositions(positions);
force->updateParametersInContext(context);
double e3 = context.getState(State::Energy).getPotentialEnergy();
ASSERT_EQUAL_TOL(2*e1, e2, 1e-6);
ASSERT_EQUAL_TOL(0.0, e3, 1e-6);
}
void testEnergyConservation() {
const int numParticles = 50;
System system;
vector<Vec3> referencePos(numParticles);
vector<Vec3> positions(numParticles);
vector<int> particles;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
NonbondedForce* nb = new NonbondedForce(); // Add a nonbonded force to activate reordering on the GPU
nb->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
system.addForce(nb);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
nb->addParticle(0.0, 0.1, 0.01);
positions[i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*5;
referencePos[i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*5;
if (genrand_real2(sfmt) < 0.5)
particles.push_back(i);
}
OrientationRestraintForce* force = new OrientationRestraintForce(10.0, referencePos, particles);
system.addForce(force);
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0, 0);
integrator.step(5);
State initialState = context.getState(State::Energy);
double energy = initialState.getPotentialEnergy()+initialState.getKineticEnergy();
for (int i = 0; i < 100; i++) {
integrator.step(5);
State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy()+state.getKineticEnergy(), 1e-4);
}
// If we modify the reference positions, the energy should change.
for (int i = 0; i < numParticles; ++i)
referencePos[i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*5;
force->setReferencePositions(referencePos);
force->updateParametersInContext(context);
State state2 = context.getState(State::Energy);
double energy2 = state2.getPotentialEnergy()+state2.getKineticEnergy();
ASSERT(fabs(energy-energy2) > 1e-3);
// Make sure it's still conserved.
for (int i = 0; i < 100; i++) {
integrator.step(5);
State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(energy2, state.getPotentialEnergy()+state.getKineticEnergy(), 1e-4);
}
}
void runPlatformTests();
int main(int argc, char* argv[]) {
try {
initializeTests(argc, argv);
testOrientationRestraint();
testEnergyConservation();
runPlatformTests();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
...@@ -32,6 +32,7 @@ ...@@ -32,6 +32,7 @@
#include "openmm/internal/AssertionUtilities.h" #include "openmm/internal/AssertionUtilities.h"
#include "openmm/RMSDForce.h" #include "openmm/RMSDForce.h"
#include "openmm/Context.h" #include "openmm/Context.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h" #include "openmm/System.h"
#include "openmm/VerletIntegrator.h" #include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h" #include "sfmt/SFMT.h"
...@@ -155,12 +156,66 @@ void testRMSD() { ...@@ -155,12 +156,66 @@ void testRMSD() {
ASSERT(rmsd1 > 0.9*estimate); ASSERT(rmsd1 > 0.9*estimate);
} }
void testEnergyConservation() {
const int numParticles = 50;
System system;
vector<Vec3> referencePos(numParticles);
vector<Vec3> positions(numParticles);
vector<int> particles;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
NonbondedForce* nb = new NonbondedForce(); // Add a nonbonded force to activate reordering on the GPU
nb->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
system.addForce(nb);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
nb->addParticle(0.0, 0.1, 0.01);
positions[i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*5;
referencePos[i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*5;
if (genrand_real2(sfmt) < 0.5)
particles.push_back(i);
}
RMSDForce* force = new RMSDForce(referencePos, particles);
system.addForce(force);
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0, 0);
integrator.step(5);
State initialState = context.getState(State::Energy);
double energy = initialState.getPotentialEnergy()+initialState.getKineticEnergy();
for (int i = 0; i < 100; i++) {
integrator.step(5);
State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy()+state.getKineticEnergy(), 1e-4);
}
// If we modify the reference positions, the energy should change.
for (int i = 0; i < numParticles; ++i)
referencePos[i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*5;
force->setReferencePositions(referencePos);
force->updateParametersInContext(context);
State state2 = context.getState(State::Energy);
double energy2 = state2.getPotentialEnergy()+state2.getKineticEnergy();
ASSERT(fabs(energy-energy2) > 1e-3);
// Make sure it's still conserved.
for (int i = 0; i < 100; i++) {
integrator.step(5);
State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(energy2, state.getPotentialEnergy()+state.getKineticEnergy(), 1e-4);
}
}
void runPlatformTests(); void runPlatformTests();
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
try { try {
initializeTests(argc, argv); initializeTests(argc, argv);
testRMSD(); testRMSD();
testEnergyConservation();
runPlatformTests(); runPlatformTests();
} }
catch(const exception& e) { catch(const exception& e) {
......
...@@ -491,6 +491,11 @@ UNITS = { ...@@ -491,6 +491,11 @@ UNITS = {
("RMSDForce", "getReferencePositions") : ("unit.nanometer", ()), ("RMSDForce", "getReferencePositions") : ("unit.nanometer", ()),
("RMSDForce", "setReferencePositions") : (None, ("unit.nanometer",)), ("RMSDForce", "setReferencePositions") : (None, ("unit.nanometer",)),
("RMSDForce", "getParticles") : (None, ()), ("RMSDForce", "getParticles") : (None, ()),
("OrientationRestraintForce", "getReferencePositions") : ("unit.nanometer", ()),
("OrientationRestraintForce", "setReferencePositions") : (None, ("unit.nanometer",)),
("OrientationRestraintForce", "getParticles") : (None, ()),
("OrientationRestraintForce", "getK") : ("unit.kilojoules_per_mole", ()),
("OrientationRestraintForce", "setK") : (None, ("unit.kilojoules_per_mole",)),
("RGForce", "getParticles") : (None, ()), ("RGForce", "getParticles") : (None, ()),
("NoseHooverChain", "NoseHooverChain") : (None, ("unit.kelvin", "unit.kelvin", "unit.picosecond**-1", "unit.picosecond**-1", None, None, None, None, None, None, None)), ("NoseHooverChain", "NoseHooverChain") : (None, ("unit.kelvin", "unit.kelvin", "unit.picosecond**-1", "unit.picosecond**-1", None, None, None, None, None, None, None)),
("NoseHooverChain", "getThermostatedPairs") : (None, ()), ("NoseHooverChain", "getThermostatedPairs") : (None, ()),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment