Skip to content

GitLab

"wrappers/vscode:/vscode.git/clone" did not exist on "d85ccb06e0c8c242560913efdd882f9474ecc774"
Commit 1dd61bab authored by peastman's avatar peastman
Browse files

Reference implementation of parameter offsets for NonbondedForce

parent f9106ddb
Show whitespace changes
Inline Side-by-side
Showing with 357 additions and 33 deletions
openmmapi/include/openmm/NonbondedForce.h
View file @ 1dd61bab
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Portions copyright (c) 2008-2018 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -76,6 +76,32 @@ namespace OpenMM {
* the effect of all Lennard-Jones interactions beyond the cutoff in a periodic system. When running a simulation
* at constant pressure, this can improve the quality of the result. Call setUseDispersionCorrection() to set whether
* this should be used.
*
* In some applications, it is useful to be able to inexpensively change the parameters of small groups of particles.
* Usually this is done to interpolate between two sets of parameters. For example, a titratable group might have
* two states it can exist in, each described by a different set of parameters for the atoms that make up the
* group. You might then want to smoothly interpolate between the two states. This is done by calling
* addParticleParameterOffset() to create a "parameter offset". Each offset defines the following:
*
* <ul>
* <li>A Context parameter used to interpolate between the states.</li>
* <li>A single particle whose parameters are influenced by the Context parameter.</li>
* <li>Three scale factors (chargeScale, sigmaScale, and epsilonScale) that specify how the Context parameter
* affects the particle.</li>
* </ul>
*
* The "effective" parameters for a particle (those used to compute forces) are given by
*
* <tt><pre>
* charge = baseCharge + param*chargeScale
* sigma = baseSigma + param*sigmaScale
* epsilon = baseEpsilon + param*epsilonScale
* </pre></tt>
*
* where the "base" values are the ones specified by addParticle() and "oaram" is the current value
* of the Context parameter. A single Context parameter can apply offsets to multiple particles,
* and multiple parameters can be used to apply offsets to the same particle. Parameters can also be used
* to modify exceptions in exactly the same way by calling addExceptionParameterOffset().
*/
class OPENMM_EXPORT NonbondedForce : public Force {
......@@ -132,6 +158,18 @@ public:
int getNumExceptions() const {
return exceptions.size();
}
/**
* Get the number of particles parameter offsets that have been added.
*/
int getNumParticleParameterOffsets() const {
return particleOffsets.size();
}
/**
* Get the number of exception parameter offsets that have been added.
*/
int getNumExceptionParameterOffsets() const {
return exceptionOffsets.size();
}
/**
* Get the method used for handling long range nonbonded interactions.
*/
......@@ -355,6 +393,78 @@ public:
* multiplied by this factor
*/
void createExceptionsFromBonds(const std::vector<std::pair<int, int> >& bonds, double coulomb14Scale, double lj14Scale);
/**
* Add an offset to the per-particle parameters of a particular particle, based on a global parameter.
*
* @param parameter the name of the global parameter. Its value is initially set to 0, and can be modified
* at any time by calling Context::setParameter().
* @param particleIndex the index of the particle whose parameters are affected
* @param chargeScale this value multiplied by the parameter value is added to the particle's charge
* @param sigmaScale this value multiplied by the parameter value is added to the particle's sigma
* @param epsilonScale this value multiplied by the parameter value is added to the particle's epsilon
* @return the index of the offset that was added
*/
int addParticleParameterOffset(const std::string& parameter, int particleIndex, double chargeScale, double sigmaScale, double epsilonScale);
/**
* Get the offset added to the per-particle parameters of a particular particle, based on a global parameter.
*
* @param index the index of the offset to query, as returned by addParticleParameterOffset()
* @param parameter the name of the global parameter. Its value is initially set to 0, and can be modified
* at any time by calling Context::setParameter().
* @param particleIndex the index of the particle whose parameters are affected
* @param chargeScale this value multiplied by the parameter value is added to the particle's charge
* @param sigmaScale this value multiplied by the parameter value is added to the particle's sigma
* @param epsilonScale this value multiplied by the parameter value is added to the particle's epsilon
*/
void getParticleParameterOffset(int index, std::string& parameter, int& particleIndex, double& chargeScale, double& sigmaScale, double& epsilonScale) const;
/**
* Set the offset added to the per-particle parameters of a particular particle, based on a global parameter.
*
* @param index the index of the offset to modify, as returned by addParticleParameterOffset()
* @param parameter the name of the global parameter. Its value is initially set to 0, and can be modified
* at any time by calling Context::setParameter().
* @param particleIndex the index of the particle whose parameters are affected
* @param chargeScale this value multiplied by the parameter value is added to the particle's charge
* @param sigmaScale this value multiplied by the parameter value is added to the particle's sigma
* @param epsilonScale this value multiplied by the parameter value is added to the particle's epsilon
*/
void setParticleParameterOffset(int index, const std::string& parameter, int particleIndex, double chargeScale, double sigmaScale, double epsilonScale);
/**
* Add an offset to the parameters of a particular exception, based on a global parameter.
*
* @param parameter the name of the global parameter. Its value is initially set to 0, and can be modified
* at any time by calling Context::setParameter().
* @param exceptionIndex the index of the exception whose parameters are affected
* @param chargeProdScale this value multiplied by the parameter value is added to the exception's charge product
* @param sigmaScale this value multiplied by the parameter value is added to the exception's sigma
* @param epsilonScale this value multiplied by the parameter value is added to the exception's epsilon
* @return the index of the offset that was added
*/
int addExceptionParameterOffset(const std::string& parameter, int exceptionIndex, double chargeProdScale, double sigmaScale, double epsilonScale);
/**
* Get the offset added to the parameters of a particular exception, based on a global parameter.
*
* @param index the index of the offset to query, as returned by addExceptionParameterOffset()
* @param parameter the name of the global parameter. Its value is initially set to 0, and can be modified
* at any time by calling Context::setParameter().
* @param exceptionIndex the index of the exception whose parameters are affected
* @param chargeProdScale this value multiplied by the parameter value is added to the exception's charge product
* @param sigmaScale this value multiplied by the parameter value is added to the exception's sigma
* @param epsilonScale this value multiplied by the parameter value is added to the exception's epsilon
*/
void getExceptionParameterOffset(int index, std::string& parameter, int& exceptionIndex, double& chargeProdScale, double& sigmaScale, double& epsilonScale) const;
/**
* Set the offset added to the parameters of a particular exception, based on a global parameter.
*
* @param index the index of the offset to modify, as returned by addExceptionParameterOffset()
* @param parameter the name of the global parameter. Its value is initially set to 0, and can be modified
* at any time by calling Context::setParameter().
* @param exceptionIndex the index of the exception whose parameters are affected
* @param chargeProdScale this value multiplied by the parameter value is added to the exception's charge product
* @param sigmaScale this value multiplied by the parameter value is added to the exception's sigma
* @param epsilonScale this value multiplied by the parameter value is added to the exception's epsilon
*/
void setExceptionParameterOffset(int index, const std::string& parameter, int exceptionIndex, double chargeProdScale, double sigmaScale, double epsilonScale);
/**
* Get whether to add a contribution to the energy that approximately represents the effect of Lennard-Jones
* interactions beyond the cutoff distance. The energy depends on the volume of the periodic box, and is only
......@@ -420,6 +530,8 @@ protected:
private:
class ParticleInfo;
class ExceptionInfo;
class ParticleOffsetInfo;
class ExceptionOffsetInfo;
NonbondedMethod nonbondedMethod;
double cutoffDistance, switchingDistance, rfDielectric, ewaldErrorTol, alpha, dalpha;
bool useSwitchingFunction, useDispersionCorrection;
......@@ -427,6 +539,8 @@ private:
void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
std::vector<ParticleInfo> particles;
std::vector<ExceptionInfo> exceptions;
std::vector<ParticleOffsetInfo> particleOffsets;
std::vector<ExceptionOffsetInfo> exceptionOffsets;
std::map<std::pair<int, int>, int> exceptionMap;
};
......@@ -462,6 +576,42 @@ public:
}
};
/**
* This is an internal class used to record information about a particle parameter offset.
* @private
*/
class NonbondedForce::ParticleOffsetInfo {
public:
std::string parameter;
int particle;
double chargeScale, sigmaScale, epsilonScale;
ParticleOffsetInfo() {
particle = -1;
chargeScale = sigmaScale = epsilonScale = 0.0;
}
ParticleOffsetInfo(const std::string& parameter, int particle, double chargeScale, double sigmaScale, double epsilonScale) :
parameter(parameter), particle(particle), chargeScale(chargeScale), sigmaScale(sigmaScale), epsilonScale(epsilonScale) {
}
};
/**
* This is an internal class used to record information about an exception parameter offset.
* @private
*/
class NonbondedForce::ExceptionOffsetInfo {
public:
std::string parameter;
int exception;
double chargeProdScale, sigmaScale, epsilonScale;
ExceptionOffsetInfo() {
exception = -1;
chargeProdScale = sigmaScale = epsilonScale = 0.0;
}
ExceptionOffsetInfo(const std::string& parameter, int exception, double chargeProdScale, double sigmaScale, double epsilonScale) :
parameter(parameter), exception(exception), chargeProdScale(chargeProdScale), sigmaScale(sigmaScale), epsilonScale(epsilonScale) {
}
};
} // namespace OpenMM
#endif /*OPENMM_NONBONDEDFORCE_H_*/
openmmapi/include/openmm/internal/NonbondedForceImpl.h
View file @ 1dd61bab
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2010 Stanford University and the Authors. *
* Portions copyright (c) 2008-2018 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -59,9 +59,7 @@ public:
// This force field doesn't update the state directly.
}
double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
std::map<std::string, double> getDefaultParameters() {
return std::map<std::string, double>(); // This force field doesn't define any parameters.
}
std::map<std::string, double> getDefaultParameters();
std::vector<std::string> getKernelNames();
void updateParametersInContext(ContextImpl& context);
void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
......
openmmapi/src/NonbondedForce.cpp
View file @ 1dd61bab
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Portions copyright (c) 2008-2018 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -256,6 +256,52 @@ void NonbondedForce::addExclusionsToSet(const vector<set<int> >& bonded12, set<i
}
}
int NonbondedForce::addParticleParameterOffset(const std::string& parameter, int particleIndex, double chargeScale, double sigmaScale, double epsilonScale) {
particleOffsets.push_back(ParticleOffsetInfo(parameter, particleIndex, chargeScale, sigmaScale, epsilonScale));
return particleOffsets.size()-1;
}
void NonbondedForce::getParticleParameterOffset(int index, std::string& parameter, int& particleIndex, double& chargeScale, double& sigmaScale, double& epsilonScale) const {
ASSERT_VALID_INDEX(index, particleOffsets);
parameter = particleOffsets[index].parameter;
particleIndex = particleOffsets[index].particle;
chargeScale = particleOffsets[index].chargeScale;
sigmaScale = particleOffsets[index].sigmaScale;
epsilonScale = particleOffsets[index].epsilonScale;
}
void NonbondedForce::setParticleParameterOffset(int index, const std::string& parameter, int particleIndex, double chargeScale, double sigmaScale, double epsilonScale) {
ASSERT_VALID_INDEX(index, particleOffsets);
particleOffsets[index].parameter = parameter;
particleOffsets[index].particle = particleIndex;
particleOffsets[index].chargeScale = chargeScale;
particleOffsets[index].sigmaScale = sigmaScale;
particleOffsets[index].epsilonScale = epsilonScale;
}
int NonbondedForce::addExceptionParameterOffset(const std::string& parameter, int exceptionIndex, double chargeProdScale, double sigmaScale, double epsilonScale) {
exceptionOffsets.push_back(ExceptionOffsetInfo(parameter, exceptionIndex, chargeProdScale, sigmaScale, epsilonScale));
return exceptionOffsets.size()-1;
}
void NonbondedForce::getExceptionParameterOffset(int index, std::string& parameter, int& exceptionIndex, double& chargeProdScale, double& sigmaScale, double& epsilonScale) const {
ASSERT_VALID_INDEX(index, exceptionOffsets);
parameter = exceptionOffsets[index].parameter;
exceptionIndex = exceptionOffsets[index].exception;
chargeProdScale = exceptionOffsets[index].chargeProdScale;
sigmaScale = exceptionOffsets[index].sigmaScale;
epsilonScale = exceptionOffsets[index].epsilonScale;
}
void NonbondedForce::setExceptionParameterOffset(int index, const std::string& parameter, int exceptionIndex, double chargeProdScale, double sigmaScale, double epsilonScale) {
ASSERT_VALID_INDEX(index, exceptionOffsets);
exceptionOffsets[index].parameter = parameter;
exceptionOffsets[index].exception = exceptionIndex;
exceptionOffsets[index].chargeProdScale = chargeProdScale;
exceptionOffsets[index].sigmaScale = sigmaScale;
exceptionOffsets[index].epsilonScale = epsilonScale;
}
int NonbondedForce::getReciprocalSpaceForceGroup() const {
return recipForceGroup;
}
......
openmmapi/src/NonbondedForceImpl.cpp
View file @ 1dd61bab
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2015 Stanford University and the Authors. *
* Portions copyright (c) 2008-2018 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -110,6 +110,27 @@ double NonbondedForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includ
return kernel.getAs<CalcNonbondedForceKernel>().execute(context, includeForces, includeEnergy, includeDirect, includeReciprocal);
}
map<string, double> NonbondedForceImpl::getDefaultParameters() {
map<string, double> params;
for (int i = 0; i < owner.getNumParticleParameterOffsets(); i++) {
string parameter;
int particle;
double charge, sigma, epsilon;
owner.getParticleParameterOffset(i, parameter, particle, charge, sigma, epsilon);
if (params.find(parameter) == params.end())
params[parameter] = 0.0;
}
for (int i = 0; i < owner.getNumExceptionParameterOffsets(); i++) {
string parameter;
int exception;
double chargeProd, sigma, epsilon;
owner.getExceptionParameterOffset(i, parameter, exception, chargeProd, sigma, epsilon);
if (params.find(parameter) == params.end())
params[parameter] = 0.0;
}
return params;
}
std::vector<std::string> NonbondedForceImpl::getKernelNames() {
std::vector<std::string> names;
names.push_back(CalcNonbondedForceKernel::Name());
......
platforms/reference/include/ReferenceKernels.h
View file @ 1dd61bab
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Portions copyright (c) 2008-2018 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -38,6 +38,8 @@
#include "ReferenceNeighborList.h"
#include "lepton/CompiledExpression.h"
#include "lepton/CustomFunction.h"
#include <array>
#include <utility>
namespace OpenMM {
......@@ -615,9 +617,12 @@ public:
*/
void getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
private:
void computeParameters(ContextImpl& context);
int numParticles, num14;
int **bonded14IndexArray;
double **particleParamArray, **bonded14ParamArray;
std::vector<std::array<double, 3> > baseParticleParams, baseExceptionParams;
std::map<std::pair<std::string, int>, std::array<double, 3> > particleParamOffsets, exceptionParamOffsets;
double nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, ewaldDispersionAlpha, dispersionCoefficient;
int kmax[3], gridSize[3], dispersionGridSize[3];
bool useSwitchingFunction;
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 1dd61bab
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Portions copyright (c) 2008-2018 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -935,23 +935,30 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
bonded14IndexArray = allocateIntArray(num14, 2);
bonded14ParamArray = allocateRealArray(num14, 3);
particleParamArray = allocateRealArray(numParticles, 3);
for (int i = 0; i < numParticles; ++i) {
double charge, radius, depth;
force.getParticleParameters(i, charge, radius, depth);
particleParamArray[i][0] = 0.5*radius;
particleParamArray[i][1] = 2.0*sqrt(depth);
particleParamArray[i][2] = charge;
}
baseParticleParams.resize(numParticles);
baseExceptionParams.resize(num14);
for (int i = 0; i < numParticles; ++i)
force.getParticleParameters(i, baseParticleParams[i][0], baseParticleParams[i][1], baseParticleParams[i][2]);
this->exclusions = exclusions;
for (int i = 0; i < num14; ++i) {
int particle1, particle2;
double charge, radius, depth;
force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
force.getExceptionParameters(nb14s[i], particle1, particle2, baseExceptionParams[i][0], baseExceptionParams[i][1], baseExceptionParams[i][2]);
bonded14IndexArray[i][0] = particle1;
bonded14IndexArray[i][1] = particle2;
bonded14ParamArray[i][0] = radius;
bonded14ParamArray[i][1] = 4.0*depth;
bonded14ParamArray[i][2] = charge;
}
for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
string param;
int particle;
double charge, sigma, epsilon;
force.getParticleParameterOffset(i, param, particle, charge, sigma, epsilon);
particleParamOffsets[make_pair(param, particle)] = {charge, sigma, epsilon};
}
for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
string param;
int exception;
double charge, sigma, epsilon;
force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
exceptionParamOffsets[make_pair(param, exception)] = {charge, sigma, epsilon};
}
nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
nonbondedCutoff = force.getCutoffDistance();
......@@ -990,6 +997,7 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
}
double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
computeParameters(context);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
double energy = 0;
......@@ -1048,22 +1056,13 @@ void ReferenceCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& con
// Record the values.
for (int i = 0; i < numParticles; ++i) {
double charge, radius, depth;
force.getParticleParameters(i, charge, radius, depth);
particleParamArray[i][0] = 0.5*radius;
particleParamArray[i][1] = 2.0*sqrt(depth);
particleParamArray[i][2] = charge;
}
for (int i = 0; i < numParticles; ++i)
force.getParticleParameters(i, baseParticleParams[i][0], baseParticleParams[i][1], baseParticleParams[i][2]);
for (int i = 0; i < num14; ++i) {
int particle1, particle2;
double charge, radius, depth;
force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
force.getExceptionParameters(nb14s[i], particle1, particle2, baseExceptionParams[i][0], baseExceptionParams[i][1], baseExceptionParams[i][2]);
bonded14IndexArray[i][0] = particle1;
bonded14IndexArray[i][1] = particle2;
bonded14ParamArray[i][0] = radius;
bonded14ParamArray[i][1] = 4.0*depth;
bonded14ParamArray[i][2] = charge;
}
// Recompute the coefficient for the dispersion correction.
......@@ -1091,6 +1090,52 @@ void ReferenceCalcNonbondedForceKernel::getLJPMEParameters(double& alpha, int& n
nz = dispersionGridSize[2];
}
void ReferenceCalcNonbondedForceKernel::computeParameters(ContextImpl& context) {
// Compute particle parameters.
vector<double> charges(numParticles), sigmas(numParticles), epsilons(numParticles);
for (int i = 0; i < numParticles; i++) {
charges[i] = baseParticleParams[i][0];
sigmas[i] = baseParticleParams[i][1];
epsilons[i] = baseParticleParams[i][2];
}
for (auto& offset : particleParamOffsets) {
double value = context.getParameter(offset.first.first);
int index = offset.first.second;
charges[index] += value*offset.second[0];
sigmas[index] += value*offset.second[1];
epsilons[index] += value*offset.second[2];
}
for (int i = 0; i < numParticles; i++) {
particleParamArray[i][0] = 0.5*sigmas[i];
particleParamArray[i][1] = 2.0*sqrt(epsilons[i]);
particleParamArray[i][2] = charges[i];
}
// Compute exception parameters.
charges.resize(num14);
sigmas.resize(num14);
epsilons.resize(num14);
for (int i = 0; i < num14; i++) {
charges[i] = baseExceptionParams[i][0];
sigmas[i] = baseExceptionParams[i][1];
epsilons[i] = baseExceptionParams[i][2];
}
for (auto& offset : exceptionParamOffsets) {
double value = context.getParameter(offset.first.first);
int index = offset.first.second;
charges[index] += value*offset.second[0];
sigmas[index] += value*offset.second[1];
epsilons[index] += value*offset.second[2];
}
for (int i = 0; i < num14; i++) {
bonded14ParamArray[i][0] = sigmas[i];
bonded14ParamArray[i][1] = 4.0*epsilons[i];
bonded14ParamArray[i][2] = charges[i];
}
}
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
disposeRealArray(particleParamArray, numParticles);
if (neighborList != NULL)
......
tests/TestNonbondedForce.h
View file @ 1dd61bab
......@@ -734,6 +734,64 @@ void testTwoForces() {
ASSERT_EQUAL_TOL(state1.getPotentialEnergy()+state2.getPotentialEnergy(), state.getPotentialEnergy(), TOL);
}
void testParameterOffsets() {
System system;
for (int i = 0; i < 4; i++)
system.addParticle(1.0);
NonbondedForce* force = new NonbondedForce();
force->addParticle(0.0, 1.0, 0.5);
force->addParticle(1.0, 0.5, 0.6);
force->addParticle(-1.0, 2.0, 0.7);
force->addParticle(0.5, 2.0, 0.8);
force->addException(0, 1, 1.0, 1.5, 1.0);
force->addException(2, 3, 0.5, 1.0, 1.5);
force->addParticleParameterOffset("p1", 0, 3.0, 0.5, 0.5);
force->addParticleParameterOffset("p2", 1, 1.0, 1.0, 2.0);
force->addExceptionParameterOffset("p1", 1, 0.5, 0.5, 1.5);
system.addForce(force);
vector<Vec3> positions(4);
for (int i = 0; i < 4; i++)
positions[i] = Vec3(i, 0, 0);
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
ASSERT_EQUAL(2, context.getParameters().size());
ASSERT_EQUAL(0.0, context.getParameter("p1"));
ASSERT_EQUAL(0.0, context.getParameter("p2"));
context.setParameter("p1", 0.5);
context.setParameter("p2", 1.5);
// Compute the expected parameters for the six interactions.
vector<double> particleCharge = {0.0+3.0*0.5, 1.0+1.0*1.5, -1.0, 0.5};
vector<double> particleSigma = {1.0+0.5*0.5, 0.5+1.0*1.5, 2.0, 2.0};
vector<double> particleEpsilon = {0.5+0.5*0.5, 0.6+2.0*1.5, 0.7, 0.8};
double pairChargeProd[4][4], pairSigma[4][4], pairEpsilon[4][4];
for (int i = 0; i < 4; i++)
for (int j = i+1; j < 4; j++) {
pairChargeProd[i][j] = particleCharge[i]*particleCharge[j];
pairSigma[i][j] = 0.5*(particleSigma[i]+particleSigma[j]);
pairEpsilon[i][j] = sqrt(particleEpsilon[i]*particleEpsilon[j]);
}
pairChargeProd[0][1] = 1.0;
pairSigma[0][1] = 1.5;
pairEpsilon[0][1] = 1.0;
pairChargeProd[2][3] = 0.5+0.5*0.5;
pairSigma[2][3] = 1.0+0.5*0.5;
pairEpsilon[2][3] = 1.5+1.5*0.5;
// Compute the expected energy.
double energy = 0.0;
for (int i = 0; i < 4; i++)
for (int j = i+1; j < 4; j++) {
double dist = j-i;
double x = pairSigma[i][j]/dist;
energy += ONE_4PI_EPS0*pairChargeProd[i][j]/dist + 4.0*pairEpsilon[i][j]*(pow(x, 12.0)-pow(x, 6.0));
}
ASSERT_EQUAL_TOL(energy, context.getState(State::Energy).getPotentialEnergy(), 1e-5);
}
void runPlatformTests();
int main(int argc, char* argv[]) {
......@@ -752,6 +810,7 @@ int main(int argc, char* argv[]) {
testSwitchingFunction(NonbondedForce::CutoffNonPeriodic);
testSwitchingFunction(NonbondedForce::PME);
testTwoForces();
testParameterOffsets();
runPlatformTests();
}
catch(const exception& e) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment