Skip to content

GitLab

Commit 56d95193 authored by peastman's avatar peastman
Browse files

Created reference implementation of CustomCentroidBondForce

parent f1603534
Show whitespace changes
Inline Side-by-side
Showing with 861 additions and 10 deletions
openmmapi/include/openmm/CustomCentroidBondForce.h
View file @ 56d95193
......@@ -320,7 +320,8 @@ public:
*
* This method has several limitations. The only information it updates is the values of per-bond parameters.
* All other aspects of the Force (such as the energy function) are unaffected and can only be changed by reinitializing
* the Context. The set of groups involved in a bond cannot be changed, nor can new bonds be added.
* the Context. Neither the definitions of groups nor the set of groups involved in a bond can be changed, nor can new
* bonds be added.
*/
void updateParametersInContext(Context& context);
/**
......
openmmapi/include/openmm/internal/CustomCentroidBondForceImpl.h
View file @ 56d95193
......@@ -35,6 +35,7 @@
#include "ForceImpl.h"
#include "openmm/CustomCentroidBondForce.h"
#include "openmm/Kernel.h"
#include "openmm/System.h"
#include "lepton/CustomFunction.h"
#include "lepton/ExpressionTreeNode.h"
#include "lepton/ParsedExpression.h"
......@@ -62,6 +63,7 @@ public:
double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
std::map<std::string, double> getDefaultParameters();
std::vector<std::string> getKernelNames();
std::vector<std::pair<int, int> > getBondedParticles() const;
void updateParametersInContext(ContextImpl& context);
/**
* This is a utility routine that parses the energy expression, identifies the angles and dihedrals
......@@ -79,11 +81,20 @@ public:
*/
static Lepton::ParsedExpression prepareExpression(const CustomCentroidBondForce& force, const std::map<std::string, Lepton::CustomFunction*>& functions, std::map<std::string, std::vector<int> >& distances,
std::map<std::string, std::vector<int> >& angles, std::map<std::string, std::vector<int> >& dihedrals);
/**
* Compute the normalized weights to use for each particle in each group.
*
* @param force the CustomCentroidBondForce to process
* @param system the System it is part of
* @param weights on exit, weights[i][j] contains the normalized weight for particle j in group i.
*/
static void computeNormalizedWeights(const CustomCentroidBondForce& force, const System& system, std::vector<std::vector<double> >& weights);
private:
class FunctionPlaceholder;
static Lepton::ExpressionTreeNode replaceFunctions(const Lepton::ExpressionTreeNode& node, std::map<std::string, int> atoms,
std::map<std::string, std::vector<int> >& distances, std::map<std::string, std::vector<int> >& angles,
std::map<std::string, std::vector<int> >& dihedrals);
void addBondsBetweenGroups(int group1, int group2, std::vector<std::pair<int, int> >& bonds) const;
const CustomCentroidBondForce& owner;
Kernel kernel;
};
......
openmmapi/src/CustomCentroidBondForceImpl.cpp
View file @ 56d95193
......@@ -36,18 +36,14 @@
#include "lepton/Operation.h"
#include "lepton/Parser.h"
#include <sstream>
#include <utility>
using namespace OpenMM;
using namespace std;
using Lepton::CustomFunction;
using Lepton::ExpressionTreeNode;
using Lepton::Operation;
using Lepton::ParsedExpression;
using std::map;
using std::pair;
using std::vector;
using std::set;
using std::string;
using std::stringstream;
/**
* This class serves as a placeholder for angles and dihedrals in expressions.
......@@ -211,6 +207,65 @@ ExpressionTreeNode CustomCentroidBondForceImpl::replaceFunctions(const Expressio
return ExpressionTreeNode(new Operation::Variable(name));
}
vector<pair<int, int> > CustomCentroidBondForceImpl::getBondedParticles() const {
vector<pair<int, int> > bonds;
for (int i = 0; i < owner.getNumBonds(); i++) {
vector<int> groups;
vector<double> parameters;
owner.getBondParameters(i, groups, parameters);
for (int j = 1; j < groups.size(); j++)
for (int k = 0; k < j; k++)
addBondsBetweenGroups(j, k, bonds);
}
return bonds;
}
void CustomCentroidBondForceImpl::addBondsBetweenGroups(int group1, int group2, vector<pair<int, int> >& bonds) const {
vector<int> atoms1;
vector<int> atoms2;
vector<double> weights;
owner.getGroupParameters(group1, atoms1, weights);
owner.getGroupParameters(group2, atoms2, weights);
for (int i = 0; i < atoms1.size(); i++)
for (int j = 0; j < atoms2.size(); j++)
bonds.push_back(make_pair(atoms1[i], atoms2[j]));
}
void CustomCentroidBondForceImpl::updateParametersInContext(ContextImpl& context) {
kernel.getAs<CalcCustomCentroidBondForceKernel>().copyParametersToContext(context, owner);
}
void CustomCentroidBondForceImpl::computeNormalizedWeights(const CustomCentroidBondForce& force, const System& system, vector<vector<double> >& weights) {
int numGroups = force.getNumGroups();
weights.resize(numGroups);
for (int i = 0; i < numGroups; i++) {
vector<int> particles;
vector<double> groupWeights;
force.getGroupParameters(i, particles, groupWeights);
int numParticles = particles.size();
// If weights were not specified, use particle masses.
if (groupWeights.size() == 0) {
groupWeights.resize(numParticles);
for (int j = 0; j < numParticles; j++)
groupWeights[j] = system.getParticleMass(particles[j]);
}
// Normalize the weights.
double total = 0;
for (int j = 0; j < numParticles; j++)
total += groupWeights[j];
if (total == 0.0) {
stringstream msg;
msg << "CustomCentroidBondForce: Weights for group ";
msg << i;
msg << " add to 0";
throw OpenMMException(msg.str());
}
weights[i].resize(numParticles);
for (int j = 0; j < numParticles; j++)
weights[i][j] = groupWeights[j]/total;
}
}
platforms/reference/include/ReferenceCustomCentroidBondIxn.h 0 → 100644
View file @ 56d95193
/* Portions copyright (c) 2009-2015 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __ReferenceCustomCentroidBondIxn_H__
#define __ReferenceCustomCentroidBondIxn_H__
#include "ReferenceBondIxn.h"
#include "lepton/ExpressionProgram.h"
#include "lepton/ParsedExpression.h"
#include <map>
#include <vector>
namespace OpenMM {
class ReferenceCustomCentroidBondIxn : public ReferenceBondIxn {
private:
class PositionTermInfo;
class DistanceTermInfo;
class AngleTermInfo;
class DihedralTermInfo;
std::vector<std::vector<int> > groupAtoms;
std::vector<std::vector<double> > normalizedWeights;
std::vector<std::vector<int> > bondGroups;
Lepton::ExpressionProgram energyExpression;
std::vector<std::string> bondParamNames;
std::vector<PositionTermInfo> positionTerms;
std::vector<DistanceTermInfo> distanceTerms;
std::vector<AngleTermInfo> angleTerms;
std::vector<DihedralTermInfo> dihedralTerms;
/**---------------------------------------------------------------------------------------
Calculate custom interaction for one bond
@param bond the index of the bond
@param groupCenters group center coordinates
@param variables the values of variables that may appear in expressions
@param forces force array (forces added)
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculateOneIxn(int bond, std::vector<OpenMM::RealVec>& groupCenters,
std::map<std::string, double>& variables, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* totalEnergy) const;
void computeDelta(int group1, int group2, RealOpenMM* delta, std::vector<OpenMM::RealVec>& groupCenters) const;
static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2);
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
ReferenceCustomCentroidBondIxn(int numGroupsPerBond, const std::vector<std::vector<int> >& groupAtoms,
const std::vector<std::vector<double> >& normalizedWeights, const std::vector<std::vector<int> >& bondGroups, const Lepton::ParsedExpression& energyExpression,
const std::vector<std::string>& bondParameterNames, const std::map<std::string, std::vector<int> >& distances,
const std::map<std::string, std::vector<int> >& angles, const std::map<std::string, std::vector<int> >& dihedrals);
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceCustomCentroidBondIxn();
/**---------------------------------------------------------------------------------------
Get the list of groups in each bond.
--------------------------------------------------------------------------------------- */
const std::vector<std::vector<int> >& getBondGroups() const {
return bondGroups;
}
/**---------------------------------------------------------------------------------------
Calculate custom compound bond interaction
@param atomCoordinates atom coordinates
@param bondParameters bond parameters values bondParameters[bondIndex][parameterIndex]
@param globalParameters the values of global parameters
@param forces force array (forces added)
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculatePairIxn(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** bondParameters,
const std::map<std::string, double>& globalParameters,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy) const;
// ---------------------------------------------------------------------------------------
};
class ReferenceCustomCentroidBondIxn::PositionTermInfo {
public:
std::string name;
int group, component;
Lepton::ExpressionProgram forceExpression;
PositionTermInfo(const std::string& name, int group, int component, const Lepton::ExpressionProgram& forceExpression) :
name(name), group(group), component(component), forceExpression(forceExpression) {
}
};
class ReferenceCustomCentroidBondIxn::DistanceTermInfo {
public:
std::string name;
int g1, g2;
Lepton::ExpressionProgram forceExpression;
mutable RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
DistanceTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::ExpressionProgram& forceExpression) :
name(name), g1(groups[0]), g2(groups[1]), forceExpression(forceExpression) {
}
};
class ReferenceCustomCentroidBondIxn::AngleTermInfo {
public:
std::string name;
int g1, g2, g3;
Lepton::ExpressionProgram forceExpression;
mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
AngleTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::ExpressionProgram& forceExpression) :
name(name), g1(groups[0]), g2(groups[1]), g3(groups[2]), forceExpression(forceExpression) {
}
};
class ReferenceCustomCentroidBondIxn::DihedralTermInfo {
public:
std::string name;
int g1, g2, g3, g4;
Lepton::ExpressionProgram forceExpression;
mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
mutable RealOpenMM delta3[ReferenceForce::LastDeltaRIndex];
mutable RealOpenMM cross1[3];
mutable RealOpenMM cross2[3];
DihedralTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::ExpressionProgram& forceExpression) :
name(name), g1(groups[0]), g2(groups[1]), g3(groups[2]), g4(groups[3]), forceExpression(forceExpression) {
}
};
} // namespace OpenMM
#endif // __ReferenceCustomCentroidBondIxn_H__
platforms/reference/include/ReferenceKernels.h
View file @ 56d95193
......@@ -44,6 +44,7 @@ namespace OpenMM {
class ReferenceObc;
class ReferenceGBVI;
class ReferenceAndersenThermostat;
class ReferenceCustomCentroidBondIxn;
class ReferenceCustomCompoundBondIxn;
class ReferenceCustomHbondIxn;
class ReferenceCustomManyParticleIxn;
......@@ -833,6 +834,44 @@ private:
std::vector<std::string> globalParameterNames;
};
/**
* This kernel is invoked by CustomCentroidBondForce to calculate the forces acting on the system.
*/
class ReferenceCalcCustomCentroidBondForceKernel : public CalcCustomCentroidBondForceKernel {
public:
ReferenceCalcCustomCentroidBondForceKernel(std::string name, const Platform& platform) : CalcCustomCentroidBondForceKernel(name, platform), ixn(NULL) {
}
~ReferenceCalcCustomCentroidBondForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the CustomCentroidBondForce this kernel will be used for
*/
void initialize(const System& system, const CustomCentroidBondForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the CustomCentroidBondForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force);
private:
int numBonds, numParticles;
RealOpenMM **bondParamArray;
ReferenceCustomCentroidBondIxn* ixn;
std::vector<std::string> globalParameterNames;
};
/**
* This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system.
*/
......@@ -865,7 +904,7 @@ public:
*/
void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force);
private:
int numBonds, numParticles;
int numBonds;
RealOpenMM **bondParamArray;
ReferenceCustomCompoundBondIxn* ixn;
std::vector<std::string> globalParameterNames;
......
platforms/reference/src/ReferenceKernelFactory.cpp
View file @ 56d95193
......@@ -78,6 +78,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
return new ReferenceCalcCustomExternalForceKernel(name, platform);
if (name == CalcCustomHbondForceKernel::Name())
return new ReferenceCalcCustomHbondForceKernel(name, platform);
if (name == CalcCustomCentroidBondForceKernel::Name())
return new ReferenceCalcCustomCentroidBondForceKernel(name, platform);
if (name == CalcCustomCompoundBondForceKernel::Name())
return new ReferenceCalcCustomCompoundBondForceKernel(name, platform);
if (name == CalcCustomManyParticleForceKernel::Name())
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 56d95193
......@@ -41,6 +41,7 @@
#include "ReferenceConstraints.h"
#include "ReferenceCustomAngleIxn.h"
#include "ReferenceCustomBondIxn.h"
#include "ReferenceCustomCentroidBondIxn.h"
#include "ReferenceCustomCompoundBondIxn.h"
#include "ReferenceCustomDynamics.h"
#include "ReferenceCustomExternalIxn.h"
......@@ -66,6 +67,7 @@
#include "openmm/System.h"
#include "openmm/internal/AndersenThermostatImpl.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/CustomCentroidBondForceImpl.h"
#include "openmm/internal/CustomCompoundBondForceImpl.h"
#include "openmm/internal/CustomHbondForceImpl.h"
#include "openmm/internal/CustomNonbondedForceImpl.h"
......@@ -1582,6 +1584,90 @@ void ReferenceCalcCustomHbondForceKernel::copyParametersToContext(ContextImpl& c
}
}
ReferenceCalcCustomCentroidBondForceKernel::~ReferenceCalcCustomCentroidBondForceKernel() {
disposeRealArray(bondParamArray, numBonds);
if (ixn != NULL)
delete ixn;
}
void ReferenceCalcCustomCentroidBondForceKernel::initialize(const System& system, const CustomCentroidBondForce& force) {
// Build the arrays.
int numGroups = force.getNumGroups();
vector<vector<int> > groupAtoms(numGroups);
vector<double> ignored;
for (int i = 0; i < numGroups; i++)
force.getGroupParameters(i, groupAtoms[i], ignored);
vector<vector<double> > normalizedWeights;
CustomCentroidBondForceImpl::computeNormalizedWeights(force, system, normalizedWeights);
numBonds = force.getNumBonds();
vector<vector<int> > bondGroups(numBonds);
int numBondParameters = force.getNumPerBondParameters();
bondParamArray = allocateRealArray(numBonds, numBondParameters);
for (int i = 0; i < numBonds; ++i) {
vector<double> parameters;
force.getBondParameters(i, bondGroups[i], parameters);
for (int j = 0; j < numBondParameters; j++)
bondParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
}
// Create custom functions for the tabulated functions.
map<string, Lepton::CustomFunction*> functions;
for (int i = 0; i < force.getNumFunctions(); i++)
functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
// Parse the expression and create the object used to calculate the interaction.
map<string, vector<int> > distances;
map<string, vector<int> > angles;
map<string, vector<int> > dihedrals;
Lepton::ParsedExpression energyExpression = CustomCentroidBondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
vector<string> bondParameterNames;
for (int i = 0; i < numBondParameters; i++)
bondParameterNames.push_back(force.getPerBondParameterName(i));
for (int i = 0; i < force.getNumGlobalParameters(); i++)
globalParameterNames.push_back(force.getGlobalParameterName(i));
ixn = new ReferenceCustomCentroidBondIxn(force.getNumGroupsPerBond(), groupAtoms, normalizedWeights, bondGroups, energyExpression, bondParameterNames, distances, angles, dihedrals);
// Delete the custom functions.
for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
delete iter->second;
}
double ReferenceCalcCustomCentroidBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
RealOpenMM energy = 0;
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
return energy;
}
void ReferenceCalcCustomCentroidBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force) {
if (numBonds != force.getNumBonds())
throw OpenMMException("updateParametersInContext: The number of bonds has changed");
// Record the values.
int numParameters = force.getNumPerBondParameters();
const vector<vector<int> >& bondGroups = ixn->getBondGroups();
vector<int> groups;
vector<double> params;
for (int i = 0; i < numBonds; ++i) {
force.getBondParameters(i, groups, params);
for (int j = 0; j < groups.size(); j++)
if (groups[j] != bondGroups[i][j])
throw OpenMMException("updateParametersInContext: The set of groups in a bond has changed");
for (int j = 0; j < numParameters; j++)
bondParamArray[i][j] = (RealOpenMM) params[j];
}
}
ReferenceCalcCustomCompoundBondForceKernel::~ReferenceCalcCustomCompoundBondForceKernel() {
disposeRealArray(bondParamArray, numBonds);
if (ixn != NULL)
......@@ -1593,7 +1679,6 @@ void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system
// Build the arrays.
numBonds = force.getNumBonds();
numParticles = system.getNumParticles();
vector<vector<int> > bondParticles(numBonds);
int numBondParameters = force.getNumPerBondParameters();
bondParamArray = allocateRealArray(numBonds, numBondParameters);
......@@ -1652,7 +1737,7 @@ void ReferenceCalcCustomCompoundBondForceKernel::copyParametersToContext(Context
vector<double> params;
for (int i = 0; i < numBonds; ++i) {
force.getBondParameters(i, particles, params);
for (int j = 0; j < numParticles; j++)
for (int j = 0; j < particles.size(); j++)
if (particles[j] != bondAtoms[i][j])
throw OpenMMException("updateParametersInContext: The set of particles in a bond has changed");
for (int j = 0; j < numParameters; j++)
......
platforms/reference/src/ReferencePlatform.cpp
View file @ 56d95193
......@@ -62,6 +62,7 @@ ReferencePlatform::ReferencePlatform() {
registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
registerKernelFactory(CalcCustomCentroidBondForceKernel::Name(), factory);
registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
......
platforms/reference/src/SimTKReference/ReferenceCustomCentroidBondIxn.cpp 0 → 100644
View file @ 56d95193
/* Portions copyright (c) 2009-2015 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include <string.h>
#include <sstream>
#include <utility>
#include "SimTKOpenMMUtilities.h"
#include "ReferenceForce.h"
#include "ReferenceCustomCentroidBondIxn.h"
using std::map;
using std::pair;
using std::string;
using std::stringstream;
using std::vector;
using namespace OpenMM;
ReferenceCustomCentroidBondIxn::ReferenceCustomCentroidBondIxn(int numGroupsPerBond, const vector<vector<int> >& groupAtoms,
const vector<vector<double> >& normalizedWeights, const vector<vector<int> >& bondGroups,
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) :
groupAtoms(groupAtoms), normalizedWeights(normalizedWeights), bondGroups(bondGroups), energyExpression(energyExpression.createProgram()), bondParamNames(bondParameterNames) {
for (int i = 0; i < numGroupsPerBond; i++) {
stringstream xname, yname, zname;
xname << 'x' << (i+1);
yname << 'y' << (i+1);
zname << 'z' << (i+1);
positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(xname.str(), i, 0, energyExpression.differentiate(xname.str()).optimize().createProgram()));
positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).optimize().createProgram()));
positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).optimize().createProgram()));
}
for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter)
distanceTerms.push_back(ReferenceCustomCentroidBondIxn::DistanceTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter)
angleTerms.push_back(ReferenceCustomCentroidBondIxn::AngleTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter)
dihedralTerms.push_back(ReferenceCustomCentroidBondIxn::DihedralTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
}
ReferenceCustomCentroidBondIxn::~ReferenceCustomCentroidBondIxn() {
}
void ReferenceCustomCentroidBondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates, RealOpenMM** bondParameters,
const map<string, double>& globalParameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy) const {
// First compute the center of each group.
int numGroups = groupAtoms.size();
vector<RealVec> groupCenters(numGroups);
for (int group = 0; group < numGroups; group++) {
for (int i = 0; i < groupAtoms[group].size(); i++)
groupCenters[group] += atomCoordinates[groupAtoms[group][i]]*normalizedWeights[group][i];
}
// Compute the forces on groups.
map<string, double> variables = globalParameters;
vector<RealVec> groupForces(numGroups);
int numBonds = bondGroups.size();
for (int bond = 0; bond < numBonds; bond++) {
for (int j = 0; j < (int) bondParamNames.size(); j++)
variables[bondParamNames[j]] = bondParameters[bond][j];
calculateOneIxn(bond, groupCenters, variables, groupForces, totalEnergy);
}
// Apply the forces to the individual atoms.
for (int group = 0; group < numGroups; group++) {
for (int i = 0; i < groupAtoms[group].size(); i++)
forces[groupAtoms[group][i]] += groupForces[group]*normalizedWeights[group][i];
}
}
void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<RealVec>& groupCenters,
map<string, double>& variables, vector<RealVec>& forces, RealOpenMM* totalEnergy) const {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\nReferenceCustomCentroidBondIxn::calculateOneIxn";
// ---------------------------------------------------------------------------------------
// Compute all of the variables the energy can depend on.
const vector<int>& groups = bondGroups[bond];
for (int i = 0; i < (int) positionTerms.size(); i++) {
const PositionTermInfo& term = positionTerms[i];
variables[term.name] = groupCenters[groups[term.group]][term.component];
}
for (int i = 0; i < (int) distanceTerms.size(); i++) {
const DistanceTermInfo& term = distanceTerms[i];
computeDelta(groups[term.g1], groups[term.g2], term.delta, groupCenters);
variables[term.name] = term.delta[ReferenceForce::RIndex];
}
for (int i = 0; i < (int) angleTerms.size(); i++) {
const AngleTermInfo& term = angleTerms[i];
computeDelta(groups[term.g1], groups[term.g2], term.delta1, groupCenters);
computeDelta(groups[term.g3], groups[term.g2], term.delta2, groupCenters);
variables[term.name] = computeAngle(term.delta1, term.delta2);
}
for (int i = 0; i < (int) dihedralTerms.size(); i++) {
const DihedralTermInfo& term = dihedralTerms[i];
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);
RealOpenMM dotDihedral, signOfDihedral;
RealOpenMM* crossProduct[] = {term.cross1, term.cross2};
variables[term.name] = getDihedralAngleBetweenThreeVectors(term.delta1, term.delta2, term.delta3, crossProduct, &dotDihedral, term.delta1, &signOfDihedral, 1);
}
// Apply forces based on individual particle coordinates.
for (int i = 0; i < (int) positionTerms.size(); i++) {
const PositionTermInfo& term = positionTerms[i];
forces[groups[term.group]][term.component] -= term.forceExpression.evaluate(variables);
}
// Apply forces based on distances.
for (int i = 0; i < (int) distanceTerms.size(); i++) {
const DistanceTermInfo& term = distanceTerms[i];
RealOpenMM dEdR = (RealOpenMM) (term.forceExpression.evaluate(variables)/(term.delta[ReferenceForce::RIndex]));
for (int i = 0; i < 3; i++) {
RealOpenMM force = -dEdR*term.delta[i];
forces[groups[term.g1]][i] -= force;
forces[groups[term.g2]][i] += force;
}
}
// Apply forces based on angles.
for (int i = 0; i < (int) angleTerms.size(); i++) {
const AngleTermInfo& term = angleTerms[i];
RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate(variables);
RealOpenMM thetaCross[ReferenceForce::LastDeltaRIndex];
SimTKOpenMMUtilities::crossProductVector3(term.delta1, term.delta2, thetaCross);
RealOpenMM lengthThetaCross = SQRT(DOT3(thetaCross, thetaCross));
if (lengthThetaCross < 1.0e-06)
lengthThetaCross = (RealOpenMM) 1.0e-06;
RealOpenMM termA = dEdTheta/(term.delta1[ReferenceForce::R2Index]*lengthThetaCross);
RealOpenMM termC = -dEdTheta/(term.delta2[ReferenceForce::R2Index]*lengthThetaCross);
RealOpenMM 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 (int i = 0; i < (int) dihedralTerms.size(); i++) {
const DihedralTermInfo& term = dihedralTerms[i];
RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate(variables);
RealOpenMM internalF[4][3];
RealOpenMM forceFactors[4];
RealOpenMM normCross1 = DOT3(term.cross1, term.cross1);
RealOpenMM normBC = term.delta2[ReferenceForce::RIndex];
forceFactors[0] = (-dEdTheta*normBC)/normCross1;
RealOpenMM 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];
RealOpenMM 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
if (totalEnergy)
*totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
}
void ReferenceCustomCentroidBondIxn::computeDelta(int group1, int group2, RealOpenMM* delta, vector<RealVec>& groupCenters) const {
ReferenceForce::getDeltaR(groupCenters[group1], groupCenters[group2], delta);
}
RealOpenMM ReferenceCustomCentroidBondIxn::computeAngle(RealOpenMM* vec1, RealOpenMM* vec2) {
RealOpenMM dot = DOT3(vec1, vec2);
RealOpenMM cosine = dot/SQRT((vec1[ReferenceForce::R2Index]*vec2[ReferenceForce::R2Index]));
RealOpenMM angle;
if (cosine >= 1)
angle = 0;
else if (cosine <= -1)
angle = PI_M;
else
angle = ACOS(cosine);
return angle;
}
platforms/reference/tests/TestReferenceCustomCentroidBondForce.cpp 0 → 100644
View file @ 56d95193
/* -------------------------------------------------------------------------- *
* 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) 2015 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. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of CustomCompoundBondForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomCentroidBondForce.h"
#include "openmm/CustomCompoundBondForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testHarmonicBond() {
System system;
system.addParticle(1.0);
system.addParticle(2.0);
system.addParticle(3.0);
system.addParticle(4.0);
system.addParticle(5.0);
CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "k*distance(g1,g2)^2");
force->addPerBondParameter("k");
vector<int> particles1;
particles1.push_back(0);
particles1.push_back(1);
vector<int> particles2;
particles2.push_back(2);
particles2.push_back(3);
particles2.push_back(4);
force->addGroup(particles1);
force->addGroup(particles2);
vector<int> groups;
groups.push_back(0);
groups.push_back(1);
vector<double> parameters;
parameters.push_back(1.0);
force->addBond(groups, parameters);
system.addForce(force);
ASSERT(!system.usesPeriodicBoundaryConditions());
// The center of mass of group 0 is (1.5, 0, 0).
vector<Vec3> positions(5);
positions[0] = Vec3(2.5, 0, 0);
positions[1] = Vec3(1, 0, 0);
// The center of mass of group 1 is (-1, 0, 0).
positions[2] = Vec3(-6, 0, 0);
positions[3] = Vec3(-1, 0, 0);
positions[4] = Vec3(2, 0, 0);
// Check the forces and energy.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(2.5*2.5, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(Vec3(-2*2.5*(1.0/3.0), 0, 0), state.getForces()[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-2*2.5*(2.0/3.0), 0, 0), state.getForces()[1], TOL);
ASSERT_EQUAL_VEC(Vec3(2*2.5*(3.0/12.0), 0, 0), state.getForces()[2], TOL);
ASSERT_EQUAL_VEC(Vec3(2*2.5*(4.0/12.0), 0, 0), state.getForces()[3], TOL);
ASSERT_EQUAL_VEC(Vec3(2*2.5*(5.0/12.0), 0, 0), state.getForces()[4], TOL);
// Update the per-bond parameter and see if the results change.
parameters[0] = 2.0;
force->setBondParameters(0, groups, parameters);
force->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(2*2.5*2.5, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(Vec3(-4*2.5*(1.0/3.0), 0, 0), state.getForces()[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-4*2.5*(2.0/3.0), 0, 0), state.getForces()[1], TOL);
ASSERT_EQUAL_VEC(Vec3(4*2.5*(3.0/12.0), 0, 0), state.getForces()[2], TOL);
ASSERT_EQUAL_VEC(Vec3(4*2.5*(4.0/12.0), 0, 0), state.getForces()[3], TOL);
ASSERT_EQUAL_VEC(Vec3(4*2.5*(5.0/12.0), 0, 0), state.getForces()[4], TOL);
}
void testComplexFunction() {
int numParticles = 4;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(2.0);
// When every group contains only one particle, a CustomCentroidBondForce is identical to a
// CustomCompoundBondForce. Use that to test a complicated energy function with lots of terms.
CustomCompoundBondForce* compound = new CustomCompoundBondForce(4, "x1+y2+z4+distance(p1,p2)*angle(p3,p2,p4)+0.5*dihedral(p2,p1,p4,p3)");
CustomCentroidBondForce* centroid = new CustomCentroidBondForce(4, "x1+y2+z4+distance(g1,g2)*angle(g3,g2,g4)+0.5*dihedral(g2,g1,g4,g3)");
// Add a single bond to the CustomCompoundBondForce.
vector<int> particles(4);
vector<double> parameters;
particles[0] = 0;
particles[1] = 1;
particles[2] = 2;
particles[3] = 3;
compound->addBond(particles, parameters);
// Add an identical bond to the CustomCentroidBondForce. As a stronger test, make sure that
// group number is different from particle number.
vector<int> groupMembers(1);
groupMembers[0] = 3;
centroid->addGroup(groupMembers);
groupMembers[0] = 0;
centroid->addGroup(groupMembers);
groupMembers[0] = 1;
centroid->addGroup(groupMembers);
groupMembers[0] = 2;
centroid->addGroup(groupMembers);
vector<int> groups(4);
groups[0] = 1;
groups[1] = 2;
groups[2] = 3;
groups[3] = 0;
centroid->addBond(groups, parameters);
// Add both forces as different force groups, and create a context.
centroid->setForceGroup(1);
system.addForce(compound);
system.addForce(centroid);
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
// Evaluate the force and energy for various positions and see if they match.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < 10; i++) {
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
context.setPositions(positions);
State state1 = context.getState(State::Forces | State::Energy, false, 1<<0);
State state2 = context.getState(State::Forces | State::Energy, false, 1<<1);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), TOL);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], TOL);
}
}
void testCustomWeights() {
System system;
system.addParticle(1.0);
system.addParticle(2.0);
system.addParticle(3.0);
system.addParticle(4.0);
CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "distance(g1,g2)^2");
vector<int> particles(2);
vector<double> weights(2);
particles[0] = 0;
particles[1] = 1;
weights[0] = 0.5;
weights[1] = 1.5;
force->addGroup(particles, weights);
particles[0] = 2;
particles[1] = 3;
weights[0] = 2.0;
weights[1] = 1.0;
force->addGroup(particles, weights);
vector<int> groups;
groups.push_back(0);
groups.push_back(1);
vector<double> parameters;
force->addBond(groups, parameters);
system.addForce(force);
// The center of mass of group 0 is (0, 1, 0).
vector<Vec3> positions(4);
positions[0] = Vec3(0, 4, 0);
positions[1] = Vec3(0, 0, 0);
// The center of mass of group 1 is (0, 10, 0).
positions[2] = Vec3(0, 9, 0);
positions[3] = Vec3(0, 12, 0);
// Check the forces and energy.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(9*9, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(Vec3(0, 2*9*(0.5/2.0), 0), state.getForces()[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 2*9*(1.5/2.0), 0), state.getForces()[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, -2*9*(2.0/3.0), 0), state.getForces()[2], TOL);
ASSERT_EQUAL_VEC(Vec3(0, -2*9*(1.0/3.0), 0), state.getForces()[3], TOL);
}
int main() {
try {
testHarmonicBond();
testComplexFunction();
testCustomWeights();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
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