Skip to content

GitLab

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

Point based geometric functions for custom forces (#3037) 

* Began implementing geometric functions on points

* Started common implementation of point functions

* Completed implementation of point functions for CustomCompoundBondForce

* Implemented point functions for CustomCentroidBondForce

* Implemented point functions for CustomManyParticleForce

* Use point functions to simplify implementation of custom forces

* Removed unnecessary code

* Fixed typo
parent 2da337e9
Show whitespace changes
Inline Side-by-side
Showing with 430 additions and 533 deletions
platforms/reference/include/ReferenceCustomCompoundBondIxn.h
View file @ a85c2428
...@@ -38,18 +38,12 @@ class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn { ...@@ -38,18 +38,12 @@ class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn {
private: private:
class ParticleTermInfo; class ParticleTermInfo;
class DistanceTermInfo;
class AngleTermInfo;
class DihedralTermInfo;
std::vector<std::vector<int> > bondAtoms; std::vector<std::vector<int> > bondAtoms;
CompiledExpressionSet expressionSet; CompiledExpressionSet expressionSet;
Lepton::CompiledExpression energyExpression; Lepton::CompiledExpression energyExpression;
std::vector<Lepton::CompiledExpression> energyParamDerivExpressions; std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
std::vector<int> bondParamIndex; std::vector<int> bondParamIndex;
std::vector<ParticleTermInfo> particleTerms; std::vector<ParticleTermInfo> particleTerms;
std::vector<DistanceTermInfo> distanceTerms;
std::vector<AngleTermInfo> angleTerms;
std::vector<DihedralTermInfo> dihedralTerms;
int numParameters; int numParameters;
bool usePeriodic; bool usePeriodic;
Vec3 boxVectors[3]; Vec3 boxVectors[3];
...@@ -83,9 +77,7 @@ class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn { ...@@ -83,9 +77,7 @@ class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
ReferenceCustomCompoundBondIxn(int numParticlesPerBond, const std::vector<std::vector<int> >& bondAtoms, const Lepton::ParsedExpression& energyExpression, ReferenceCustomCompoundBondIxn(int numParticlesPerBond, const std::vector<std::vector<int> >& bondAtoms, const Lepton::ParsedExpression& energyExpression,
const std::vector<std::string>& bondParameterNames, const std::map<std::string, std::vector<int> >& distances, const std::vector<std::string>& bondParameterNames, const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions);
const std::map<std::string, std::vector<int> >& angles, const std::map<std::string, std::vector<int> >& dihedrals,
const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -145,44 +137,6 @@ public: ...@@ -145,44 +137,6 @@ public:
} }
}; };
class ReferenceCustomCompoundBondIxn::DistanceTermInfo {
public:
std::string name;
int p1, p2, index;
Lepton::CompiledExpression forceExpression;
mutable double delta[ReferenceForce::LastDeltaRIndex];
DistanceTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression) :
name(name), p1(atoms[0]), p2(atoms[1]), forceExpression(forceExpression) {
}
};
class ReferenceCustomCompoundBondIxn::AngleTermInfo {
public:
std::string name;
int p1, p2, p3, index;
Lepton::CompiledExpression forceExpression;
mutable double delta1[ReferenceForce::LastDeltaRIndex];
mutable double delta2[ReferenceForce::LastDeltaRIndex];
AngleTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression) :
name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), forceExpression(forceExpression) {
}
};
class ReferenceCustomCompoundBondIxn::DihedralTermInfo {
public:
std::string name;
int p1, p2, p3, p4, index;
Lepton::CompiledExpression forceExpression;
mutable double delta1[ReferenceForce::LastDeltaRIndex];
mutable double delta2[ReferenceForce::LastDeltaRIndex];
mutable double delta3[ReferenceForce::LastDeltaRIndex];
mutable double cross1[3];
mutable double cross2[3];
DihedralTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression) :
name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), p4(atoms[3]), forceExpression(forceExpression) {
}
};
} // namespace OpenMM } // namespace OpenMM
#endif // __ReferenceCustomCompoundBondIxn_H__ #endif // __ReferenceCustomCompoundBondIxn_H__
platforms/reference/include/ReferenceCustomManyParticleIxn.h
View file @ a85c2428
/* Portions copyright (c) 2009-2021 Stanford University and Simbios.
/* Portions copyright (c) 2009-2018 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -40,13 +39,10 @@ class ReferenceCustomManyParticleIxn { ...@@ -40,13 +39,10 @@ class ReferenceCustomManyParticleIxn {
private: private:
class ParticleTermInfo; class ParticleTermInfo;
class DistanceTermInfo;
class AngleTermInfo;
class DihedralTermInfo;
int numParticlesPerSet, numPerParticleParameters, numTypes; int numParticlesPerSet, numPerParticleParameters, numTypes;
bool useCutoff, usePeriodic, centralParticleMode; bool useCutoff, usePeriodic, centralParticleMode;
double cutoffDistance; double cutoffDistance;
OpenMM::Vec3 periodicBoxVectors[3]; OpenMM::Vec3* periodicBoxVectors;
Lepton::ExpressionProgram energyExpression; Lepton::ExpressionProgram energyExpression;
std::vector<std::vector<std::string> > particleParamNames; std::vector<std::vector<std::string> > particleParamNames;
std::vector<std::set<int> > exclusions; std::vector<std::set<int> > exclusions;
...@@ -54,9 +50,6 @@ class ReferenceCustomManyParticleIxn { ...@@ -54,9 +50,6 @@ class ReferenceCustomManyParticleIxn {
std::vector<int> orderIndex; std::vector<int> orderIndex;
std::vector<std::vector<int> > particleOrder; std::vector<std::vector<int> > particleOrder;
std::vector<ParticleTermInfo> particleTerms; std::vector<ParticleTermInfo> particleTerms;
std::vector<DistanceTermInfo> distanceTerms;
std::vector<AngleTermInfo> angleTerms;
std::vector<DihedralTermInfo> dihedralTerms;
void loopOverInteractions(std::vector<int>& particles, int loopIndex, std::vector<OpenMM::Vec3>& atomCoordinates, void loopOverInteractions(std::vector<int>& particles, int loopIndex, std::vector<OpenMM::Vec3>& atomCoordinates,
std::vector<std::vector<double> >& particleParameters, std::map<std::string, double>& variables, std::vector<std::vector<double> >& particleParameters, std::map<std::string, double>& variables,
...@@ -151,44 +144,6 @@ public: ...@@ -151,44 +144,6 @@ public:
} }
}; };
class ReferenceCustomManyParticleIxn::DistanceTermInfo {
public:
std::string name;
int p1, p2;
Lepton::ExpressionProgram forceExpression;
mutable double delta[ReferenceForce::LastDeltaRIndex];
DistanceTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
name(name), p1(atoms[0]), p2(atoms[1]), forceExpression(forceExpression) {
}
};
class ReferenceCustomManyParticleIxn::AngleTermInfo {
public:
std::string name;
int p1, p2, p3;
Lepton::ExpressionProgram forceExpression;
mutable double delta1[ReferenceForce::LastDeltaRIndex];
mutable double delta2[ReferenceForce::LastDeltaRIndex];
AngleTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), forceExpression(forceExpression) {
}
};
class ReferenceCustomManyParticleIxn::DihedralTermInfo {
public:
std::string name;
int p1, p2, p3, p4;
Lepton::ExpressionProgram forceExpression;
mutable double delta1[ReferenceForce::LastDeltaRIndex];
mutable double delta2[ReferenceForce::LastDeltaRIndex];
mutable double delta3[ReferenceForce::LastDeltaRIndex];
mutable double cross1[3];
mutable double cross2[3];
DihedralTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), p4(atoms[3]), forceExpression(forceExpression) {
}
};
} // namespace OpenMM } // namespace OpenMM
#endif // __ReferenceCustomManyParticleIxn_H__ #endif // __ReferenceCustomManyParticleIxn_H__
platforms/reference/include/ReferenceKernels.h
View file @ a85c2428
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008-2020 Stanford University and the Authors. * * Portions copyright (c) 2008-2021 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -806,7 +806,6 @@ public: ...@@ -806,7 +806,6 @@ public:
*/ */
void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force); void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force);
private: private:
class PeriodicDistanceFunction;
int numParticles; int numParticles;
ReferenceCustomExternalIxn* ixn; ReferenceCustomExternalIxn* ixn;
std::vector<int> particles; std::vector<int> particles;
...@@ -816,16 +815,6 @@ private: ...@@ -816,16 +815,6 @@ private:
Vec3* boxVectors; Vec3* boxVectors;
}; };
class ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction : public Lepton::CustomFunction {
public:
Vec3** boxVectorHandle;
PeriodicDistanceFunction(Vec3** boxVectorHandle);
int getNumArguments() const;
double evaluate(const double* arguments) const;
double evaluateDerivative(const double* arguments, const int* derivOrder) const;
Lepton::CustomFunction* clone() const;
};
/** /**
* This kernel is invoked by CustomHbondForce to calculate the forces acting on the system. * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
*/ */
...@@ -904,6 +893,7 @@ private: ...@@ -904,6 +893,7 @@ private:
ReferenceCustomCentroidBondIxn* ixn; ReferenceCustomCentroidBondIxn* ixn;
std::vector<std::string> globalParameterNames, energyParamDerivNames; std::vector<std::string> globalParameterNames, energyParamDerivNames;
bool usePeriodic; bool usePeriodic;
Vec3* boxVectors;
}; };
/** /**
...@@ -943,6 +933,7 @@ private: ...@@ -943,6 +933,7 @@ private:
ReferenceCustomCompoundBondIxn* ixn; ReferenceCustomCompoundBondIxn* ixn;
std::vector<std::string> globalParameterNames, energyParamDerivNames; std::vector<std::string> globalParameterNames, energyParamDerivNames;
bool usePeriodic; bool usePeriodic;
Vec3* boxVectors;
}; };
/** /**
......
platforms/reference/include/ReferencePointFunctions.h 0 → 100644
View file @ a85c2428
#ifndef OPENMM_REFERENCEPOINTFUNCTIONS_H_
#define OPENMM_REFERENCEPOINTFUNCTIONS_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) 2021 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/Vec3.h"
#include "openmm/internal/windowsExport.h"
#include "lepton/CustomFunction.h"
namespace OpenMM {
/**
* This implements the pointdistance() function used in custom forces.
*/
class OPENMM_EXPORT ReferencePointDistanceFunction : public Lepton::CustomFunction {
public:
ReferencePointDistanceFunction(bool periodic, Vec3** boxVectorHandle);
int getNumArguments() const;
double evaluate(const double* arguments) const;
double evaluateDerivative(const double* arguments, const int* derivOrder) const;
Lepton::CustomFunction* clone() const;
private:
bool periodic;
Vec3** boxVectorHandle;
};
/**
* This implements the pointangle() function used in custom forces.
*/
class OPENMM_EXPORT ReferencePointAngleFunction : public Lepton::CustomFunction {
public:
ReferencePointAngleFunction(bool periodic, Vec3** boxVectorHandle);
int getNumArguments() const;
double evaluate(const double* arguments) const;
double evaluateDerivative(const double* arguments, const int* derivOrder) const;
Lepton::CustomFunction* clone() const;
private:
bool periodic;
Vec3** boxVectorHandle;
};
/**
* This implements the pointdihedral() function used in custom forces.
*/
class OPENMM_EXPORT ReferencePointDihedralFunction : public Lepton::CustomFunction {
public:
ReferencePointDihedralFunction(bool periodic, Vec3** boxVectorHandle);
int getNumArguments() const;
double evaluate(const double* arguments) const;
double evaluateDerivative(const double* arguments, const int* derivOrder) const;
Lepton::CustomFunction* clone() const;
private:
bool periodic;
Vec3** boxVectorHandle;
};
} // namespace OpenMM
#endif /*OPENMM_REFERENCEPOINTFUNCTIONS_H_*/
platforms/reference/src/ReferenceKernels.cpp
View file @ a85c2428
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008-2020 Stanford University and the Authors. * * Portions copyright (c) 2008-2021 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -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 "ReferencePointFunctions.h"
#include "ReferenceProperDihedralBond.h" #include "ReferenceProperDihedralBond.h"
#include "ReferenceRbDihedralBond.h" #include "ReferenceRbDihedralBond.h"
#include "ReferenceRMSDForce.h" #include "ReferenceRMSDForce.h"
...@@ -903,7 +904,6 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N ...@@ -903,7 +904,6 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
baseExceptionParams.resize(num14); baseExceptionParams.resize(num14);
for (int i = 0; i < numParticles; ++i) for (int i = 0; i < numParticles; ++i)
force.getParticleParameters(i, baseParticleParams[i][0], baseParticleParams[i][1], baseParticleParams[i][2]); force.getParticleParameters(i, baseParticleParams[i][0], baseParticleParams[i][1], baseParticleParams[i][2]);
this->exclusions = exclusions;
for (int i = 0; i < num14; ++i) { for (int i = 0; i < num14; ++i) {
int particle1, particle2; int particle1, particle2;
force.getExceptionParameters(nb14s[i], particle1, particle2, baseExceptionParams[i][0], baseExceptionParams[i][1], baseExceptionParams[i][2]); force.getExceptionParameters(nb14s[i], particle1, particle2, baseExceptionParams[i][0], baseExceptionParams[i][1], baseExceptionParams[i][2]);
...@@ -1525,48 +1525,6 @@ void ReferenceCalcCustomGBForceKernel::copyParametersToContext(ContextImpl& cont ...@@ -1525,48 +1525,6 @@ void ReferenceCalcCustomGBForceKernel::copyParametersToContext(ContextImpl& cont
} }
} }
ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::PeriodicDistanceFunction(Vec3** boxVectorHandle) : boxVectorHandle(boxVectorHandle) {
}
int ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::getNumArguments() const {
return 6;
}
double ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::evaluate(const double* arguments) const {
Vec3* boxVectors = *boxVectorHandle;
Vec3 delta = Vec3(arguments[0], arguments[1], arguments[2])-Vec3(arguments[3], arguments[4], arguments[5]);
delta -= boxVectors[2]*floor(delta[2]/boxVectors[2][2]+0.5);
delta -= boxVectors[1]*floor(delta[1]/boxVectors[1][1]+0.5);
delta -= boxVectors[0]*floor(delta[0]/boxVectors[0][0]+0.5);
return sqrt(delta.dot(delta));
}
double ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
int argIndex = -1;
for (int i = 0; i < 6; i++) {
if (derivOrder[i] > 0) {
if (derivOrder[i] > 1 || argIndex != -1)
throw OpenMMException("Unsupported derivative of periodicdistance"); // Should be impossible for this to happen.
argIndex = i;
}
}
Vec3* boxVectors = *boxVectorHandle;
Vec3 delta = Vec3(arguments[0], arguments[1], arguments[2])-Vec3(arguments[3], arguments[4], arguments[5]);
delta -= boxVectors[2]*floor(delta[2]/boxVectors[2][2]+0.5);
delta -= boxVectors[1]*floor(delta[1]/boxVectors[1][1]+0.5);
delta -= boxVectors[0]*floor(delta[0]/boxVectors[0][0]+0.5);
double r = sqrt(delta.dot(delta));
if (r == 0)
return 0.0;
if (argIndex < 3)
return delta[argIndex]/r;
return -delta[argIndex-3]/r;
}
Lepton::CustomFunction* ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::clone() const {
return new PeriodicDistanceFunction(boxVectorHandle);
}
ReferenceCalcCustomExternalForceKernel::~ReferenceCalcCustomExternalForceKernel() { ReferenceCalcCustomExternalForceKernel::~ReferenceCalcCustomExternalForceKernel() {
if (ixn != NULL) if (ixn != NULL)
delete ixn; delete ixn;
...@@ -1586,7 +1544,7 @@ void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, co ...@@ -1586,7 +1544,7 @@ void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, co
// Parse the expression used to calculate the force. // Parse the expression used to calculate the force.
map<string, Lepton::CustomFunction*> functions; map<string, Lepton::CustomFunction*> functions;
PeriodicDistanceFunction periodicDistance(&boxVectors); ReferencePointDistanceFunction periodicDistance(true, &boxVectors);
functions["periodicdistance"] = &periodicDistance; functions["periodicdistance"] = &periodicDistance;
Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize(); Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
energyExpression = expression.createCompiledExpression(); energyExpression = expression.createCompiledExpression();
...@@ -1790,12 +1748,15 @@ void ReferenceCalcCustomCentroidBondForceKernel::initialize(const System& system ...@@ -1790,12 +1748,15 @@ void ReferenceCalcCustomCentroidBondForceKernel::initialize(const System& system
for (int i = 0; i < force.getNumFunctions(); i++) for (int i = 0; i < force.getNumFunctions(); i++)
functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i)); functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
// Create implementations of point functions.
functions["pointdistance"] = new ReferencePointDistanceFunction(usePeriodic, &boxVectors);
functions["pointangle"] = new ReferencePointAngleFunction(usePeriodic, &boxVectors);
functions["pointdihedral"] = new ReferencePointDihedralFunction(usePeriodic, &boxVectors);
// Parse the expression and create the object used to calculate the interaction. // Parse the expression and create the object used to calculate the interaction.
map<string, vector<int> > distances; Lepton::ParsedExpression energyExpression = CustomCentroidBondForceImpl::prepareExpression(force, functions);
map<string, vector<int> > angles;
map<string, vector<int> > dihedrals;
Lepton::ParsedExpression energyExpression = CustomCentroidBondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
vector<string> bondParameterNames; vector<string> bondParameterNames;
for (int i = 0; i < numBondParameters; i++) for (int i = 0; i < numBondParameters; i++)
bondParameterNames.push_back(force.getPerBondParameterName(i)); bondParameterNames.push_back(force.getPerBondParameterName(i));
...@@ -1807,7 +1768,7 @@ void ReferenceCalcCustomCentroidBondForceKernel::initialize(const System& system ...@@ -1807,7 +1768,7 @@ void ReferenceCalcCustomCentroidBondForceKernel::initialize(const System& system
energyParamDerivNames.push_back(param); energyParamDerivNames.push_back(param);
energyParamDerivExpressions.push_back(energyExpression.differentiate(param).createCompiledExpression()); energyParamDerivExpressions.push_back(energyExpression.differentiate(param).createCompiledExpression());
} }
ixn = new ReferenceCustomCentroidBondIxn(force.getNumGroupsPerBond(), groupAtoms, normalizedWeights, bondGroups, energyExpression, bondParameterNames, distances, angles, dihedrals, energyParamDerivExpressions); ixn = new ReferenceCustomCentroidBondIxn(force.getNumGroupsPerBond(), groupAtoms, normalizedWeights, bondGroups, energyExpression, bondParameterNames, energyParamDerivExpressions);
// Delete the custom functions. // Delete the custom functions.
...@@ -1822,8 +1783,10 @@ double ReferenceCalcCustomCentroidBondForceKernel::execute(ContextImpl& context, ...@@ -1822,8 +1783,10 @@ double ReferenceCalcCustomCentroidBondForceKernel::execute(ContextImpl& context,
map<string, double> globalParameters; map<string, double> globalParameters;
for (auto& name : globalParameterNames) for (auto& name : globalParameterNames)
globalParameters[name] = context.getParameter(name); globalParameters[name] = context.getParameter(name);
if (usePeriodic) if (usePeriodic) {
ixn->setPeriodic(extractBoxVectors(context)); boxVectors = extractBoxVectors(context);
ixn->setPeriodic(boxVectors);
}
vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0); vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]); ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context); map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
...@@ -1875,12 +1838,15 @@ void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system ...@@ -1875,12 +1838,15 @@ void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system
for (int i = 0; i < force.getNumFunctions(); i++) for (int i = 0; i < force.getNumFunctions(); i++)
functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i)); functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
// Create implementations of point functions.
functions["pointdistance"] = new ReferencePointDistanceFunction(usePeriodic, &boxVectors);
functions["pointangle"] = new ReferencePointAngleFunction(usePeriodic, &boxVectors);
functions["pointdihedral"] = new ReferencePointDihedralFunction(usePeriodic, &boxVectors);
// Parse the expression and create the object used to calculate the interaction. // Parse the expression and create the object used to calculate the interaction.
map<string, vector<int> > distances; Lepton::ParsedExpression energyExpression = CustomCompoundBondForceImpl::prepareExpression(force, functions);
map<string, vector<int> > angles;
map<string, vector<int> > dihedrals;
Lepton::ParsedExpression energyExpression = CustomCompoundBondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
vector<string> bondParameterNames; vector<string> bondParameterNames;
for (int i = 0; i < numBondParameters; i++) for (int i = 0; i < numBondParameters; i++)
bondParameterNames.push_back(force.getPerBondParameterName(i)); bondParameterNames.push_back(force.getPerBondParameterName(i));
...@@ -1892,7 +1858,7 @@ void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system ...@@ -1892,7 +1858,7 @@ void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system
energyParamDerivNames.push_back(param); energyParamDerivNames.push_back(param);
energyParamDerivExpressions.push_back(energyExpression.differentiate(param).createCompiledExpression()); energyParamDerivExpressions.push_back(energyExpression.differentiate(param).createCompiledExpression());
} }
ixn = new ReferenceCustomCompoundBondIxn(force.getNumParticlesPerBond(), bondParticles, energyExpression, bondParameterNames, distances, angles, dihedrals, energyParamDerivExpressions); ixn = new ReferenceCustomCompoundBondIxn(force.getNumParticlesPerBond(), bondParticles, energyExpression, bondParameterNames, energyParamDerivExpressions);
// Delete the custom functions. // Delete the custom functions.
...@@ -1907,8 +1873,10 @@ double ReferenceCalcCustomCompoundBondForceKernel::execute(ContextImpl& context, ...@@ -1907,8 +1873,10 @@ double ReferenceCalcCustomCompoundBondForceKernel::execute(ContextImpl& context,
map<string, double> globalParameters; map<string, double> globalParameters;
for (auto& name : globalParameterNames) for (auto& name : globalParameterNames)
globalParameters[name] = context.getParameter(name); globalParameters[name] = context.getParameter(name);
if (usePeriodic) if (usePeriodic) {
ixn->setPeriodic(extractBoxVectors(context)); boxVectors = extractBoxVectors(context);
ixn->setPeriodic(boxVectors);
}
vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0); vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]); ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context); map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
......
platforms/reference/src/ReferencePointFunctions.cpp 0 → 100644
View file @ a85c2428
/* -------------------------------------------------------------------------- *
* 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) 2021 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 "ReferencePointFunctions.h"
#include "ReferenceBondIxn.h"
#include "openmm/OpenMMException.h"
#include <cmath>
using namespace OpenMM;
ReferencePointDistanceFunction::ReferencePointDistanceFunction(bool periodic, Vec3** boxVectorHandle) : periodic(periodic), boxVectorHandle(boxVectorHandle) {
}
int ReferencePointDistanceFunction::getNumArguments() const {
return 6;
}
double ReferencePointDistanceFunction::evaluate(const double* arguments) const {
Vec3 delta = Vec3(arguments[0], arguments[1], arguments[2])-Vec3(arguments[3], arguments[4], arguments[5]);
if (periodic) {
Vec3* boxVectors = *boxVectorHandle;
delta -= boxVectors[2]*floor(delta[2]/boxVectors[2][2]+0.5);
delta -= boxVectors[1]*floor(delta[1]/boxVectors[1][1]+0.5);
delta -= boxVectors[0]*floor(delta[0]/boxVectors[0][0]+0.5);
}
return sqrt(delta.dot(delta));
}
double ReferencePointDistanceFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
int argIndex = -1;
for (int i = 0; i < 6; i++) {
if (derivOrder[i] > 0) {
if (derivOrder[i] > 1 || argIndex != -1)
throw OpenMMException("Unsupported derivative of pointdistance"); // Should be impossible for this to happen.
argIndex = i;
}
}
Vec3 delta = Vec3(arguments[0], arguments[1], arguments[2])-Vec3(arguments[3], arguments[4], arguments[5]);
if (periodic) {
Vec3* boxVectors = *boxVectorHandle;
delta -= boxVectors[2]*floor(delta[2]/boxVectors[2][2]+0.5);
delta -= boxVectors[1]*floor(delta[1]/boxVectors[1][1]+0.5);
delta -= boxVectors[0]*floor(delta[0]/boxVectors[0][0]+0.5);
}
double r = sqrt(delta.dot(delta));
if (r == 0)
return 0.0;
if (argIndex < 3)
return delta[argIndex]/r;
return -delta[argIndex-3]/r;
}
Lepton::CustomFunction* ReferencePointDistanceFunction::clone() const {
return new ReferencePointDistanceFunction(periodic, boxVectorHandle);
}
ReferencePointAngleFunction::ReferencePointAngleFunction(bool periodic, Vec3** boxVectorHandle) : periodic(periodic), boxVectorHandle(boxVectorHandle) {
}
int ReferencePointAngleFunction::getNumArguments() const {
return 9;
}
double ReferencePointAngleFunction::evaluate(const double* arguments) const {
Vec3 delta12 = Vec3(arguments[3], arguments[4], arguments[5])-Vec3(arguments[0], arguments[1], arguments[2]);
Vec3 delta32 = Vec3(arguments[3], arguments[4], arguments[5])-Vec3(arguments[6], arguments[7], arguments[8]);
if (periodic) {
Vec3* boxVectors = *boxVectorHandle;
delta12 -= boxVectors[2]*floor(delta12[2]/boxVectors[2][2]+0.5);
delta12 -= boxVectors[1]*floor(delta12[1]/boxVectors[1][1]+0.5);
delta12 -= boxVectors[0]*floor(delta12[0]/boxVectors[0][0]+0.5);
delta32 -= boxVectors[2]*floor(delta32[2]/boxVectors[2][2]+0.5);
delta32 -= boxVectors[1]*floor(delta32[1]/boxVectors[1][1]+0.5);
delta32 -= boxVectors[0]*floor(delta32[0]/boxVectors[0][0]+0.5);
}
return ReferenceBondIxn::getAngleBetweenTwoVectors(&delta12[0], &delta32[0]);
}
double ReferencePointAngleFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
int argIndex = -1;
for (int i = 0; i < 9; i++) {
if (derivOrder[i] > 0) {
if (derivOrder[i] > 1 || argIndex != -1)
throw OpenMMException("Unsupported derivative of pointangle"); // Should be impossible for this to happen.
argIndex = i;
}
}
Vec3 delta12 = Vec3(arguments[3], arguments[4], arguments[5])-Vec3(arguments[0], arguments[1], arguments[2]);
Vec3 delta32 = Vec3(arguments[3], arguments[4], arguments[5])-Vec3(arguments[6], arguments[7], arguments[8]);
if (periodic) {
Vec3* boxVectors = *boxVectorHandle;
delta12 -= boxVectors[2]*floor(delta12[2]/boxVectors[2][2]+0.5);
delta12 -= boxVectors[1]*floor(delta12[1]/boxVectors[1][1]+0.5);
delta12 -= boxVectors[0]*floor(delta12[0]/boxVectors[0][0]+0.5);
delta32 -= boxVectors[2]*floor(delta32[2]/boxVectors[2][2]+0.5);
delta32 -= boxVectors[1]*floor(delta32[1]/boxVectors[1][1]+0.5);
delta32 -= boxVectors[0]*floor(delta32[0]/boxVectors[0][0]+0.5);
}
Vec3 thetaCross = delta12.cross(delta32);
double lengthThetaCross = sqrt(thetaCross.dot(thetaCross));
if (lengthThetaCross < 1.0e-6)
lengthThetaCross = 1.0e-6;
Vec3 deltaCrossP[3];
deltaCrossP[0] = delta12.cross(thetaCross)/(delta12.dot(delta12)*lengthThetaCross);
deltaCrossP[2] = -delta32.cross(thetaCross)/(delta32.dot(delta32)*lengthThetaCross);
deltaCrossP[1] = -(deltaCrossP[0]+deltaCrossP[2]);
return -deltaCrossP[argIndex/3][argIndex%3];
}
Lepton::CustomFunction* ReferencePointAngleFunction::clone() const {
return new ReferencePointAngleFunction(periodic, boxVectorHandle);
}
ReferencePointDihedralFunction::ReferencePointDihedralFunction(bool periodic, Vec3** boxVectorHandle) : periodic(periodic), boxVectorHandle(boxVectorHandle) {
}
int ReferencePointDihedralFunction::getNumArguments() const {
return 12;
}
double ReferencePointDihedralFunction::evaluate(const double* arguments) const {
Vec3 delta12 = Vec3(arguments[0], arguments[1], arguments[2])-Vec3(arguments[3], arguments[4], arguments[5]);
Vec3 delta32 = Vec3(arguments[6], arguments[7], arguments[8])-Vec3(arguments[3], arguments[4], arguments[5]);
Vec3 delta34 = Vec3(arguments[6], arguments[7], arguments[8])-Vec3(arguments[9], arguments[10], arguments[11]);
if (periodic) {
Vec3* boxVectors = *boxVectorHandle;
delta12 -= boxVectors[2]*floor(delta12[2]/boxVectors[2][2]+0.5);
delta12 -= boxVectors[1]*floor(delta12[1]/boxVectors[1][1]+0.5);
delta12 -= boxVectors[0]*floor(delta12[0]/boxVectors[0][0]+0.5);
delta32 -= boxVectors[2]*floor(delta32[2]/boxVectors[2][2]+0.5);
delta32 -= boxVectors[1]*floor(delta32[1]/boxVectors[1][1]+0.5);
delta32 -= boxVectors[0]*floor(delta32[0]/boxVectors[0][0]+0.5);
delta34 -= boxVectors[2]*floor(delta34[2]/boxVectors[2][2]+0.5);
delta34 -= boxVectors[1]*floor(delta34[1]/boxVectors[1][1]+0.5);
delta34 -= boxVectors[0]*floor(delta34[0]/boxVectors[0][0]+0.5);
}
return ReferenceBondIxn::getDihedralAngleBetweenThreeVectors(&delta12[0], &delta32[0], &delta34[0], NULL, NULL, &delta12[0]);
}
double ReferencePointDihedralFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
int argIndex = -1;
for (int i = 0; i < 12; i++) {
if (derivOrder[i] > 0) {
if (derivOrder[i] > 1 || argIndex != -1)
throw OpenMMException("Unsupported derivative of pointdihedral"); // Should be impossible for this to happen.
argIndex = i;
}
}
Vec3 delta12 = Vec3(arguments[0], arguments[1], arguments[2])-Vec3(arguments[3], arguments[4], arguments[5]);
Vec3 delta32 = Vec3(arguments[6], arguments[7], arguments[8])-Vec3(arguments[3], arguments[4], arguments[5]);
Vec3 delta34 = Vec3(arguments[6], arguments[7], arguments[8])-Vec3(arguments[9], arguments[10], arguments[11]);
if (periodic) {
Vec3* boxVectors = *boxVectorHandle;
delta12 -= boxVectors[2]*floor(delta12[2]/boxVectors[2][2]+0.5);
delta12 -= boxVectors[1]*floor(delta12[1]/boxVectors[1][1]+0.5);
delta12 -= boxVectors[0]*floor(delta12[0]/boxVectors[0][0]+0.5);
delta32 -= boxVectors[2]*floor(delta32[2]/boxVectors[2][2]+0.5);
delta32 -= boxVectors[1]*floor(delta32[1]/boxVectors[1][1]+0.5);
delta32 -= boxVectors[0]*floor(delta32[0]/boxVectors[0][0]+0.5);
delta34 -= boxVectors[2]*floor(delta34[2]/boxVectors[2][2]+0.5);
delta34 -= boxVectors[1]*floor(delta34[1]/boxVectors[1][1]+0.5);
delta34 -= boxVectors[0]*floor(delta34[0]/boxVectors[0][0]+0.5);
}
Vec3 cross1 = delta12.cross(delta32);
Vec3 cross2 = delta32.cross(delta34);
double norm2Cross1 = cross1.dot(cross1);
double norm2Cross2 = cross2.dot(cross2);
double norm2Delta32 = delta32.dot(delta32);
double normDelta32 = sqrt(norm2Delta32);
double forceFactors[4];
forceFactors[0] = -normDelta32/norm2Cross1;
forceFactors[3] = normDelta32/norm2Cross2;
forceFactors[1] = delta12.dot(delta32);
forceFactors[1] /= norm2Delta32;
forceFactors[2] = delta34.dot(delta32);
forceFactors[2] /= norm2Delta32;
Vec3 internalF[4];
internalF[0] = forceFactors[0]*cross1;
internalF[3] = forceFactors[3]*cross2;
Vec3 s = forceFactors[1]*internalF[0] - forceFactors[2]*internalF[3];
internalF[1] = internalF[0] - s;
internalF[2] = internalF[3] + s;
internalF[0] *= -1;
internalF[3] *= -1;
return internalF[argIndex/3][argIndex%3];
}
Lepton::CustomFunction* ReferencePointDihedralFunction::clone() const {
return new ReferencePointDihedralFunction(periodic, boxVectorHandle);
}
platforms/reference/src/SimTKReference/ReferenceBondIxn.cpp
View file @ a85c2428
...@@ -83,7 +83,7 @@ ReferenceBondIxn::~ReferenceBondIxn() { ...@@ -83,7 +83,7 @@ ReferenceBondIxn::~ReferenceBondIxn() {
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
double ReferenceBondIxn::getNormedDotProduct(double* vector1, double* vector2, double ReferenceBondIxn::getNormedDotProduct(double* vector1, double* vector2,
int hasREntry = 0) { int hasREntry) {
double dotProduct = DOT3(vector1, vector2); double dotProduct = DOT3(vector1, vector2);
if (dotProduct != 0.0) { if (dotProduct != 0.0) {
...@@ -122,9 +122,7 @@ double ReferenceBondIxn::getNormedDotProduct(double* vector1, double* vector2, ...@@ -122,9 +122,7 @@ double ReferenceBondIxn::getNormedDotProduct(double* vector1, double* vector2,
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
double ReferenceBondIxn::getAngleBetweenTwoVectors(double* vector1, double* vector2, double ReferenceBondIxn::getAngleBetweenTwoVectors(double* vector1, double* vector2, double* outputDotProduct, int hasREntry) {
double* outputDotProduct = NULL,
int hasREntry = 0) {
// get dot product betweenn vectors and then angle // get dot product betweenn vectors and then angle
...@@ -173,11 +171,11 @@ double ReferenceBondIxn::getAngleBetweenTwoVectors(double* vector1, double* vect ...@@ -173,11 +171,11 @@ double ReferenceBondIxn::getAngleBetweenTwoVectors(double* vector1, double* vect
double ReferenceBondIxn::getDihedralAngleBetweenThreeVectors(double* vector1, double ReferenceBondIxn::getDihedralAngleBetweenThreeVectors(double* vector1,
double* vector2, double* vector2,
double* vector3, double* vector3,
double** outputCrossProduct = NULL, double** outputCrossProduct,
double* cosineOfAngle = NULL, double* cosineOfAngle,
double* signVector = NULL, double* signVector,
double* signOfAngle = NULL, double* signOfAngle,
int hasREntry = 0) { int hasREntry) {
double tempVectors[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; double tempVectors[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
......
platforms/reference/src/SimTKReference/ReferenceCustomCentroidBondIxn.cpp
View file @ a85c2428
...@@ -40,7 +40,6 @@ using namespace OpenMM; ...@@ -40,7 +40,6 @@ using namespace OpenMM;
ReferenceCustomCentroidBondIxn::ReferenceCustomCentroidBondIxn(int numGroupsPerBond, const vector<vector<int> >& groupAtoms, ReferenceCustomCentroidBondIxn::ReferenceCustomCentroidBondIxn(int numGroupsPerBond, const vector<vector<int> >& groupAtoms,
const vector<vector<double> >& normalizedWeights, const vector<vector<int> >& bondGroups, const vector<vector<double> >& normalizedWeights, const vector<vector<int> >& bondGroups,
const Lepton::ParsedExpression& energyExpression, const vector<string>& bondParameterNames, const Lepton::ParsedExpression& energyExpression, const vector<string>& bondParameterNames,
const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals,
const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions) : const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions) :
groupAtoms(groupAtoms), normalizedWeights(normalizedWeights), bondGroups(bondGroups), energyExpression(energyExpression.createCompiledExpression()), groupAtoms(groupAtoms), normalizedWeights(normalizedWeights), bondGroups(bondGroups), energyExpression(energyExpression.createCompiledExpression()),
usePeriodic(false), energyParamDerivExpressions(energyParamDerivExpressions) { usePeriodic(false), energyParamDerivExpressions(energyParamDerivExpressions) {
...@@ -56,28 +55,10 @@ ReferenceCustomCentroidBondIxn::ReferenceCustomCentroidBondIxn(int numGroupsPerB ...@@ -56,28 +55,10 @@ ReferenceCustomCentroidBondIxn::ReferenceCustomCentroidBondIxn(int numGroupsPerB
positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).createCompiledExpression())); positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).createCompiledExpression()));
positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).createCompiledExpression())); positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).createCompiledExpression()));
} }
for (auto& term : distances)
distanceTerms.push_back(ReferenceCustomCentroidBondIxn::DistanceTermInfo(term.first, term.second, energyExpression.differentiate(term.first).createCompiledExpression()));
for (auto& term : angles)
angleTerms.push_back(ReferenceCustomCentroidBondIxn::AngleTermInfo(term.first, term.second, energyExpression.differentiate(term.first).createCompiledExpression()));
for (auto& term : dihedrals)
dihedralTerms.push_back(ReferenceCustomCentroidBondIxn::DihedralTermInfo(term.first, term.second, energyExpression.differentiate(term.first).createCompiledExpression()));
for (int i = 0; i < positionTerms.size(); i++) { for (int i = 0; i < positionTerms.size(); i++) {
expressionSet.registerExpression(positionTerms[i].forceExpression); expressionSet.registerExpression(positionTerms[i].forceExpression);
positionTerms[i].index = expressionSet.getVariableIndex(positionTerms[i].name); positionTerms[i].index = expressionSet.getVariableIndex(positionTerms[i].name);
} }
for (int i = 0; i < distanceTerms.size(); i++) {
expressionSet.registerExpression(distanceTerms[i].forceExpression);
distanceTerms[i].index = expressionSet.getVariableIndex(distanceTerms[i].name);
}
for (int i = 0; i < angleTerms.size(); i++) {
expressionSet.registerExpression(angleTerms[i].forceExpression);
angleTerms[i].index = expressionSet.getVariableIndex(angleTerms[i].name);
}
for (int i = 0; i < dihedralTerms.size(); i++) {
expressionSet.registerExpression(dihedralTerms[i].forceExpression);
dihedralTerms[i].index = expressionSet.getVariableIndex(dihedralTerms[i].name);
}
numParameters = bondParameterNames.size(); numParameters = bondParameterNames.size();
for (int i = 0; i < numParameters; i++) for (int i = 0; i < numParameters; i++)
bondParamIndex.push_back(expressionSet.getVariableIndex(bondParameterNames[i])); bondParamIndex.push_back(expressionSet.getVariableIndex(bondParameterNames[i]));
...@@ -133,96 +114,12 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<Vec3>& gro ...@@ -133,96 +114,12 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<Vec3>& gro
const vector<int>& groups = bondGroups[bond]; const vector<int>& groups = bondGroups[bond];
for (auto& term : positionTerms) for (auto& term : positionTerms)
expressionSet.setVariable(term.index, groupCenters[groups[term.group]][term.component]); expressionSet.setVariable(term.index, groupCenters[groups[term.group]][term.component]);
for (auto& term : distanceTerms) {
computeDelta(groups[term.g1], groups[term.g2], term.delta, groupCenters);
expressionSet.setVariable(term.index, term.delta[ReferenceForce::RIndex]);
}
for (auto& term : angleTerms) {
computeDelta(groups[term.g1], groups[term.g2], term.delta1, groupCenters);
computeDelta(groups[term.g3], groups[term.g2], term.delta2, groupCenters);
expressionSet.setVariable(term.index, computeAngle(term.delta1, term.delta2));
}
for (auto& term : dihedralTerms) {
computeDelta(groups[term.g2], groups[term.g1], term.delta1, groupCenters);
computeDelta(groups[term.g2], groups[term.g3], term.delta2, groupCenters);
computeDelta(groups[term.g4], groups[term.g3], term.delta3, groupCenters);
double dotDihedral, signOfDihedral;
double* crossProduct[] = {term.cross1, term.cross2};
expressionSet.setVariable(term.index, getDihedralAngleBetweenThreeVectors(term.delta1, term.delta2, term.delta3, crossProduct, &dotDihedral, term.delta1, &signOfDihedral, 1));
}
// Apply forces based on individual particle coordinates. // Apply forces based on particle coordinates.
for (auto& term : positionTerms) for (auto& term : positionTerms)
forces[groups[term.group]][term.component] -= term.forceExpression.evaluate(); forces[groups[term.group]][term.component] -= term.forceExpression.evaluate();
// Apply forces based on distances.
for (auto& term : distanceTerms) {
double dEdR = term.forceExpression.evaluate()/(term.delta[ReferenceForce::RIndex]);
for (int i = 0; i < 3; i++) {
double force = -dEdR*term.delta[i];
forces[groups[term.g1]][i] -= force;
forces[groups[term.g2]][i] += force;
}
}
// Apply forces based on angles.
for (auto& term : angleTerms) {
double dEdTheta = term.forceExpression.evaluate();
double thetaCross[ReferenceForce::LastDeltaRIndex];
SimTKOpenMMUtilities::crossProductVector3(term.delta1, term.delta2, thetaCross);
double lengthThetaCross = sqrt(DOT3(thetaCross, thetaCross));
if (lengthThetaCross < 1.0e-06)
lengthThetaCross = 1.0e-06;
double termA = dEdTheta/(term.delta1[ReferenceForce::R2Index]*lengthThetaCross);
double termC = -dEdTheta/(term.delta2[ReferenceForce::R2Index]*lengthThetaCross);
double deltaCrossP[3][3];
SimTKOpenMMUtilities::crossProductVector3(term.delta1, thetaCross, deltaCrossP[0]);
SimTKOpenMMUtilities::crossProductVector3(term.delta2, thetaCross, deltaCrossP[2]);
for (int i = 0; i < 3; i++) {
deltaCrossP[0][i] *= termA;
deltaCrossP[2][i] *= termC;
deltaCrossP[1][i] = -(deltaCrossP[0][i]+deltaCrossP[2][i]);
}
for (int i = 0; i < 3; i++) {
forces[groups[term.g1]][i] += deltaCrossP[0][i];
forces[groups[term.g2]][i] += deltaCrossP[1][i];
forces[groups[term.g3]][i] += deltaCrossP[2][i];
}
}
// Apply forces based on dihedrals.
for (auto& term : dihedralTerms) {
double dEdTheta = term.forceExpression.evaluate();
double internalF[4][3];
double forceFactors[4];
double normCross1 = DOT3(term.cross1, term.cross1);
double normBC = term.delta2[ReferenceForce::RIndex];
forceFactors[0] = (-dEdTheta*normBC)/normCross1;
double normCross2 = DOT3(term.cross2, term.cross2);
forceFactors[3] = (dEdTheta*normBC)/normCross2;
forceFactors[1] = DOT3(term.delta1, term.delta2);
forceFactors[1] /= term.delta2[ReferenceForce::R2Index];
forceFactors[2] = DOT3(term.delta3, term.delta2);
forceFactors[2] /= term.delta2[ReferenceForce::R2Index];
for (int i = 0; i < 3; i++) {
internalF[0][i] = forceFactors[0]*term.cross1[i];
internalF[3][i] = forceFactors[3]*term.cross2[i];
double s = forceFactors[1]*internalF[0][i] - forceFactors[2]*internalF[3][i];
internalF[1][i] = internalF[0][i] - s;
internalF[2][i] = internalF[3][i] + s;
}
for (int i = 0; i < 3; i++) {
forces[groups[term.g1]][i] += internalF[0][i];
forces[groups[term.g2]][i] -= internalF[1][i];
forces[groups[term.g3]][i] -= internalF[2][i];
forces[groups[term.g4]][i] += internalF[3][i];
}
}
// Add the energy // Add the energy
if (totalEnergy) if (totalEnergy)
......
platforms/reference/src/SimTKReference/ReferenceCustomCompoundBondIxn.cpp
View file @ a85c2428
...@@ -45,7 +45,6 @@ using namespace OpenMM; ...@@ -45,7 +45,6 @@ using namespace OpenMM;
ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesPerBond, const vector<vector<int> >& bondAtoms, ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesPerBond, const vector<vector<int> >& bondAtoms,
const Lepton::ParsedExpression& energyExpression, const vector<string>& bondParameterNames, const Lepton::ParsedExpression& energyExpression, const vector<string>& bondParameterNames,
const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals,
const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions) : const std::vector<Lepton::CompiledExpression> energyParamDerivExpressions) :
bondAtoms(bondAtoms), energyExpression(energyExpression.createCompiledExpression()), usePeriodic(false), bondAtoms(bondAtoms), energyExpression(energyExpression.createCompiledExpression()), usePeriodic(false),
energyParamDerivExpressions(energyParamDerivExpressions) { energyParamDerivExpressions(energyParamDerivExpressions) {
...@@ -61,28 +60,10 @@ ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesP ...@@ -61,28 +60,10 @@ ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesP
particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).createCompiledExpression())); particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).createCompiledExpression()));
particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).createCompiledExpression())); particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).createCompiledExpression()));
} }
for (auto& term : distances)
distanceTerms.push_back(ReferenceCustomCompoundBondIxn::DistanceTermInfo(term.first, term.second, energyExpression.differentiate(term.first).createCompiledExpression()));
for (auto& term : angles)
angleTerms.push_back(ReferenceCustomCompoundBondIxn::AngleTermInfo(term.first, term.second, energyExpression.differentiate(term.first).createCompiledExpression()));
for (auto& term : dihedrals)
dihedralTerms.push_back(ReferenceCustomCompoundBondIxn::DihedralTermInfo(term.first, term.second, energyExpression.differentiate(term.first).createCompiledExpression()));
for (int i = 0; i < particleTerms.size(); i++) { for (int i = 0; i < particleTerms.size(); i++) {
expressionSet.registerExpression(particleTerms[i].forceExpression); expressionSet.registerExpression(particleTerms[i].forceExpression);
particleTerms[i].index = expressionSet.getVariableIndex(particleTerms[i].name); particleTerms[i].index = expressionSet.getVariableIndex(particleTerms[i].name);
} }
for (int i = 0; i < distanceTerms.size(); i++) {
expressionSet.registerExpression(distanceTerms[i].forceExpression);
distanceTerms[i].index = expressionSet.getVariableIndex(distanceTerms[i].name);
}
for (int i = 0; i < angleTerms.size(); i++) {
expressionSet.registerExpression(angleTerms[i].forceExpression);
angleTerms[i].index = expressionSet.getVariableIndex(angleTerms[i].name);
}
for (int i = 0; i < dihedralTerms.size(); i++) {
expressionSet.registerExpression(dihedralTerms[i].forceExpression);
dihedralTerms[i].index = expressionSet.getVariableIndex(dihedralTerms[i].name);
}
numParameters = bondParameterNames.size(); numParameters = bondParameterNames.size();
for (int i = 0; i < numParameters; i++) for (int i = 0; i < numParameters; i++)
bondParamIndex.push_back(expressionSet.getVariableIndex(bondParameterNames[i])); bondParamIndex.push_back(expressionSet.getVariableIndex(bondParameterNames[i]));
...@@ -148,96 +129,12 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<Vec3>& ato ...@@ -148,96 +129,12 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<Vec3>& ato
const vector<int>& atoms = bondAtoms[bond]; const vector<int>& atoms = bondAtoms[bond];
for (auto& term : particleTerms) for (auto& term : particleTerms)
expressionSet.setVariable(term.index, atomCoordinates[atoms[term.atom]][term.component]); expressionSet.setVariable(term.index, atomCoordinates[atoms[term.atom]][term.component]);
for (auto& term : distanceTerms) {
computeDelta(atoms[term.p1], atoms[term.p2], term.delta, atomCoordinates);
expressionSet.setVariable(term.index, term.delta[ReferenceForce::RIndex]);
}
for (auto& term : angleTerms) {
computeDelta(atoms[term.p1], atoms[term.p2], term.delta1, atomCoordinates);
computeDelta(atoms[term.p3], atoms[term.p2], term.delta2, atomCoordinates);
expressionSet.setVariable(term.index, computeAngle(term.delta1, term.delta2));
}
for (auto& term : dihedralTerms) {
computeDelta(atoms[term.p2], atoms[term.p1], term.delta1, atomCoordinates);
computeDelta(atoms[term.p2], atoms[term.p3], term.delta2, atomCoordinates);
computeDelta(atoms[term.p4], atoms[term.p3], term.delta3, atomCoordinates);
double dotDihedral, signOfDihedral;
double* crossProduct[] = {term.cross1, term.cross2};
expressionSet.setVariable(term.index,getDihedralAngleBetweenThreeVectors(term.delta1, term.delta2, term.delta3, crossProduct, &dotDihedral, term.delta1, &signOfDihedral, 1));
}
// Apply forces based on individual particle coordinates. // Apply forces based on particle coordinates.
for (auto& term : particleTerms) for (auto& term : particleTerms)
forces[atoms[term.atom]][term.component] -= term.forceExpression.evaluate(); forces[atoms[term.atom]][term.component] -= term.forceExpression.evaluate();
// Apply forces based on distances.
for (auto& term : distanceTerms) {
double dEdR = term.forceExpression.evaluate()/(term.delta[ReferenceForce::RIndex]);
for (int i = 0; i < 3; i++) {
double force = -dEdR*term.delta[i];
forces[atoms[term.p1]][i] -= force;
forces[atoms[term.p2]][i] += force;
}
}
// Apply forces based on angles.
for (auto& term : angleTerms) {
double dEdTheta = term.forceExpression.evaluate();
double thetaCross[ReferenceForce::LastDeltaRIndex];
SimTKOpenMMUtilities::crossProductVector3(term.delta1, term.delta2, thetaCross);
double lengthThetaCross = sqrt(DOT3(thetaCross, thetaCross));
if (lengthThetaCross < 1.0e-06)
lengthThetaCross = 1.0e-06;
double termA = dEdTheta/(term.delta1[ReferenceForce::R2Index]*lengthThetaCross);
double termC = -dEdTheta/(term.delta2[ReferenceForce::R2Index]*lengthThetaCross);
double deltaCrossP[3][3];
SimTKOpenMMUtilities::crossProductVector3(term.delta1, thetaCross, deltaCrossP[0]);
SimTKOpenMMUtilities::crossProductVector3(term.delta2, thetaCross, deltaCrossP[2]);
for (int i = 0; i < 3; i++) {
deltaCrossP[0][i] *= termA;
deltaCrossP[2][i] *= termC;
deltaCrossP[1][i] = -(deltaCrossP[0][i]+deltaCrossP[2][i]);
}
for (int i = 0; i < 3; i++) {
forces[atoms[term.p1]][i] += deltaCrossP[0][i];
forces[atoms[term.p2]][i] += deltaCrossP[1][i];
forces[atoms[term.p3]][i] += deltaCrossP[2][i];
}
}
// Apply forces based on dihedrals.
for (auto& term : dihedralTerms) {
double dEdTheta = term.forceExpression.evaluate();
double internalF[4][3];
double forceFactors[4];
double normCross1 = DOT3(term.cross1, term.cross1);
double normBC = term.delta2[ReferenceForce::RIndex];
forceFactors[0] = (-dEdTheta*normBC)/normCross1;
double normCross2 = DOT3(term.cross2, term.cross2);
forceFactors[3] = (dEdTheta*normBC)/normCross2;
forceFactors[1] = DOT3(term.delta1, term.delta2);
forceFactors[1] /= term.delta2[ReferenceForce::R2Index];
forceFactors[2] = DOT3(term.delta3, term.delta2);
forceFactors[2] /= term.delta2[ReferenceForce::R2Index];
for (int i = 0; i < 3; i++) {
internalF[0][i] = forceFactors[0]*term.cross1[i];
internalF[3][i] = forceFactors[3]*term.cross2[i];
double s = forceFactors[1]*internalF[0][i] - forceFactors[2]*internalF[3][i];
internalF[1][i] = internalF[0][i] - s;
internalF[2][i] = internalF[3][i] + s;
}
for (int i = 0; i < 3; i++) {
forces[atoms[term.p1]][i] += internalF[0][i];
forces[atoms[term.p2]][i] -= internalF[1][i];
forces[atoms[term.p3]][i] -= internalF[2][i];
forces[atoms[term.p4]][i] += internalF[3][i];
}
}
// Add the energy // Add the energy
if (totalEnergy) if (totalEnergy)
......
platforms/reference/src/SimTKReference/ReferenceCustomManyParticleIxn.cpp
View file @ a85c2428
/* Portions copyright (c) 2009-2021 Stanford University and Simbios.
/* Portions copyright (c) 2009-2018 Stanford University and Simbios.
* Contributors: Peter Eastman * Contributors: Peter Eastman
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -29,6 +28,7 @@ ...@@ -29,6 +28,7 @@
#include "SimTKOpenMMUtilities.h" #include "SimTKOpenMMUtilities.h"
#include "ReferenceForce.h" #include "ReferenceForce.h"
#include "ReferenceCustomManyParticleIxn.h" #include "ReferenceCustomManyParticleIxn.h"
#include "ReferencePointFunctions.h"
#include "ReferenceTabulatedFunction.h" #include "ReferenceTabulatedFunction.h"
#include "openmm/internal/CustomManyParticleForceImpl.h" #include "openmm/internal/CustomManyParticleForceImpl.h"
#include "lepton/CustomFunction.h" #include "lepton/CustomFunction.h"
...@@ -47,14 +47,16 @@ ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(const CustomManyP ...@@ -47,14 +47,16 @@ ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(const CustomManyP
for (int i = 0; i < (int) force.getNumTabulatedFunctions(); i++) for (int i = 0; i < (int) force.getNumTabulatedFunctions(); i++)
functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i)); functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
// Create implementations of point functions.
functions["pointdistance"] = new ReferencePointDistanceFunction(force.usesPeriodicBoundaryConditions(), &periodicBoxVectors);
functions["pointangle"] = new ReferencePointAngleFunction(force.usesPeriodicBoundaryConditions(), &periodicBoxVectors);
functions["pointdihedral"] = new ReferencePointDihedralFunction(force.usesPeriodicBoundaryConditions(), &periodicBoxVectors);
// Parse the expression and create the object used to calculate the interaction. // Parse the expression and create the object used to calculate the interaction.
map<string, vector<int> > distances; Lepton::ParsedExpression energyExpr = CustomManyParticleForceImpl::prepareExpression(force, functions);
map<string, vector<int> > angles;
map<string, vector<int> > dihedrals;
Lepton::ParsedExpression energyExpr = CustomManyParticleForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
energyExpression = energyExpr.createProgram(); energyExpression = energyExpr.createProgram();
vector<string> particleParameterNames;
if (force.getNonbondedMethod() != CustomManyParticleForce::NoCutoff) if (force.getNonbondedMethod() != CustomManyParticleForce::NoCutoff)
setUseCutoff(force.getCutoffDistance()); setUseCutoff(force.getCutoffDistance());
...@@ -80,12 +82,6 @@ ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(const CustomManyP ...@@ -80,12 +82,6 @@ ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(const CustomManyP
particleParamNames[i].push_back(paramname.str()); particleParamNames[i].push_back(paramname.str());
} }
} }
for (auto& term : distances)
distanceTerms.push_back(ReferenceCustomManyParticleIxn::DistanceTermInfo(term.first, term.second, energyExpr.differentiate(term.first).optimize().createProgram()));
for (auto& term : angles)
angleTerms.push_back(ReferenceCustomManyParticleIxn::AngleTermInfo(term.first, term.second, energyExpr.differentiate(term.first).optimize().createProgram()));
for (auto& term : dihedrals)
dihedralTerms.push_back(ReferenceCustomManyParticleIxn::DihedralTermInfo(term.first, term.second, energyExpr.differentiate(term.first).optimize().createProgram()));
// Record exclusions. // Record exclusions.
...@@ -124,9 +120,7 @@ void ReferenceCustomManyParticleIxn::setPeriodic(Vec3* vectors) { ...@@ -124,9 +120,7 @@ void ReferenceCustomManyParticleIxn::setPeriodic(Vec3* vectors) {
assert(vectors[1][1] >= 2.0*cutoffDistance); assert(vectors[1][1] >= 2.0*cutoffDistance);
assert(vectors[2][2] >= 2.0*cutoffDistance); assert(vectors[2][2] >= 2.0*cutoffDistance);
usePeriodic = true; usePeriodic = true;
periodicBoxVectors[0] = vectors[0]; periodicBoxVectors = vectors;
periodicBoxVectors[1] = vectors[1];
periodicBoxVectors[2] = vectors[2];
} }
void ReferenceCustomManyParticleIxn::loopOverInteractions(vector<int>& particles, int loopIndex, vector<OpenMM::Vec3>& atomCoordinates, void ReferenceCustomManyParticleIxn::loopOverInteractions(vector<int>& particles, int loopIndex, vector<OpenMM::Vec3>& atomCoordinates,
...@@ -190,100 +184,16 @@ void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particle ...@@ -190,100 +184,16 @@ void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particle
for (int j = 0; j < numPerParticleParameters; j++) for (int j = 0; j < numPerParticleParameters; j++)
variables[particleParamNames[i][j]] = particleParameters[permutedParticles[i]][j]; variables[particleParamNames[i][j]] = particleParameters[permutedParticles[i]][j];
// Compute all of the variables the energy can depend on. // Record particle coordinates.
for (auto& term : particleTerms) for (auto& term : particleTerms)
variables[term.name] = atomCoordinates[permutedParticles[term.atom]][term.component]; variables[term.name] = atomCoordinates[permutedParticles[term.atom]][term.component];
for (auto& term : distanceTerms) {
computeDelta(permutedParticles[term.p1], permutedParticles[term.p2], term.delta, atomCoordinates);
variables[term.name] = term.delta[ReferenceForce::RIndex];
}
for (auto& term : angleTerms) {
computeDelta(permutedParticles[term.p1], permutedParticles[term.p2], term.delta1, atomCoordinates);
computeDelta(permutedParticles[term.p3], permutedParticles[term.p2], term.delta2, atomCoordinates);
variables[term.name] = computeAngle(term.delta1, term.delta2);
}
for (auto& term : dihedralTerms) {
computeDelta(permutedParticles[term.p2], permutedParticles[term.p1], term.delta1, atomCoordinates);
computeDelta(permutedParticles[term.p2], permutedParticles[term.p3], term.delta2, atomCoordinates);
computeDelta(permutedParticles[term.p4], permutedParticles[term.p3], term.delta3, atomCoordinates);
double dotDihedral, signOfDihedral;
double* crossProduct[] = {term.cross1, term.cross2};
variables[term.name] = ReferenceBondIxn::getDihedralAngleBetweenThreeVectors(term.delta1, term.delta2, term.delta3, crossProduct, &dotDihedral, term.delta1, &signOfDihedral, 1);
}
// Apply forces based on individual particle coordinates. // Apply forces based on particle coordinates.
for (auto& term : particleTerms) for (auto& term : particleTerms)
forces[permutedParticles[term.atom]][term.component] -= term.forceExpression.evaluate(variables); forces[permutedParticles[term.atom]][term.component] -= term.forceExpression.evaluate(variables);
// Apply forces based on distances.
for (auto& term : distanceTerms) {
double dEdR = term.forceExpression.evaluate(variables)/(term.delta[ReferenceForce::RIndex]);
for (int i = 0; i < 3; i++) {
double force = -dEdR*term.delta[i];
forces[permutedParticles[term.p1]][i] -= force;
forces[permutedParticles[term.p2]][i] += force;
}
}
// Apply forces based on angles.
for (auto& term : angleTerms) {
double dEdTheta = term.forceExpression.evaluate(variables);
double thetaCross[ReferenceForce::LastDeltaRIndex];
SimTKOpenMMUtilities::crossProductVector3(term.delta1, term.delta2, thetaCross);
double lengthThetaCross = sqrt(DOT3(thetaCross, thetaCross));
if (lengthThetaCross < 1.0e-06)
lengthThetaCross = 1.0e-06;
double termA = dEdTheta/(term.delta1[ReferenceForce::R2Index]*lengthThetaCross);
double termC = -dEdTheta/(term.delta2[ReferenceForce::R2Index]*lengthThetaCross);
double deltaCrossP[3][3];
SimTKOpenMMUtilities::crossProductVector3(term.delta1, thetaCross, deltaCrossP[0]);
SimTKOpenMMUtilities::crossProductVector3(term.delta2, thetaCross, deltaCrossP[2]);
for (int i = 0; i < 3; i++) {
deltaCrossP[0][i] *= termA;
deltaCrossP[2][i] *= termC;
deltaCrossP[1][i] = -(deltaCrossP[0][i]+deltaCrossP[2][i]);
}
for (int i = 0; i < 3; i++) {
forces[permutedParticles[term.p1]][i] += deltaCrossP[0][i];
forces[permutedParticles[term.p2]][i] += deltaCrossP[1][i];
forces[permutedParticles[term.p3]][i] += deltaCrossP[2][i];
}
}
// Apply forces based on dihedrals.
for (auto& term : dihedralTerms) {
double dEdTheta = term.forceExpression.evaluate(variables);
double internalF[4][3];
double forceFactors[4];
double normCross1 = DOT3(term.cross1, term.cross1);
double normBC = term.delta2[ReferenceForce::RIndex];
forceFactors[0] = (-dEdTheta*normBC)/normCross1;
double normCross2 = DOT3(term.cross2, term.cross2);
forceFactors[3] = (dEdTheta*normBC)/normCross2;
forceFactors[1] = DOT3(term.delta1, term.delta2);
forceFactors[1] /= term.delta2[ReferenceForce::R2Index];
forceFactors[2] = DOT3(term.delta3, term.delta2);
forceFactors[2] /= term.delta2[ReferenceForce::R2Index];
for (int i = 0; i < 3; i++) {
internalF[0][i] = forceFactors[0]*term.cross1[i];
internalF[3][i] = forceFactors[3]*term.cross2[i];
double s = forceFactors[1]*internalF[0][i] - forceFactors[2]*internalF[3][i];
internalF[1][i] = internalF[0][i] - s;
internalF[2][i] = internalF[3][i] + s;
}
for (int i = 0; i < 3; i++) {
forces[permutedParticles[term.p1]][i] += internalF[0][i];
forces[permutedParticles[term.p2]][i] -= internalF[1][i];
forces[permutedParticles[term.p3]][i] -= internalF[2][i];
forces[permutedParticles[term.p4]][i] += internalF[3][i];
}
}
// Add the energy // Add the energy
if (totalEnergy) if (totalEnergy)
......
tests/TestCustomCentroidBondForce.h
View file @ a85c2428
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2015-2016 Stanford University and the Authors. * * Portions copyright (c) 2015-2021 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -117,7 +117,7 @@ void testHarmonicBond() { ...@@ -117,7 +117,7 @@ void testHarmonicBond() {
ASSERT_EQUAL(5, molecules[0].size()); ASSERT_EQUAL(5, molecules[0].size());
} }
void testComplexFunction() { void testComplexFunction(bool byGroups) {
int numParticles = 5; int numParticles = 5;
System system; System system;
for (int i = 0; i < numParticles; i++) for (int i = 0; i < numParticles; i++)
...@@ -130,7 +130,12 @@ void testComplexFunction() { ...@@ -130,7 +130,12 @@ void testComplexFunction() {
// CustomCompoundBondForce. Use that to test a complicated energy function with lots of terms. // CustomCompoundBondForce. Use that to test a complicated energy function with lots of terms.
CustomCompoundBondForce* compound = new CustomCompoundBondForce(4, "x1+y2+z4+fn(distance(p1,p2))*angle(p3,p2,p4)+scale*dihedral(p2,p1,p4,p3)"); CustomCompoundBondForce* compound = new CustomCompoundBondForce(4, "x1+y2+z4+fn(distance(p1,p2))*angle(p3,p2,p4)+scale*dihedral(p2,p1,p4,p3)");
CustomCentroidBondForce* centroid = new CustomCentroidBondForce(4, "x1+y2+z4+fn(distance(g1,g2))*angle(g3,g2,g4)+scale*dihedral(g2,g1,g4,g3)"); string expression;
if (byGroups)
expression = "x1+y2+z4+fn(distance(g1,g2))*angle(g3,g2,g4)+scale*dihedral(g2,g1,g4,g3)";
else
expression = "x1+y2+z4+fn(pointdistance(x1,y1,z1,x2,y2,z2))*pointangle(x3,y3,z3,x2,y2,z2,x4,y4,z4)+scale*pointdihedral(x2,y2,z2,x1,y1,z1,x4,y4,z4,x3,y3,z3)";
CustomCentroidBondForce* centroid = new CustomCentroidBondForce(4, expression);
compound->addGlobalParameter("scale", 0.5); compound->addGlobalParameter("scale", 0.5);
centroid->addGlobalParameter("scale", 0.5); centroid->addGlobalParameter("scale", 0.5);
compound->addTabulatedFunction("fn", new Continuous1DFunction(table, -1, 10)); compound->addTabulatedFunction("fn", new Continuous1DFunction(table, -1, 10));
...@@ -277,7 +282,7 @@ void testIllegalVariable() { ...@@ -277,7 +282,7 @@ void testIllegalVariable() {
ASSERT(threwException); ASSERT(threwException);
} }
void testPeriodic() { void testPeriodic(bool byGroups) {
// Create a force that uses periodic boundary conditions. // Create a force that uses periodic boundary conditions.
System system; System system;
...@@ -287,7 +292,12 @@ void testPeriodic() { ...@@ -287,7 +292,12 @@ void testPeriodic() {
system.addParticle(4.0); system.addParticle(4.0);
system.addParticle(5.0); system.addParticle(5.0);
system.setDefaultPeriodicBoxVectors(Vec3(2, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3)); system.setDefaultPeriodicBoxVectors(Vec3(2, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "k*distance(g1,g2)^2"); string expression;
if (byGroups)
expression = "k*distance(g1,g2)^2";
else
expression = "k*pointdistance(x1,y1,z1,x2,y2,z2)^2";
CustomCentroidBondForce* force = new CustomCentroidBondForce(2, expression);
force->addPerBondParameter("k"); force->addPerBondParameter("k");
vector<int> particles1; vector<int> particles1;
particles1.push_back(0); particles1.push_back(0);
...@@ -398,10 +408,12 @@ int main(int argc, char* argv[]) { ...@@ -398,10 +408,12 @@ int main(int argc, char* argv[]) {
try { try {
initializeTests(argc, argv); initializeTests(argc, argv);
testHarmonicBond(); testHarmonicBond();
testComplexFunction(); testComplexFunction(true);
testComplexFunction(false);
testCustomWeights(); testCustomWeights();
testIllegalVariable(); testIllegalVariable();
testPeriodic(); testPeriodic(true);
testPeriodic(false);
testEnergyParameterDerivatives(); testEnergyParameterDerivatives();
runPlatformTests(); runPlatformTests();
} }
......
tests/TestCustomCompoundBondForce.h
View file @ a85c2428
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2012-2016 Stanford University and the Authors. * * Portions copyright (c) 2012-2021 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -50,7 +50,7 @@ using namespace std; ...@@ -50,7 +50,7 @@ using namespace std;
const double TOL = 1e-5; const double TOL = 1e-5;
void testBond() { void testBond(bool byParticles) {
// Create a system using a CustomCompoundBondForce. // Create a system using a CustomCompoundBondForce.
System customSystem; System customSystem;
...@@ -58,7 +58,12 @@ void testBond() { ...@@ -58,7 +58,12 @@ void testBond() {
customSystem.addParticle(1.0); customSystem.addParticle(1.0);
customSystem.addParticle(1.0); customSystem.addParticle(1.0);
customSystem.addParticle(1.0); customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(4, "0.5*kb*((distance(p1,p2)-b0)^2+(distance(p2,p3)-b0)^2)+0.5*ka*(angle(p2,p3,p4)-a0)^2+kt*(1+cos(dihedral(p1,p2,p3,p4)-t0))"); string expression;
if (byParticles)
expression = "0.5*kb*((distance(p1,p2)-b0)^2+(distance(p2,p3)-b0)^2)+0.5*ka*(angle(p2,p3,p4)-a0)^2+kt*(1+cos(dihedral(p1,p2,p3,p4)-t0))";
else
expression = "0.5*kb*((pointdistance(x1,y1,z1,x2,y2,z2)-b0)^2+(pointdistance(x2,y2,z2,x3,y3,z3)-b0)^2)+0.5*ka*(pointangle(x2,y2,z2,x3,y3,z3,x4,y4,z4)-a0)^2+kt*(1+cos(pointdihedral(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4)-t0))";
CustomCompoundBondForce* custom = new CustomCompoundBondForce(4, expression);
custom->addPerBondParameter("kb"); custom->addPerBondParameter("kb");
custom->addPerBondParameter("ka"); custom->addPerBondParameter("ka");
custom->addPerBondParameter("kt"); custom->addPerBondParameter("kt");
...@@ -330,7 +335,7 @@ void testIllegalVariable() { ...@@ -330,7 +335,7 @@ void testIllegalVariable() {
ASSERT(threwException); ASSERT(threwException);
} }
void testPeriodic() { void testPeriodic(bool byParticles) {
// Create a force that uses periodic boundary conditions. // Create a force that uses periodic boundary conditions.
System customSystem; System customSystem;
...@@ -339,7 +344,12 @@ void testPeriodic() { ...@@ -339,7 +344,12 @@ void testPeriodic() {
customSystem.addParticle(1.0); customSystem.addParticle(1.0);
customSystem.addParticle(1.0); customSystem.addParticle(1.0);
customSystem.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3)); customSystem.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
CustomCompoundBondForce* custom = new CustomCompoundBondForce(4, "0.5*kb*((distance(p1,p2)-b0)^2+(distance(p2,p3)-b0)^2)+0.5*ka*(angle(p2,p3,p4)-a0)^2+kt*(1+cos(dihedral(p1,p2,p3,p4)-t0))"); string expression;
if (byParticles)
expression = "0.5*kb*((distance(p1,p2)-b0)^2+(distance(p2,p3)-b0)^2)+0.5*ka*(angle(p2,p3,p4)-a0)^2+kt*(1+cos(dihedral(p1,p2,p3,p4)-t0))";
else
expression = "0.5*kb*((pointdistance(x1,y1,z1,x2,y2,z2)-b0)^2+(pointdistance(x2,y2,z2,x3,y3,z3)-b0)^2)+0.5*ka*(pointangle(x2,y2,z2,x3,y3,z3,x4,y4,z4)-a0)^2+kt*(1+cos(pointdihedral(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4)-t0))";
CustomCompoundBondForce* custom = new CustomCompoundBondForce(4, expression);
custom->addPerBondParameter("kb"); custom->addPerBondParameter("kb");
custom->addPerBondParameter("ka"); custom->addPerBondParameter("ka");
custom->addPerBondParameter("kt"); custom->addPerBondParameter("kt");
...@@ -453,13 +463,15 @@ void runPlatformTests(); ...@@ -453,13 +463,15 @@ void runPlatformTests();
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
try { try {
initializeTests(argc, argv); initializeTests(argc, argv);
testBond(); testBond(true);
testBond(false);
testPositionDependence(); testPositionDependence();
testContinuous2DFunction(); testContinuous2DFunction();
testContinuous3DFunction(); testContinuous3DFunction();
testMultipleBonds(); testMultipleBonds();
testIllegalVariable(); testIllegalVariable();
testPeriodic(); testPeriodic(true);
testPeriodic(false);
testEnergyParameterDerivatives(); testEnergyParameterDerivatives();
runPlatformTests(); runPlatformTests();
} }
......
tests/TestCustomManyParticleForce.h
View file @ a85c2428
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2014-2015 Stanford University and the Authors. * * Portions copyright (c) 2014-2021 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -48,6 +48,17 @@ using namespace std; ...@@ -48,6 +48,17 @@ using namespace std;
const double TOL = 1e-5; const double TOL = 1e-5;
string axilrodTellerParticleExpression =
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)";
string axilrodTellerPointExpression =
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=pointangle(x1,y1,z1,x2,y2,z2,x3,y3,z3); theta2=pointangle(x2,y2,z2,x3,y3,z3,x1,y1,z1); theta3=pointangle(x3,y3,z3,x1,y1,z1,x2,y2,z2);"
"r12=pointdistance(x1,y1,z1,x2,y2,z2); r13=pointdistance(x1,y1,z1,x3,y3,z3); r23=pointdistance(x2,y2,z2,x3,y3,z3)";
Vec3 computeDelta(const Vec3& pos1, const Vec3& pos2, bool periodic, const Vec3* periodicBoxVectors) { Vec3 computeDelta(const Vec3& pos1, const Vec3& pos2, bool periodic, const Vec3* periodicBoxVectors) {
Vec3 diff = pos1-pos2; Vec3 diff = pos1-pos2;
if (periodic) { if (periodic) {
...@@ -208,11 +219,8 @@ void validateStillingerWeber(CustomManyParticleForce* force, const vector<Vec3>& ...@@ -208,11 +219,8 @@ void validateStillingerWeber(CustomManyParticleForce* force, const vector<Vec3>&
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4); ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
} }
void testNoCutoff() { void testNoCutoff(bool byParticle) {
CustomManyParticleForce* force = new CustomManyParticleForce(3, CustomManyParticleForce* force = new CustomManyParticleForce(3, byParticle ? axilrodTellerParticleExpression : axilrodTellerPointExpression);
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5); force->addGlobalParameter("C", 1.5);
vector<double> params; vector<double> params;
force->addParticle(params); force->addParticle(params);
...@@ -229,11 +237,8 @@ void testNoCutoff() { ...@@ -229,11 +237,8 @@ void testNoCutoff() {
validateAxilrodTeller(force, positions, expectedSets, 2.0, false); validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
} }
void testCutoff() { void testCutoff(bool byParticle) {
CustomManyParticleForce* force = new CustomManyParticleForce(3, CustomManyParticleForce* force = new CustomManyParticleForce(3, byParticle ? axilrodTellerParticleExpression : axilrodTellerPointExpression);
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5); force->addGlobalParameter("C", 1.5);
force->setNonbondedMethod(CustomManyParticleForce::CutoffNonPeriodic); force->setNonbondedMethod(CustomManyParticleForce::CutoffNonPeriodic);
force->setCutoffDistance(1.55); force->setCutoffDistance(1.55);
...@@ -254,11 +259,8 @@ void testCutoff() { ...@@ -254,11 +259,8 @@ void testCutoff() {
validateAxilrodTeller(force, positions, expectedSets, 2.0, false); validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
} }
void testPeriodic() { void testPeriodic(bool byParticle) {
CustomManyParticleForce* force = new CustomManyParticleForce(3, CustomManyParticleForce* force = new CustomManyParticleForce(3, byParticle ? axilrodTellerParticleExpression : axilrodTellerPointExpression);
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5); force->addGlobalParameter("C", 1.5);
force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic); force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
force->setCutoffDistance(1.05); force->setCutoffDistance(1.05);
...@@ -280,11 +282,8 @@ void testPeriodic() { ...@@ -280,11 +282,8 @@ void testPeriodic() {
validateAxilrodTeller(force, positions, expectedSets, boxSize, false); validateAxilrodTeller(force, positions, expectedSets, boxSize, false);
} }
void testTriclinic() { void testTriclinic(bool byParticle) {
CustomManyParticleForce* force = new CustomManyParticleForce(3, CustomManyParticleForce* force = new CustomManyParticleForce(3, byParticle ? axilrodTellerParticleExpression : axilrodTellerPointExpression);
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5); force->addGlobalParameter("C", 1.5);
force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic); force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
force->setCutoffDistance(1.05); force->setCutoffDistance(1.05);
...@@ -307,10 +306,7 @@ void testTriclinic() { ...@@ -307,10 +306,7 @@ void testTriclinic() {
} }
void testExclusions() { void testExclusions() {
CustomManyParticleForce* force = new CustomManyParticleForce(3, CustomManyParticleForce* force = new CustomManyParticleForce(3, axilrodTellerParticleExpression);
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5); force->addGlobalParameter("C", 1.5);
vector<double> params; vector<double> params;
force->addParticle(params); force->addParticle(params);
...@@ -578,10 +574,7 @@ void testLargeSystem() { ...@@ -578,10 +574,7 @@ void testLargeSystem() {
int numParticles = gridSize*gridSize*gridSize; int numParticles = gridSize*gridSize*gridSize;
double boxSize = 3.0; double boxSize = 3.0;
double spacing = boxSize/gridSize; double spacing = boxSize/gridSize;
CustomManyParticleForce* force = new CustomManyParticleForce(3, CustomManyParticleForce* force = new CustomManyParticleForce(3, axilrodTellerParticleExpression);
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5); force->addGlobalParameter("C", 1.5);
force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic); force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
force->setCutoffDistance(0.6); force->setCutoffDistance(0.6);
...@@ -729,10 +722,14 @@ void runPlatformTests(); ...@@ -729,10 +722,14 @@ void runPlatformTests();
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
try { try {
initializeTests(argc, argv); initializeTests(argc, argv);
testNoCutoff(); testNoCutoff(true);
testCutoff(); testNoCutoff(false);
testPeriodic(); testCutoff(true);
testTriclinic(); testCutoff(false);
testPeriodic(true);
testPeriodic(false);
testTriclinic(true);
testTriclinic(false);
testExclusions(); testExclusions();
testAllTerms(); testAllTerms();
testParameters(); testParameters();
......
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