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

Drude integrators can find DrudeForce inside another force (#5054)

parent 54fa65ed
......@@ -36,7 +36,7 @@ public:
std::vector<std::pair<int, int> > getBondedParticles() const;
void updateParametersInContext(ContextImpl& context);
void getPerturbationEnergy(ContextImpl& context, double& u1, double& u0, double& energy);
std::vector<const Force*> getContainedForces() const;
private:
const ATMForce& owner;
Kernel kernel;
......
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2022 Stanford University and the Authors. *
* Portions copyright (c) 2008-2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -66,6 +66,7 @@ public:
void getCollectiveVariableValues(ContextImpl& context, std::vector<double>& values);
Context& getInnerContext();
void updateParametersInContext(ContextImpl& context);
std::vector<const Force*> getContainedForces() const;
private:
const CustomCVForce& owner;
Kernel kernel;
......
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -112,6 +112,11 @@ public:
virtual std::vector<std::pair<int, int> > getBondedParticles() const {
return std::vector<std::pair<int, int> >(0);
}
/**
* Some Forces may themselves contain other Forces. This returns any Forces contained by
* this one. The default implementation returns an empty vector.
*/
virtual std::vector<const Force*> getContainedForces() const;
protected:
/**
* Get the ContextImpl corresponding to a Context.
......
......@@ -9,7 +9,7 @@
* https://github.com/Gallicchio-Lab/openmm-atmmetaforce-plugin *
* with support from the National Science Foundation CAREER 1750511 *
* *
* Portions copyright (c) 2021-2024 by the Authors *
* Portions copyright (c) 2021-2025 by the Authors *
* Authors: Emilio Gallicchio *
* Contributors: Peter Eastman *
* *
......@@ -203,3 +203,10 @@ void ATMForceImpl::getPerturbationEnergy(ContextImpl& context, double& u1, doubl
u1 = state1Energy;
energy = combinedEnergy;
}
vector<const Force*> ATMForceImpl::getContainedForces() const {
vector<const Force*> forces;
for (int i = 0; i < owner.getNumForces(); i++)
forces.push_back(&owner.getForce(i));
return forces;
}
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2022 Stanford University and the Authors. *
* Portions copyright (c) 2008-2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -127,3 +127,10 @@ void CustomCVForceImpl::updateParametersInContext(ContextImpl& context) {
kernel.getAs<CalcCustomCVForceKernel>().copyParametersToContext(context, owner);
context.systemChanged();
}
vector<const Force*> CustomCVForceImpl::getContainedForces() const {
vector<const Force*> forces;
for (int i = 0; i < owner.getNumCollectiveVariables(); i++)
forces.push_back(&owner.getCollectiveVariable(i));
return forces;
}
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2017 Stanford University and the Authors. *
* Portions copyright (c) 2017-2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -44,3 +44,7 @@ void ForceImpl::updateContextState(ContextImpl& context, bool& forcesInvalid) {
void ForceImpl::updateContextState(ContextImpl& context) {
}
vector<const Force*> ForceImpl::getContainedForces() const {
return {};
}
#ifndef OPENMM_DRUDEHELPERS_H_
#define OPENMM_DRUDEHELPERS_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) 2008-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/DrudeForce.h"
#include "openmm/System.h"
#include "openmm/internal/ContextImpl.h"
namespace OpenMM {
/**
* Find the DrudeForce contained in the System.
*/
const DrudeForce* getDrudeForce(const ContextImpl& context);
/**
* Identify normal particles (not part of a pair) and Drude particle pairs.
*/
void findParticlesAndPairs(const ContextImpl& context, std::vector<int>& normalParticles, std::vector<std::pair<int, int> >& pairParticles);
std::vector<Vec3> assignDrudeVelocities(const ContextImpl& context, double temperature, double drudeTemperature, int randomSeed);
/**
* Computes the instantaneous temperatures of the system and the internal Drude motion and returns a pair (T_system, T_drude)
*/
std::pair<double, double> computeTemperaturesFromVelocities(const ContextImpl& context, const std::vector<Vec3>& velocities);
}
#endif /*OPENMM_DRUDEHELPERS_H_*/
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2023 Stanford University and the Authors. *
* Portions copyright (c) 2008-2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,14 +29,14 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/DrudeHelpers.h"
#include "sfmt/SFMT.h"
#include "SimTKOpenMMRealType.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/DrudeForce.h"
#include "openmm/DrudeLangevinIntegrator.h"
#include "openmm/Context.h"
#include "openmm/System.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ForceImpl.h"
#include <set>
......@@ -45,20 +45,42 @@ using namespace std;
namespace OpenMM {
/**
* Identify normal particles (not part of a pair) and Drude particle pairs.
* Find the DrudeForce contained in the System.
*/
void findParticlesAndPairs(const System &system, vector<int>& normalParticles, vector<pair<int, int> >& pairParticles) {
// Find the underlying Drude force object
const DrudeForce* drudeForce = NULL;
const DrudeForce* getDrudeForce(const ContextImpl& context) {
const DrudeForce* force = NULL;
const System& system = context.getSystem();
for (int i = 0; i < system.getNumForces(); i++)
if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
if (drudeForce == NULL)
drudeForce = dynamic_cast<const DrudeForce*>(&system.getForce(i));
if (force == NULL)
force = dynamic_cast<const DrudeForce*>(&system.getForce(i));
else
throw OpenMMException("The System contains multiple DrudeForces");
}
if (force != NULL)
return force;
// The System doesn't directly contain a DrudeForce. Look for one inside another force.
for (const ForceImpl* impl : context.getForceImpls())
for (const Force* f : impl->getContainedForces())
if (dynamic_cast<const DrudeForce*>(f) != NULL) {
if (force == NULL)
force = dynamic_cast<const DrudeForce*>(f);
else
throw OpenMMException("The System contains multiple DrudeForces");
}
if (drudeForce == NULL)
if (force != NULL)
return force;
throw OpenMMException("The System does not contain a DrudeForce");
}
/**
* Identify normal particles (not part of a pair) and Drude particle pairs.
*/
void findParticlesAndPairs(const ContextImpl& context, vector<int>& normalParticles, vector<pair<int, int> >& pairParticles) {
const System& system = context.getSystem();
const DrudeForce* drudeForce = getDrudeForce(context);
// Figure out which particles are individual and which are Drude pairs
set<int> particles;
......@@ -75,10 +97,11 @@ void findParticlesAndPairs(const System &system, vector<int>& normalParticles, v
normalParticles.insert(normalParticles.begin(), particles.begin(), particles.end());
}
vector<Vec3> assignDrudeVelocities(const System &system, double temperature, double drudeTemperature, int randomSeed) {
vector<Vec3> assignDrudeVelocities(const ContextImpl& context, double temperature, double drudeTemperature, int randomSeed) {
const System& system = context.getSystem();
vector<int> normalParticles;
vector<pair<int, int> > pairParticles;
findParticlesAndPairs(system, normalParticles, pairParticles);
findParticlesAndPairs(context, normalParticles, pairParticles);
// Generate the list of Gaussian random numbers.
OpenMM_SFMT::SFMT sfmt;
......@@ -135,10 +158,11 @@ vector<Vec3> assignDrudeVelocities(const System &system, double temperature, dou
/**
* Computes the instantaneous temperatures of the system and the internal Drude motion and returns a pair (T_system, T_drude)
*/
pair<double, double> computeTemperaturesFromVelocities(const System& system, const vector<Vec3>& velocities) {
pair<double, double> computeTemperaturesFromVelocities(const ContextImpl& context, const vector<Vec3>& velocities) {
const System& system = context.getSystem();
vector<int> normalParticles;
vector<pair<int, int> > pairParticles;
findParticlesAndPairs(system, normalParticles, pairParticles);
findParticlesAndPairs(context, normalParticles, pairParticles);
double energy = 0.0, drudeEnergy = 0;
int dof = 0, drudeDof = 0;
......
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2021 Stanford University and the Authors. *
* Portions copyright (c) 2008-2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -34,17 +34,13 @@
#include "openmm/DrudeIntegrator.h"
#include "openmm/OpenMMException.h"
#include "openmm/System.h"
#include "openmm/internal/DrudeHelpers.h"
#include <set>
using namespace OpenMM;
namespace OpenMM {
extern std::vector<Vec3> assignDrudeVelocities(const System &system, double temperature, double drudeTemperature, int randomSeed);
}
std::vector<Vec3> DrudeIntegrator::getVelocitiesForTemperature(const System &system, double temperature, int randomSeedIn) const {
return assignDrudeVelocities(system, temperature, drudeTemperature, randomSeedIn);
return assignDrudeVelocities(*context, temperature, drudeTemperature, randomSeedIn);
}
void DrudeIntegrator::setDrudeTemperature(double temp) {
......
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2022 Stanford University and the Authors. *
* Portions copyright (c) 2008-2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -32,8 +32,9 @@
#include "openmm/DrudeLangevinIntegrator.h"
#include "openmm/Context.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/DrudeKernels.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/DrudeHelpers.h"
#include <cmath>
#include <ctime>
#include <set>
......@@ -44,10 +45,6 @@ using std::string;
using std::vector;
using std::pair;
namespace OpenMM {
pair<double, double> computeTemperaturesFromVelocities(const System& system, const vector<Vec3>& velocities);
}
DrudeLangevinIntegrator::DrudeLangevinIntegrator(double temperature, double frictionCoeff, double drudeTemperature, double drudeFrictionCoeff, double stepSize) : DrudeIntegrator(stepSize) {
setTemperature(temperature);
setFriction(frictionCoeff);
......@@ -62,17 +59,7 @@ DrudeLangevinIntegrator::DrudeLangevinIntegrator(double temperature, double fric
void DrudeLangevinIntegrator::initialize(ContextImpl& contextRef) {
if (owner != NULL && &contextRef.getOwner() != owner)
throw OpenMMException("This Integrator is already bound to a context");
const DrudeForce* force = NULL;
const System& system = contextRef.getSystem();
for (int i = 0; i < system.getNumForces(); i++)
if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
if (force == NULL)
force = dynamic_cast<const DrudeForce*>(&system.getForce(i));
else
throw OpenMMException("The System contains multiple DrudeForces");
}
if (force == NULL)
throw OpenMMException("The System does not contain a DrudeForce");
const DrudeForce* force = getDrudeForce(contextRef);
context = &contextRef;
owner = &contextRef.getOwner();
kernel = context->getPlatform().createKernel(IntegrateDrudeLangevinStepKernel::Name(), contextRef);
......@@ -127,7 +114,7 @@ double DrudeLangevinIntegrator::computeSystemTemperature() {
context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
vector<Vec3> velocities;
context->computeShiftedVelocities(getVelocityTimeOffset(), velocities);
return computeTemperaturesFromVelocities(context->getSystem(), velocities).first;
return computeTemperaturesFromVelocities(*context, velocities).first;
}
double DrudeLangevinIntegrator::computeDrudeTemperature() {
......@@ -136,6 +123,6 @@ double DrudeLangevinIntegrator::computeDrudeTemperature() {
context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
vector<Vec3> velocities;
context->computeShiftedVelocities(getVelocityTimeOffset(), velocities);
return computeTemperaturesFromVelocities(context->getSystem(), velocities).second;
return computeTemperaturesFromVelocities(*context, velocities).second;
}
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2019-2022 Stanford University and the Authors. *
* Portions copyright (c) 2019-2025 Stanford University and the Authors. *
* Authors: Andreas Krämer and Andrew C. Simmonett *
* Contributors: Peter Eastman *
* *
......@@ -33,11 +33,12 @@
#include "SimTKOpenMMRealType.h"
#include "openmm/Context.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/DrudeKernels.h"
#include "openmm/DrudeForce.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/kernels.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/DrudeHelpers.h"
#include <ctime>
#include <iostream>
#include <string>
......@@ -48,12 +49,6 @@ using std::string;
using std::vector;
using std::pair;
namespace OpenMM {
extern std::vector<Vec3> assignDrudeVelocities(const System &system, double temperature, double drudeTemperature, int randomSeed);
pair<double, double> computeTemperaturesFromVelocities(const System& system, const vector<Vec3>& velocities);
}
DrudeNoseHooverIntegrator::DrudeNoseHooverIntegrator(double temperature, double collisionFrequency,
double drudeTemperature, double drudeCollisionFrequency,
double stepSize, int chainLength, int numMTS, int numYoshidaSuzuki) :
......@@ -90,22 +85,14 @@ void DrudeNoseHooverIntegrator::initialize(ContextImpl& contextRef) {
// check for drude particles and build the Nose-Hoover Chains
const System& system = context->getSystem();
const DrudeForce* drudeForce = NULL;
const DrudeForce* drudeForce = getDrudeForce(contextRef);
bool hasCMMotionRemover = false;
for (int i = 0; i < system.getNumForces(); i++){
if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
if (drudeForce == NULL)
drudeForce = dynamic_cast<const DrudeForce*>(&system.getForce(i));
else
throw OpenMMException("The System contains multiple DrudeForces");
}
if (dynamic_cast<const CMMotionRemover*>(&system.getForce(i))) {
hasCMMotionRemover = true;
}
}
if (drudeForce == NULL)
throw OpenMMException("The System does not contain a DrudeForce");
bool hasCMMotionRemover = false;
for (int i = 0; i < system.getNumForces(); i++)
if (dynamic_cast<const CMMotionRemover*>(&system.getForce(i)))
hasCMMotionRemover = true;
if (!hasCMMotionRemover) {
std::cout << "Warning: Did not find a center-of-mass motion remover in the system. "
"This is problematic when using Drude." << std::endl;
......@@ -154,7 +141,7 @@ double DrudeNoseHooverIntegrator::computeSystemTemperature() {
context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
vector<Vec3> velocities;
context->computeShiftedVelocities(getVelocityTimeOffset(), velocities);
return computeTemperaturesFromVelocities(context->getSystem(), velocities).first;
return computeTemperaturesFromVelocities(*context, velocities).first;
}
double DrudeNoseHooverIntegrator::computeDrudeTemperature() {
......@@ -170,6 +157,6 @@ double DrudeNoseHooverIntegrator::computeDrudeTemperature() {
std::vector<Vec3> DrudeNoseHooverIntegrator::getVelocitiesForTemperature(const System &system, double temperature,
int randomSeedIn) const {
return assignDrudeVelocities(system, temperature, drudeTemperature, randomSeedIn);
return assignDrudeVelocities(*context, temperature, drudeTemperature, randomSeedIn);
}
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2023 Stanford University and the Authors. *
* Portions copyright (c) 2008-2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -30,10 +30,11 @@
* -------------------------------------------------------------------------- */
#include "openmm/DrudeSCFIntegrator.h"
#include "openmm/DrudeKernels.h"
#include "openmm/Context.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/DrudeKernels.h"
#include "openmm/internal/DrudeHelpers.h"
#include <cmath>
#include <ctime>
#include <string>
......@@ -54,17 +55,7 @@ DrudeSCFIntegrator::DrudeSCFIntegrator(double stepSize) : DrudeIntegrator(stepSi
void DrudeSCFIntegrator::initialize(ContextImpl& contextRef) {
if (owner != NULL && &contextRef.getOwner() != owner)
throw OpenMMException("This Integrator is already bound to a context");
const DrudeForce* force = NULL;
const System& system = contextRef.getSystem();
for (int i = 0; i < system.getNumForces(); i++)
if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
if (force == NULL)
force = dynamic_cast<const DrudeForce*>(&system.getForce(i));
else
throw OpenMMException("The System contains multiple DrudeForces");
}
if (force == NULL)
throw OpenMMException("The System does not contain a DrudeForce");
const DrudeForce* force = getDrudeForce(contextRef);
if (getMaxDrudeDistance() != 0.0)
throw OpenMMException("DrudeSCFIntegrator does not currently support setting max Drude distance");
context = &contextRef;
......
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013-2016 Stanford University and the Authors. *
* Portions copyright (c) 2013-2025 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -31,6 +31,7 @@
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/CustomCVForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/Platform.h"
#include "openmm/System.h"
......@@ -44,7 +45,7 @@
using namespace OpenMM;
using namespace std;
void testSinglePair() {
void testSinglePair(bool contained) {
const double temperature = 300.0;
const double temperatureDrude = 10.0;
const double k = ONE_4PI_EPS0*1.5;
......@@ -60,6 +61,12 @@ void testSinglePair() {
system.addParticle(mass2);
DrudeForce* drude = new DrudeForce();
drude->addParticle(1, 0, -1, -1, -1, charge, alpha, 1, 1);
if (contained) {
CustomCVForce* cv = new CustomCVForce("x");
cv->addCollectiveVariable("x", drude);
system.addForce(cv);
}
else
system.addForce(drude);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
......@@ -297,10 +304,11 @@ void runPlatformTests();
int main(int argc, char* argv[]) {
try {
setupKernels(argc, argv);
testInitialTemperature();
testSinglePair();
testSinglePair(false);
testSinglePair(true);
testWater();
testForceEnergyConsistency();
testInitialTemperature();
runPlatformTests();
}
catch(const std::exception& e) {
......
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