/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2024 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see . *
* -------------------------------------------------------------------------- */
#include "openmm/common/CommonKernels.h"
#include "openmm/common/ContextSelector.h"
#include "openmm/common/ExpressionUtilities.h"
#include "openmm/Context.h"
#include "openmm/internal/AndersenThermostatImpl.h"
#include "openmm/internal/CMAPTorsionForceImpl.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/CustomCentroidBondForceImpl.h"
#include "openmm/internal/CustomCompoundBondForceImpl.h"
#include "openmm/internal/CustomHbondForceImpl.h"
#include "openmm/internal/CustomManyParticleForceImpl.h"
#include "openmm/internal/timer.h"
#include "CommonKernelSources.h"
#include "lepton/CustomFunction.h"
#include "lepton/ExpressionTreeNode.h"
#include "lepton/Operation.h"
#include "lepton/Parser.h"
#include "lepton/ParsedExpression.h"
#include "ReferenceTabulatedFunction.h"
#include "SimTKOpenMMRealType.h"
#include "SimTKOpenMMUtilities.h"
#include "jama_eig.h"
#include
#include
#include
#include
using namespace OpenMM;
using namespace std;
using namespace Lepton;
static void setPeriodicBoxArgs(ComputeContext& cc, ComputeKernel kernel, int index) {
Vec3 a, b, c;
cc.getPeriodicBoxVectors(a, b, c);
if (cc.getUseDoublePrecision()) {
kernel->setArg(index++, mm_double4(a[0], b[1], c[2], 0.0));
kernel->setArg(index++, mm_double4(1.0/a[0], 1.0/b[1], 1.0/c[2], 0.0));
kernel->setArg(index++, mm_double4(a[0], a[1], a[2], 0.0));
kernel->setArg(index++, mm_double4(b[0], b[1], b[2], 0.0));
kernel->setArg(index, mm_double4(c[0], c[1], c[2], 0.0));
}
else {
kernel->setArg(index++, mm_float4((float) a[0], (float) b[1], (float) c[2], 0.0f));
kernel->setArg(index++, mm_float4(1.0f/(float) a[0], 1.0f/(float) b[1], 1.0f/(float) c[2], 0.0f));
kernel->setArg(index++, mm_float4((float) a[0], (float) a[1], (float) a[2], 0.0f));
kernel->setArg(index++, mm_float4((float) b[0], (float) b[1], (float) b[2], 0.0f));
kernel->setArg(index, mm_float4((float) c[0], (float) c[1], (float) c[2], 0.0f));
}
}
static bool isZeroExpression(const Lepton::ParsedExpression& expression) {
const Lepton::Operation& op = expression.getRootNode().getOperation();
if (op.getId() != Lepton::Operation::CONSTANT)
return false;
return (dynamic_cast(op).getValue() == 0.0);
}
static bool usesVariable(const Lepton::ExpressionTreeNode& node, const string& variable) {
const Lepton::Operation& op = node.getOperation();
if (op.getId() == Lepton::Operation::VARIABLE && op.getName() == variable)
return true;
for (auto& child : node.getChildren())
if (usesVariable(child, variable))
return true;
return false;
}
static bool usesVariable(const Lepton::ParsedExpression& expression, const string& variable) {
return usesVariable(expression.getRootNode(), variable);
}
static pair makeVariable(const string& name, const string& value) {
return make_pair(ExpressionTreeNode(new Operation::Variable(name)), value);
}
static void flushPeriodically(ComputeContext& cc) {
#ifdef WIN32
// When running on Windows, we periodically flush the queue to keep the UI responsive.
static double lastTime = getCurrentTime();
double currentTime = getCurrentTime();
if (currentTime-lastTime > 0.025) {
cc.flushQueue();
lastTime = currentTime;
}
#endif
}
void CommonUpdateStateDataKernel::initialize(const System& system) {
}
double CommonUpdateStateDataKernel::getTime(const ContextImpl& context) const {
return cc.getTime();
}
void CommonUpdateStateDataKernel::setTime(ContextImpl& context, double time) {
for (auto ctx : cc.getAllContexts())
ctx->setTime(time);
}
long long CommonUpdateStateDataKernel::getStepCount(const ContextImpl& context) const {
return cc.getStepCount();
}
void CommonUpdateStateDataKernel::setStepCount(const ContextImpl& context, long long count) {
for (auto ctx : cc.getAllContexts())
ctx->setStepCount(count);
}
void CommonUpdateStateDataKernel::getPositions(ContextImpl& context, vector& positions) {
ContextSelector selector(cc);
int numParticles = context.getSystem().getNumParticles();
positions.resize(numParticles);
vector posCorrection;
if (cc.getUseDoublePrecision()) {
mm_double4* posq = (mm_double4*) cc.getPinnedBuffer();
cc.getPosq().download(posq);
}
else if (cc.getUseMixedPrecision()) {
mm_float4* posq = (mm_float4*) cc.getPinnedBuffer();
cc.getPosq().download(posq, false);
posCorrection.resize(numParticles);
cc.getPosqCorrection().download(posCorrection);
}
else {
mm_float4* posq = (mm_float4*) cc.getPinnedBuffer();
cc.getPosq().download(posq);
}
// Filling in the output array is done in parallel for speed.
cc.getThreadPool().execute([&] (ThreadPool& threads, int threadIndex) {
// Compute the position of each particle to return to the user. This is done in parallel for speed.
const vector& order = cc.getAtomIndex();
int numParticles = cc.getNumAtoms();
Vec3 boxVectors[3];
cc.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
int numThreads = threads.getNumThreads();
int start = threadIndex*numParticles/numThreads;
int end = (threadIndex+1)*numParticles/numThreads;
if (cc.getUseDoublePrecision()) {
mm_double4* posq = (mm_double4*) cc.getPinnedBuffer();
for (int i = start; i < end; ++i) {
mm_double4 pos = posq[i];
mm_int4 offset = cc.getPosCellOffsets()[i];
positions[order[i]] = Vec3(pos.x, pos.y, pos.z)-boxVectors[0]*offset.x-boxVectors[1]*offset.y-boxVectors[2]*offset.z;
}
}
else if (cc.getUseMixedPrecision()) {
mm_float4* posq = (mm_float4*) cc.getPinnedBuffer();
for (int i = start; i < end; ++i) {
mm_float4 pos1 = posq[i];
mm_float4 pos2 = posCorrection[i];
mm_int4 offset = cc.getPosCellOffsets()[i];
positions[order[i]] = Vec3((double)pos1.x+(double)pos2.x, (double)pos1.y+(double)pos2.y, (double)pos1.z+(double)pos2.z)-boxVectors[0]*offset.x-boxVectors[1]*offset.y-boxVectors[2]*offset.z;
}
}
else {
mm_float4* posq = (mm_float4*) cc.getPinnedBuffer();
for (int i = start; i < end; ++i) {
mm_float4 pos = posq[i];
mm_int4 offset = cc.getPosCellOffsets()[i];
positions[order[i]] = Vec3(pos.x, pos.y, pos.z)-boxVectors[0]*offset.x-boxVectors[1]*offset.y-boxVectors[2]*offset.z;
}
}
});
cc.getThreadPool().waitForThreads();
}
void CommonUpdateStateDataKernel::setPositions(ContextImpl& context, const vector& positions) {
ContextSelector selector(cc);
const vector& order = cc.getAtomIndex();
int numParticles = context.getSystem().getNumParticles();
if (cc.getUseDoublePrecision()) {
mm_double4* posq = (mm_double4*) cc.getPinnedBuffer();
cc.getPosq().download(posq);
for (int i = 0; i < numParticles; ++i) {
mm_double4& pos = posq[i];
const Vec3& p = positions[order[i]];
pos.x = p[0];
pos.y = p[1];
pos.z = p[2];
}
for (int i = numParticles; i < cc.getPaddedNumAtoms(); i++)
posq[i] = mm_double4(0.0, 0.0, 0.0, 0.0);
cc.getPosq().upload(posq);
}
else {
mm_float4* posq = (mm_float4*) cc.getPinnedBuffer();
cc.getPosq().download(posq);
for (int i = 0; i < numParticles; ++i) {
mm_float4& pos = posq[i];
const Vec3& p = positions[order[i]];
pos.x = (float) p[0];
pos.y = (float) p[1];
pos.z = (float) p[2];
}
for (int i = numParticles; i < cc.getPaddedNumAtoms(); i++)
posq[i] = mm_float4(0.0f, 0.0f, 0.0f, 0.0f);
cc.getPosq().upload(posq);
}
if (cc.getUseMixedPrecision()) {
mm_float4* posCorrection = (mm_float4*) cc.getPinnedBuffer();
for (int i = 0; i < numParticles; ++i) {
mm_float4& c = posCorrection[i];
const Vec3& p = positions[order[i]];
c.x = (float) (p[0]-(float)p[0]);
c.y = (float) (p[1]-(float)p[1]);
c.z = (float) (p[2]-(float)p[2]);
c.w = 0;
}
for (int i = numParticles; i < cc.getPaddedNumAtoms(); i++)
posCorrection[i] = mm_float4(0.0f, 0.0f, 0.0f, 0.0f);
cc.getPosqCorrection().upload(posCorrection);
}
for (auto& offset : cc.getPosCellOffsets())
offset = mm_int4(0, 0, 0, 0);
cc.reorderAtoms();
}
void CommonUpdateStateDataKernel::getVelocities(ContextImpl& context, vector& velocities) {
ContextSelector selector(cc);
const vector& order = cc.getAtomIndex();
int numParticles = context.getSystem().getNumParticles();
velocities.resize(numParticles);
if (cc.getUseDoublePrecision() || cc.getUseMixedPrecision()) {
mm_double4* velm = (mm_double4*) cc.getPinnedBuffer();
cc.getVelm().download(velm);
for (int i = 0; i < numParticles; ++i) {
mm_double4 vel = velm[i];
velocities[order[i]] = Vec3(vel.x, vel.y, vel.z);
}
}
else {
mm_float4* velm = (mm_float4*) cc.getPinnedBuffer();
cc.getVelm().download(velm);
for (int i = 0; i < numParticles; ++i) {
mm_float4 vel = velm[i];
velocities[order[i]] = Vec3(vel.x, vel.y, vel.z);
}
}
}
void CommonUpdateStateDataKernel::setVelocities(ContextImpl& context, const vector& velocities) {
ContextSelector selector(cc);
const vector& order = cc.getAtomIndex();
int numParticles = context.getSystem().getNumParticles();
if (cc.getUseDoublePrecision() || cc.getUseMixedPrecision()) {
mm_double4* velm = (mm_double4*) cc.getPinnedBuffer();
cc.getVelm().download(velm);
for (int i = 0; i < numParticles; ++i) {
mm_double4& vel = velm[i];
const Vec3& p = velocities[order[i]];
vel.x = p[0];
vel.y = p[1];
vel.z = p[2];
}
for (int i = numParticles; i < cc.getPaddedNumAtoms(); i++)
velm[i] = mm_double4(0.0, 0.0, 0.0, 0.0);
cc.getVelm().upload(velm);
}
else {
mm_float4* velm = (mm_float4*) cc.getPinnedBuffer();
cc.getVelm().download(velm);
for (int i = 0; i < numParticles; ++i) {
mm_float4& vel = velm[i];
const Vec3& p = velocities[order[i]];
vel.x = p[0];
vel.y = p[1];
vel.z = p[2];
}
for (int i = numParticles; i < cc.getPaddedNumAtoms(); i++)
velm[i] = mm_float4(0.0f, 0.0f, 0.0f, 0.0f);
cc.getVelm().upload(velm);
}
}
void CommonUpdateStateDataKernel::computeShiftedVelocities(ContextImpl& context, double timeShift, vector& velocities) {
cc.getIntegrationUtilities().computeShiftedVelocities(timeShift, velocities);
}
void CommonUpdateStateDataKernel::getForces(ContextImpl& context, vector& forces) {
ContextSelector selector(cc);
long long* force = (long long*) cc.getPinnedBuffer();
cc.getLongForceBuffer().download(force);
const vector& order = cc.getAtomIndex();
int numParticles = context.getSystem().getNumParticles();
int paddedNumParticles = cc.getPaddedNumAtoms();
forces.resize(numParticles);
double scale = 1.0/(double) 0x100000000LL;
for (int i = 0; i < numParticles; ++i)
forces[order[i]] = Vec3(scale*force[i], scale*force[i+paddedNumParticles], scale*force[i+paddedNumParticles*2]);
}
void CommonUpdateStateDataKernel::getEnergyParameterDerivatives(ContextImpl& context, map& derivs) {
ContextSelector selector(cc);
const vector& paramDerivNames = cc.getEnergyParamDerivNames();
int numDerivs = paramDerivNames.size();
if (numDerivs == 0)
return;
derivs = cc.getEnergyParamDerivWorkspace();
ArrayInterface& derivArray = cc.getEnergyParamDerivBuffer();
if (cc.getUseDoublePrecision() || cc.getUseMixedPrecision()) {
vector derivBuffers;
derivArray.download(derivBuffers);
for (int i = numDerivs; i < derivArray.getSize(); i += numDerivs)
for (int j = 0; j < numDerivs; j++)
derivBuffers[j] += derivBuffers[i+j];
for (int i = 0; i < numDerivs; i++)
derivs[paramDerivNames[i]] += derivBuffers[i];
}
else {
vector derivBuffers;
derivArray.download(derivBuffers);
for (int i = numDerivs; i < derivArray.getSize(); i += numDerivs)
for (int j = 0; j < numDerivs; j++)
derivBuffers[j] += derivBuffers[i+j];
for (int i = 0; i < numDerivs; i++)
derivs[paramDerivNames[i]] += derivBuffers[i];
}
}
void CommonUpdateStateDataKernel::getPeriodicBoxVectors(ContextImpl& context, Vec3& a, Vec3& b, Vec3& c) const {
cc.getPeriodicBoxVectors(a, b, c);
}
void CommonUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) {
if (!cc.getBoxIsTriclinic() && (b[0] != 0 || c[0] != 0 || c[1] != 0))
throw OpenMMException("The box shape has changed from rectangular to triclinic. To do this, you must call setDefaultPeriodicBoxVectors() on the System to specify a triclinic default box, then reinitialize the Context.");
// If any particles have been wrapped to the first periodic box, we need to unwrap them
// to avoid changing their positions.
vector positions;
for (auto offset : cc.getPosCellOffsets()) {
if (offset.x != 0 || offset.y != 0 || offset.z != 0) {
getPositions(context, positions);
break;
}
}
// Update the vectors.
for (auto ctx : cc.getAllContexts())
ctx->setPeriodicBoxVectors(a, b, c);
if (positions.size() > 0)
setPositions(context, positions);
}
void CommonUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostream& stream) {
ContextSelector selector(cc);
int version = 3;
stream.write((char*) &version, sizeof(int));
int precision = (cc.getUseDoublePrecision() ? 2 : cc.getUseMixedPrecision() ? 1 : 0);
stream.write((char*) &precision, sizeof(int));
double time = cc.getTime();
stream.write((char*) &time, sizeof(double));
long long stepCount = cc.getStepCount();
stream.write((char*) &stepCount, sizeof(long long));
int stepsSinceReorder = cc.getStepsSinceReorder();
stream.write((char*) &stepsSinceReorder, sizeof(int));
char* buffer = (char*) cc.getPinnedBuffer();
cc.getPosq().download(buffer);
stream.write(buffer, cc.getPosq().getSize()*cc.getPosq().getElementSize());
if (cc.getUseMixedPrecision()) {
cc.getPosqCorrection().download(buffer);
stream.write(buffer, cc.getPosqCorrection().getSize()*cc.getPosqCorrection().getElementSize());
}
cc.getVelm().download(buffer);
stream.write(buffer, cc.getVelm().getSize()*cc.getVelm().getElementSize());
stream.write((char*) &cc.getAtomIndex()[0], sizeof(int)*cc.getAtomIndex().size());
stream.write((char*) &cc.getPosCellOffsets()[0], sizeof(mm_int4)*cc.getPosCellOffsets().size());
Vec3 boxVectors[3];
cc.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
stream.write((char*) boxVectors, 3*sizeof(Vec3));
cc.getIntegrationUtilities().createCheckpoint(stream);
SimTKOpenMMUtilities::createCheckpoint(stream);
}
void CommonUpdateStateDataKernel::loadCheckpoint(ContextImpl& context, istream& stream) {
ContextSelector selector(cc);
int version;
stream.read((char*) &version, sizeof(int));
if (version != 3)
throw OpenMMException("Checkpoint was created with a different version of OpenMM");
int precision;
stream.read((char*) &precision, sizeof(int));
int expectedPrecision = (cc.getUseDoublePrecision() ? 2 : cc.getUseMixedPrecision() ? 1 : 0);
if (precision != expectedPrecision)
throw OpenMMException("Checkpoint was created with a different numeric precision");
double time;
stream.read((char*) &time, sizeof(double));
long long stepCount;
stream.read((char*) &stepCount, sizeof(long long));
int stepsSinceReorder;
stream.read((char*) &stepsSinceReorder, sizeof(int));
vector contexts = cc.getAllContexts();
for (auto ctx : contexts) {
ctx->setTime(time);
ctx->setStepCount(stepCount);
ctx->setStepsSinceReorder(stepsSinceReorder);
}
char* buffer = (char*) cc.getPinnedBuffer();
stream.read(buffer, cc.getPosq().getSize()*cc.getPosq().getElementSize());
cc.getPosq().upload(buffer);
if (cc.getUseMixedPrecision()) {
stream.read(buffer, cc.getPosqCorrection().getSize()*cc.getPosqCorrection().getElementSize());
cc.getPosqCorrection().upload(buffer);
}
stream.read(buffer, cc.getVelm().getSize()*cc.getVelm().getElementSize());
cc.getVelm().upload(buffer);
stream.read((char*) &cc.getAtomIndex()[0], sizeof(int)*cc.getAtomIndex().size());
cc.getAtomIndexArray().upload(cc.getAtomIndex());
stream.read((char*) &cc.getPosCellOffsets()[0], sizeof(mm_int4)*cc.getPosCellOffsets().size());
Vec3 boxVectors[3];
stream.read((char*) &boxVectors, 3*sizeof(Vec3));
for (auto ctx : contexts)
ctx->setPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
cc.getIntegrationUtilities().loadCheckpoint(stream);
SimTKOpenMMUtilities::loadCheckpoint(stream);
for (auto listener : cc.getReorderListeners())
listener->execute();
cc.validateAtomOrder();
}
void CommonApplyConstraintsKernel::initialize(const System& system) {
}
void CommonApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
ContextSelector selector(cc);
if (!hasInitializedKernel) {
hasInitializedKernel = true;
map defines;
ComputeProgram program = cc.compileProgram(CommonKernelSources::constraints, defines);
applyDeltasKernel = program->createKernel("applyPositionDeltas");
applyDeltasKernel->addArg(cc.getNumAtoms());
applyDeltasKernel->addArg(cc.getPosq());
applyDeltasKernel->addArg(cc.getIntegrationUtilities().getPosDelta());
if (cc.getUseMixedPrecision())
applyDeltasKernel->addArg(cc.getPosqCorrection());
}
IntegrationUtilities& integration = cc.getIntegrationUtilities();
cc.clearBuffer(integration.getPosDelta());
integration.applyConstraints(tol);
applyDeltasKernel->execute(cc.getNumAtoms());
integration.computeVirtualSites();
}
void CommonApplyConstraintsKernel::applyToVelocities(ContextImpl& context, double tol) {
cc.getIntegrationUtilities().applyVelocityConstraints(tol);
}
void CommonVirtualSitesKernel::initialize(const System& system) {
}
void CommonVirtualSitesKernel::computePositions(ContextImpl& context) {
cc.getIntegrationUtilities().computeVirtualSites();
}
class CommonCalcHarmonicBondForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const HarmonicBondForce& force) : force(force) {
}
int getNumParticleGroups() {
return force.getNumBonds();
}
void getParticlesInGroup(int index, vector& particles) {
int particle1, particle2;
double length, k;
force.getBondParameters(index, particle1, particle2, length, k);
particles.resize(2);
particles[0] = particle1;
particles[1] = particle2;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2;
double length1, length2, k1, k2;
force.getBondParameters(group1, particle1, particle2, length1, k1);
force.getBondParameters(group2, particle1, particle2, length2, k2);
return (length1 == length2 && k1 == k2);
}
private:
const HarmonicBondForce& force;
};
void CommonCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumBonds()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumBonds()/numContexts;
numBonds = endIndex-startIndex;
if (numBonds == 0)
return;
vector > atoms(numBonds, vector(2));
params.initialize(cc, numBonds, "bondParams");
vector paramVector(numBonds);
for (int i = 0; i < numBonds; i++) {
double length, k;
force.getBondParameters(startIndex+i, atoms[i][0], atoms[i][1], length, k);
paramVector[i] = mm_float2((float) length, (float) k);
}
params.upload(paramVector);
map replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["COMPUTE_FORCE"] = CommonKernelSources::harmonicBondForce;
replacements["PARAMS"] = cc.getBondedUtilities().addArgument(params, "float2");
cc.getBondedUtilities().addInteraction(atoms, cc.replaceStrings(CommonKernelSources::bondForce, replacements), force.getForceGroup());
info = new ForceInfo(force);
cc.addForce(info);
}
double CommonCalcHarmonicBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
return 0.0;
}
void CommonCalcHarmonicBondForceKernel::copyParametersToContext(ContextImpl& context, const HarmonicBondForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumBonds()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumBonds()/numContexts;
if (numBonds != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of bonds has changed");
if (numBonds == 0)
return;
// Record the per-bond parameters.
vector paramVector(numBonds);
for (int i = 0; i < numBonds; i++) {
int atom1, atom2;
double length, k;
force.getBondParameters(startIndex+i, atom1, atom2, length, k);
paramVector[i] = mm_float2((float) length, (float) k);
}
params.upload(paramVector);
// Mark that the current reordering may be invalid.
cc.invalidateMolecules(info);
}
class CommonCalcCustomBondForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const CustomBondForce& force) : force(force) {
}
int getNumParticleGroups() {
return force.getNumBonds();
}
void getParticlesInGroup(int index, vector& particles) {
int particle1, particle2;
thread_local static vector parameters;
force.getBondParameters(index, particle1, particle2, parameters);
particles.resize(2);
particles[0] = particle1;
particles[1] = particle2;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2;
thread_local static vector parameters1, parameters2;
force.getBondParameters(group1, particle1, particle2, parameters1);
force.getBondParameters(group2, particle1, particle2, parameters2);
for (int i = 0; i < (int) parameters1.size(); i++)
if (parameters1[i] != parameters2[i])
return false;
return true;
}
private:
const CustomBondForce& force;
};
CommonCalcCustomBondForceKernel::~CommonCalcCustomBondForceKernel() {
ContextSelector selector(cc);
if (params != NULL)
delete params;
}
void CommonCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumBonds()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumBonds()/numContexts;
numBonds = endIndex-startIndex;
if (numBonds == 0)
return;
vector > atoms(numBonds, vector(2));
params = new ComputeParameterSet(cc, force.getNumPerBondParameters(), numBonds, "customBondParams");
vector > paramVector(numBonds);
for (int i = 0; i < numBonds; i++)
force.getBondParameters(startIndex+i, atoms[i][0], atoms[i][1], paramVector[i]);
params->setParameterValues(paramVector, true);
info = new ForceInfo(force);
cc.addForce(info);
// Record information for the expressions.
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
globalParamNames[i] = force.getGlobalParameterName(i);
globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
}
Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
Lepton::ParsedExpression forceExpression = energyExpression.differentiate("r").optimize();
map expressions;
expressions["energy += "] = energyExpression;
expressions["real dEdR = "] = forceExpression;
// Create the kernels.
map variables;
variables["r"] = "r";
for (int i = 0; i < force.getNumPerBondParameters(); i++) {
const string& name = force.getPerBondParameterName(i);
variables[name] = "bondParams"+params->getParameterSuffix(i);
}
if (force.getNumGlobalParameters() > 0) {
globals.initialize(cc, force.getNumGlobalParameters(), "customBondGlobals");
globals.upload(globalParamValues);
string argName = cc.getBondedUtilities().addArgument(globals, "float");
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = argName+"["+cc.intToString(i)+"]";
variables[name] = value;
}
}
for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
string paramName = force.getEnergyParameterDerivativeName(i);
string derivVariable = cc.getBondedUtilities().addEnergyParameterDerivative(paramName);
Lepton::ParsedExpression derivExpression = energyExpression.differentiate(paramName).optimize();
expressions[derivVariable+" += "] = derivExpression;
}
stringstream compute;
for (int i = 0; i < (int) params->getParameterInfos().size(); i++) {
ComputeParameterInfo& parameter = params->getParameterInfos()[i];
string argName = cc.getBondedUtilities().addArgument(parameter.getArray(), parameter.getType());
compute< functions;
vector > functionNames;
compute << cc.getExpressionUtilities().createExpressions(expressions, variables, functions, functionNames, "temp");
map replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["COMPUTE_FORCE"] = compute.str();
cc.getBondedUtilities().addInteraction(atoms, cc.replaceStrings(CommonKernelSources::bondForce, replacements), force.getForceGroup());
}
double CommonCalcCustomBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
ContextSelector selector(cc);
if (globals.isInitialized()) {
bool changed = false;
for (int i = 0; i < (int) globalParamNames.size(); i++) {
float value = (float) context.getParameter(globalParamNames[i]);
if (value != globalParamValues[i])
changed = true;
globalParamValues[i] = value;
}
if (changed)
globals.upload(globalParamValues);
}
return 0.0;
}
void CommonCalcCustomBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomBondForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumBonds()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumBonds()/numContexts;
if (numBonds != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of bonds has changed");
if (numBonds == 0)
return;
// Record the per-bond parameters.
vector > paramVector(numBonds);
int atom1, atom2;
for (int i = 0; i < numBonds; i++)
force.getBondParameters(startIndex+i, atom1, atom2, paramVector[i]);
params->setParameterValues(paramVector, true);
// Mark that the current reordering may be invalid.
cc.invalidateMolecules(info);
}
class CommonCalcHarmonicAngleForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const HarmonicAngleForce& force) : force(force) {
}
int getNumParticleGroups() {
return force.getNumAngles();
}
void getParticlesInGroup(int index, vector& particles) {
int particle1, particle2, particle3;
double angle, k;
force.getAngleParameters(index, particle1, particle2, particle3, angle, k);
particles.resize(3);
particles[0] = particle1;
particles[1] = particle2;
particles[2] = particle3;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2, particle3;
double angle1, angle2, k1, k2;
force.getAngleParameters(group1, particle1, particle2, particle3, angle1, k1);
force.getAngleParameters(group2, particle1, particle2, particle3, angle2, k2);
return (angle1 == angle2 && k1 == k2);
}
private:
const HarmonicAngleForce& force;
};
void CommonCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumAngles()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumAngles()/numContexts;
numAngles = endIndex-startIndex;
if (numAngles == 0)
return;
vector > atoms(numAngles, vector(3));
params.initialize(cc, numAngles, "angleParams");
vector paramVector(numAngles);
for (int i = 0; i < numAngles; i++) {
double angle, k;
force.getAngleParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], angle, k);
paramVector[i] = mm_float2((float) angle, (float) k);
}
params.upload(paramVector);
map replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["COMPUTE_FORCE"] = CommonKernelSources::harmonicAngleForce;
replacements["PARAMS"] = cc.getBondedUtilities().addArgument(params, "float2");
cc.getBondedUtilities().addInteraction(atoms, cc.replaceStrings(CommonKernelSources::angleForce, replacements), force.getForceGroup());
info = new ForceInfo(force);
cc.addForce(info);
}
double CommonCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
return 0.0;
}
void CommonCalcHarmonicAngleForceKernel::copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumAngles()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumAngles()/numContexts;
if (numAngles != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of angles has changed");
if (numAngles == 0)
return;
// Record the per-angle parameters.
vector paramVector(numAngles);
for (int i = 0; i < numAngles; i++) {
int atom1, atom2, atom3;
double angle, k;
force.getAngleParameters(startIndex+i, atom1, atom2, atom3, angle, k);
paramVector[i] = mm_float2((float) angle, (float) k);
}
params.upload(paramVector);
// Mark that the current reordering may be invalid.
cc.invalidateMolecules();
}
class CommonCalcCustomAngleForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const CustomAngleForce& force) : force(force) {
}
int getNumParticleGroups() {
return force.getNumAngles();
}
void getParticlesInGroup(int index, vector& particles) {
int particle1, particle2, particle3;
thread_local static vector parameters;
force.getAngleParameters(index, particle1, particle2, particle3, parameters);
particles.resize(3);
particles[0] = particle1;
particles[1] = particle2;
particles[2] = particle3;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2, particle3;
thread_local static vector parameters1, parameters2;
force.getAngleParameters(group1, particle1, particle2, particle3, parameters1);
force.getAngleParameters(group2, particle1, particle2, particle3, parameters2);
for (int i = 0; i < (int) parameters1.size(); i++)
if (parameters1[i] != parameters2[i])
return false;
return true;
}
private:
const CustomAngleForce& force;
};
CommonCalcCustomAngleForceKernel::~CommonCalcCustomAngleForceKernel() {
ContextSelector selector(cc);
if (params != NULL)
delete params;
}
void CommonCalcCustomAngleForceKernel::initialize(const System& system, const CustomAngleForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumAngles()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumAngles()/numContexts;
numAngles = endIndex-startIndex;
if (numAngles == 0)
return;
vector > atoms(numAngles, vector(3));
params = new ComputeParameterSet(cc, force.getNumPerAngleParameters(), numAngles, "customAngleParams");
vector > paramVector(numAngles);
for (int i = 0; i < numAngles; i++)
force.getAngleParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], paramVector[i]);
params->setParameterValues(paramVector, true);
info = new ForceInfo(force);
cc.addForce(info);
// Record information for the expressions.
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
globalParamNames[i] = force.getGlobalParameterName(i);
globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
}
Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
Lepton::ParsedExpression forceExpression = energyExpression.differentiate("theta").optimize();
map expressions;
expressions["energy += "] = energyExpression;
expressions["real dEdAngle = "] = forceExpression;
// Create the kernels.
map variables;
variables["theta"] = "theta";
for (int i = 0; i < force.getNumPerAngleParameters(); i++) {
const string& name = force.getPerAngleParameterName(i);
variables[name] = "angleParams"+params->getParameterSuffix(i);
}
if (force.getNumGlobalParameters() > 0) {
globals.initialize(cc, force.getNumGlobalParameters(), "customAngleGlobals");
globals.upload(globalParamValues);
string argName = cc.getBondedUtilities().addArgument(globals, "float");
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = argName+"["+cc.intToString(i)+"]";
variables[name] = value;
}
}
for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
string paramName = force.getEnergyParameterDerivativeName(i);
string derivVariable = cc.getBondedUtilities().addEnergyParameterDerivative(paramName);
Lepton::ParsedExpression derivExpression = energyExpression.differentiate(paramName).optimize();
expressions[derivVariable+" += "] = derivExpression;
}
stringstream compute;
for (int i = 0; i < (int) params->getParameterInfos().size(); i++) {
ComputeParameterInfo& parameter = params->getParameterInfos()[i];
string argName = cc.getBondedUtilities().addArgument(parameter.getArray(), parameter.getType());
compute< functions;
vector > functionNames;
compute << cc.getExpressionUtilities().createExpressions(expressions, variables, functions, functionNames, "temp");
map replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["COMPUTE_FORCE"] = compute.str();
cc.getBondedUtilities().addInteraction(atoms, cc.replaceStrings(CommonKernelSources::angleForce, replacements), force.getForceGroup());
}
double CommonCalcCustomAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
ContextSelector selector(cc);
if (globals.isInitialized()) {
bool changed = false;
for (int i = 0; i < (int) globalParamNames.size(); i++) {
float value = (float) context.getParameter(globalParamNames[i]);
if (value != globalParamValues[i])
changed = true;
globalParamValues[i] = value;
}
if (changed)
globals.upload(globalParamValues);
}
return 0.0;
}
void CommonCalcCustomAngleForceKernel::copyParametersToContext(ContextImpl& context, const CustomAngleForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumAngles()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumAngles()/numContexts;
if (numAngles != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of angles has changed");
if (numAngles == 0)
return;
// Record the per-angle parameters.
vector > paramVector(numAngles);
int atom1, atom2, atom3;
for (int i = 0; i < numAngles; i++)
force.getAngleParameters(startIndex+i, atom1, atom2, atom3, paramVector[i]);
params->setParameterValues(paramVector, true);
// Mark that the current reordering may be invalid.
cc.invalidateMolecules(info);
}
class CommonCalcPeriodicTorsionForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const PeriodicTorsionForce& force) : force(force) {
}
int getNumParticleGroups() {
return force.getNumTorsions();
}
void getParticlesInGroup(int index, vector& particles) {
int particle1, particle2, particle3, particle4, periodicity;
double phase, k;
force.getTorsionParameters(index, particle1, particle2, particle3, particle4, periodicity, phase, k);
particles.resize(4);
particles[0] = particle1;
particles[1] = particle2;
particles[2] = particle3;
particles[3] = particle4;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2, particle3, particle4, periodicity1, periodicity2;
double phase1, phase2, k1, k2;
force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, periodicity1, phase1, k1);
force.getTorsionParameters(group2, particle1, particle2, particle3, particle4, periodicity2, phase2, k2);
return (periodicity1 == periodicity2 && phase1 == phase2 && k1 == k2);
}
private:
const PeriodicTorsionForce& force;
};
void CommonCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumTorsions()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumTorsions()/numContexts;
numTorsions = endIndex-startIndex;
if (numTorsions == 0)
return;
vector > atoms(numTorsions, vector(4));
params.initialize(cc, numTorsions, "periodicTorsionParams");
vector paramVector(numTorsions);
for (int i = 0; i < numTorsions; i++) {
int periodicity;
double phase, k;
force.getTorsionParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], atoms[i][3], periodicity, phase, k);
paramVector[i] = mm_float4((float) k, (float) phase, (float) periodicity, 0.0f);
}
params.upload(paramVector);
map replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["COMPUTE_FORCE"] = CommonKernelSources::periodicTorsionForce;
replacements["PARAMS"] = cc.getBondedUtilities().addArgument(params, "float4");
cc.getBondedUtilities().addInteraction(atoms, cc.replaceStrings(CommonKernelSources::torsionForce, replacements), force.getForceGroup());
info = new ForceInfo(force);
cc.addForce(info);
}
double CommonCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
return 0.0;
}
void CommonCalcPeriodicTorsionForceKernel::copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumTorsions()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumTorsions()/numContexts;
if (numTorsions != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
if (numTorsions == 0)
return;
// Record the per-torsion parameters.
vector paramVector(numTorsions);
for (int i = 0; i < numTorsions; i++) {
int atom1, atom2, atom3, atom4, periodicity;
double phase, k;
force.getTorsionParameters(startIndex+i, atom1, atom2, atom3, atom4, periodicity, phase, k);
paramVector[i] = mm_float4((float) k, (float) phase, (float) periodicity, 0.0f);
}
params.upload(paramVector);
// Mark that the current reordering may be invalid.
cc.invalidateMolecules();
}
class CommonCalcRBTorsionForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const RBTorsionForce& force) : force(force) {
}
int getNumParticleGroups() {
return force.getNumTorsions();
}
void getParticlesInGroup(int index, vector& particles) {
int particle1, particle2, particle3, particle4;
double c0, c1, c2, c3, c4, c5;
force.getTorsionParameters(index, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
particles.resize(4);
particles[0] = particle1;
particles[1] = particle2;
particles[2] = particle3;
particles[3] = particle4;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2, particle3, particle4;
double c0a, c0b, c1a, c1b, c2a, c2b, c3a, c3b, c4a, c4b, c5a, c5b;
force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, c0a, c1a, c2a, c3a, c4a, c5a);
force.getTorsionParameters(group2, particle1, particle2, particle3, particle4, c0b, c1b, c2b, c3b, c4b, c5b);
return (c0a == c0b && c1a == c1b && c2a == c2b && c3a == c3b && c4a == c4b && c5a == c5b);
}
private:
const RBTorsionForce& force;
};
void CommonCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumTorsions()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumTorsions()/numContexts;
numTorsions = endIndex-startIndex;
if (numTorsions == 0)
return;
vector > atoms(numTorsions, vector(4));
params1.initialize(cc, numTorsions, "rbTorsionParams1");
params2.initialize(cc, numTorsions, "rbTorsionParams2");
vector paramVector1(numTorsions);
vector paramVector2(numTorsions);
for (int i = 0; i < numTorsions; i++) {
double c0, c1, c2, c3, c4, c5;
force.getTorsionParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], atoms[i][3], c0, c1, c2, c3, c4, c5);
paramVector1[i] = mm_float4((float) c0, (float) c1, (float) c2, (float) c3);
paramVector2[i] = mm_float2((float) c4, (float) c5);
}
params1.upload(paramVector1);
params2.upload(paramVector2);
map replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["COMPUTE_FORCE"] = CommonKernelSources::rbTorsionForce;
replacements["PARAMS1"] = cc.getBondedUtilities().addArgument(params1, "float4");
replacements["PARAMS2"] = cc.getBondedUtilities().addArgument(params2, "float2");
cc.getBondedUtilities().addInteraction(atoms, cc.replaceStrings(CommonKernelSources::torsionForce, replacements), force.getForceGroup());
info = new ForceInfo(force);
cc.addForce(info);
}
double CommonCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
return 0.0;
}
void CommonCalcRBTorsionForceKernel::copyParametersToContext(ContextImpl& context, const RBTorsionForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumTorsions()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumTorsions()/numContexts;
if (numTorsions != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
if (numTorsions == 0)
return;
// Record the per-torsion parameters.
vector paramVector1(numTorsions);
vector paramVector2(numTorsions);
for (int i = 0; i < numTorsions; i++) {
int atom1, atom2, atom3, atom4;
double c0, c1, c2, c3, c4, c5;
force.getTorsionParameters(startIndex+i, atom1, atom2, atom3, atom4, c0, c1, c2, c3, c4, c5);
paramVector1[i] = mm_float4((float) c0, (float) c1, (float) c2, (float) c3);
paramVector2[i] = mm_float2((float) c4, (float) c5);
}
params1.upload(paramVector1);
params2.upload(paramVector2);
// Mark that the current reordering may be invalid.
cc.invalidateMolecules();
}
class CommonCalcCustomTorsionForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const CustomTorsionForce& force) : force(force) {
}
int getNumParticleGroups() {
return force.getNumTorsions();
}
void getParticlesInGroup(int index, vector& particles) {
int particle1, particle2, particle3, particle4;
thread_local static vector parameters;
force.getTorsionParameters(index, particle1, particle2, particle3, particle4, parameters);
particles.resize(4);
particles[0] = particle1;
particles[1] = particle2;
particles[2] = particle3;
particles[3] = particle4;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2, particle3, particle4;
thread_local static vector parameters1, parameters2;
force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, parameters1);
force.getTorsionParameters(group2, particle1, particle2, particle3, particle4, parameters2);
for (int i = 0; i < (int) parameters1.size(); i++)
if (parameters1[i] != parameters2[i])
return false;
return true;
}
private:
const CustomTorsionForce& force;
};
CommonCalcCustomTorsionForceKernel::~CommonCalcCustomTorsionForceKernel() {
if (params != NULL)
delete params;
}
void CommonCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumTorsions()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumTorsions()/numContexts;
numTorsions = endIndex-startIndex;
if (numTorsions == 0)
return;
vector > atoms(numTorsions, vector(4));
params = new ComputeParameterSet(cc, force.getNumPerTorsionParameters(), numTorsions, "customTorsionParams");
vector > paramVector(numTorsions);
for (int i = 0; i < numTorsions; i++)
force.getTorsionParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], atoms[i][3], paramVector[i]);
params->setParameterValues(paramVector, true);
info = new ForceInfo(force);
cc.addForce(info);
// Record information for the expressions.
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
globalParamNames[i] = force.getGlobalParameterName(i);
globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
}
Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
Lepton::ParsedExpression forceExpression = energyExpression.differentiate("theta").optimize();
map expressions;
expressions["energy += "] = energyExpression;
expressions["real dEdAngle = "] = forceExpression;
// Create the kernels.
map variables;
variables["theta"] = "theta";
for (int i = 0; i < force.getNumPerTorsionParameters(); i++) {
const string& name = force.getPerTorsionParameterName(i);
variables[name] = "torsionParams"+params->getParameterSuffix(i);
}
if (force.getNumGlobalParameters() > 0) {
globals.initialize(cc, force.getNumGlobalParameters(), "customTorsionGlobals");
globals.upload(globalParamValues);
string argName = cc.getBondedUtilities().addArgument(globals, "float");
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = argName+"["+cc.intToString(i)+"]";
variables[name] = value;
}
}
for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
string paramName = force.getEnergyParameterDerivativeName(i);
string derivVariable = cc.getBondedUtilities().addEnergyParameterDerivative(paramName);
Lepton::ParsedExpression derivExpression = energyExpression.differentiate(paramName).optimize();
expressions[derivVariable+" += "] = derivExpression;
}
stringstream compute;
for (int i = 0; i < (int) params->getParameterInfos().size(); i++) {
ComputeParameterInfo& parameter = params->getParameterInfos()[i];
string argName = cc.getBondedUtilities().addArgument(parameter.getArray(), parameter.getType());
compute< functions;
vector > functionNames;
compute << cc.getExpressionUtilities().createExpressions(expressions, variables, functions, functionNames, "temp");
map replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["COMPUTE_FORCE"] = compute.str();
cc.getBondedUtilities().addInteraction(atoms, cc.replaceStrings(CommonKernelSources::torsionForce, replacements), force.getForceGroup());
}
double CommonCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
ContextSelector selector(cc);
if (globals.isInitialized()) {
bool changed = false;
for (int i = 0; i < (int) globalParamNames.size(); i++) {
float value = (float) context.getParameter(globalParamNames[i]);
if (value != globalParamValues[i])
changed = true;
globalParamValues[i] = value;
}
if (changed)
globals.upload(globalParamValues);
}
return 0.0;
}
void CommonCalcCustomTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumTorsions()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumTorsions()/numContexts;
if (numTorsions != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
if (numTorsions == 0)
return;
// Record the per-torsion parameters.
vector > paramVector(numTorsions);
int atom1, atom2, atom3, atom4;
for (int i = 0; i < numTorsions; i++)
force.getTorsionParameters(startIndex+i, atom1, atom2, atom3, atom4, paramVector[i]);
params->setParameterValues(paramVector, true);
// Mark that the current reordering may be invalid.
cc.invalidateMolecules(info);
}
class CommonCalcCMAPTorsionForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const CMAPTorsionForce& force) : force(force) {
}
int getNumParticleGroups() {
return force.getNumTorsions();
}
void getParticlesInGroup(int index, vector& particles) {
int map, a1, a2, a3, a4, b1, b2, b3, b4;
force.getTorsionParameters(index, map, a1, a2, a3, a4, b1, b2, b3, b4);
particles.resize(8);
particles[0] = a1;
particles[1] = a2;
particles[2] = a3;
particles[3] = a4;
particles[4] = b1;
particles[5] = b2;
particles[6] = b3;
particles[7] = b4;
}
bool areGroupsIdentical(int group1, int group2) {
int map1, map2, a1, a2, a3, a4, b1, b2, b3, b4;
force.getTorsionParameters(group1, map1, a1, a2, a3, a4, b1, b2, b3, b4);
force.getTorsionParameters(group2, map2, a1, a2, a3, a4, b1, b2, b3, b4);
return (map1 == map2);
}
private:
const CMAPTorsionForce& force;
};
void CommonCalcCMAPTorsionForceKernel::initialize(const System& system, const CMAPTorsionForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumTorsions()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumTorsions()/numContexts;
numTorsions = endIndex-startIndex;
if (numTorsions == 0)
return;
int numMaps = force.getNumMaps();
vector coeffVec;
mapPositionsVec.resize(numMaps);
vector energy;
vector > c;
int currentPosition = 0;
for (int i = 0; i < numMaps; i++) {
int size;
force.getMapParameters(i, size, energy);
CMAPTorsionForceImpl::calcMapDerivatives(size, energy, c);
mapPositionsVec[i] = mm_int2(currentPosition, size);
currentPosition += 4*size*size;
for (int j = 0; j < size*size; j++) {
coeffVec.push_back(mm_float4((float) c[j][0], (float) c[j][1], (float) c[j][2], (float) c[j][3]));
coeffVec.push_back(mm_float4((float) c[j][4], (float) c[j][5], (float) c[j][6], (float) c[j][7]));
coeffVec.push_back(mm_float4((float) c[j][8], (float) c[j][9], (float) c[j][10], (float) c[j][11]));
coeffVec.push_back(mm_float4((float) c[j][12], (float) c[j][13], (float) c[j][14], (float) c[j][15]));
}
}
vector > atoms(numTorsions, vector(8));
vector torsionMapsVec(numTorsions);
for (int i = 0; i < numTorsions; i++)
force.getTorsionParameters(startIndex+i, torsionMapsVec[i], atoms[i][0], atoms[i][1], atoms[i][2], atoms[i][3], atoms[i][4], atoms[i][5], atoms[i][6], atoms[i][7]);
coefficients.initialize(cc, coeffVec.size(), "cmapTorsionCoefficients");
mapPositions.initialize(cc, numMaps, "cmapTorsionMapPositions");
torsionMaps.initialize(cc, numTorsions, "cmapTorsionMaps");
coefficients.upload(coeffVec);
mapPositions.upload(mapPositionsVec);
torsionMaps.upload(torsionMapsVec);
map replacements;
replacements["APPLY_PERIODIC"] = (force.usesPeriodicBoundaryConditions() ? "1" : "0");
replacements["COEFF"] = cc.getBondedUtilities().addArgument(coefficients, "float4");
replacements["MAP_POS"] = cc.getBondedUtilities().addArgument(mapPositions, "int2");
replacements["MAPS"] = cc.getBondedUtilities().addArgument(torsionMaps, "int");
cc.getBondedUtilities().addInteraction(atoms, cc.replaceStrings(CommonKernelSources::cmapTorsionForce, replacements), force.getForceGroup());
info = new ForceInfo(force);
cc.addForce(info);
}
double CommonCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
return 0.0;
}
void CommonCalcCMAPTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force) {
int numMaps = force.getNumMaps();
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumTorsions()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumTorsions()/numContexts;
numTorsions = endIndex-startIndex;
if (mapPositions.getSize() != numMaps)
throw OpenMMException("updateParametersInContext: The number of maps has changed");
if (torsionMaps.getSize() != numTorsions)
throw OpenMMException("updateParametersInContext: The number of CMAP torsions has changed");
// Update the maps.
ContextSelector selector(cc);
vector coeffVec;
vector energy;
vector > c;
int currentPosition = 0;
for (int i = 0; i < numMaps; i++) {
int size;
force.getMapParameters(i, size, energy);
if (size != mapPositionsVec[i].y)
throw OpenMMException("updateParametersInContext: The size of a map has changed");
CMAPTorsionForceImpl::calcMapDerivatives(size, energy, c);
currentPosition += 4*size*size;
for (int j = 0; j < size*size; j++) {
coeffVec.push_back(mm_float4((float) c[j][0], (float) c[j][1], (float) c[j][2], (float) c[j][3]));
coeffVec.push_back(mm_float4((float) c[j][4], (float) c[j][5], (float) c[j][6], (float) c[j][7]));
coeffVec.push_back(mm_float4((float) c[j][8], (float) c[j][9], (float) c[j][10], (float) c[j][11]));
coeffVec.push_back(mm_float4((float) c[j][12], (float) c[j][13], (float) c[j][14], (float) c[j][15]));
}
}
coefficients.upload(coeffVec);
// Update the indices.
vector torsionMapsVec(numTorsions);
for (int i = 0; i < numTorsions; i++) {
int index[8];
force.getTorsionParameters(i, torsionMapsVec[i], index[0], index[1], index[2], index[3], index[4], index[5], index[6], index[7]);
}
torsionMaps.upload(torsionMapsVec);
}
class CommonCalcCustomExternalForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const CustomExternalForce& force, int numParticles) : force(force), indices(numParticles, -1) {
vector params;
for (int i = 0; i < force.getNumParticles(); i++) {
int particle;
force.getParticleParameters(i, particle, params);
indices[particle] = i;
}
}
bool areParticlesIdentical(int particle1, int particle2) {
particle1 = indices[particle1];
particle2 = indices[particle2];
if (particle1 == -1 && particle2 == -1)
return true;
if (particle1 == -1 || particle2 == -1)
return false;
int temp;
thread_local static vector params1, params2;
force.getParticleParameters(particle1, temp, params1);
force.getParticleParameters(particle2, temp, params2);
for (int i = 0; i < (int) params1.size(); i++)
if (params1[i] != params2[i])
return false;
return true;
}
private:
const CustomExternalForce& force;
vector indices;
};
CommonCalcCustomExternalForceKernel::~CommonCalcCustomExternalForceKernel() {
ContextSelector selector(cc);
if (params != NULL)
delete params;
}
void CommonCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumParticles()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumParticles()/numContexts;
numParticles = endIndex-startIndex;
if (numParticles == 0)
return;
vector > atoms(numParticles, vector(1));
params = new ComputeParameterSet(cc, force.getNumPerParticleParameters(), numParticles, "customExternalParams");
vector > paramVector(numParticles);
for (int i = 0; i < numParticles; i++)
force.getParticleParameters(startIndex+i, atoms[i][0], paramVector[i]);
params->setParameterValues(paramVector, true);
info = new ForceInfo(force, system.getNumParticles());
cc.addForce(info);
// Record information for the expressions.
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
globalParamNames[i] = force.getGlobalParameterName(i);
globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
}
map customFunctions;
customFunctions["periodicdistance"] = cc.getExpressionUtilities().getPeriodicDistancePlaceholder();
Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction(), customFunctions).optimize();
Lepton::ParsedExpression forceExpressionX = energyExpression.differentiate("x").optimize();
Lepton::ParsedExpression forceExpressionY = energyExpression.differentiate("y").optimize();
Lepton::ParsedExpression forceExpressionZ = energyExpression.differentiate("z").optimize();
map expressions;
expressions["energy += "] = energyExpression;
expressions["real dEdX = "] = forceExpressionX;
expressions["real dEdY = "] = forceExpressionY;
expressions["real dEdZ = "] = forceExpressionZ;
// Create the kernels.
map variables;
variables["x"] = "pos1.x";
variables["y"] = "pos1.y";
variables["z"] = "pos1.z";
for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
const string& name = force.getPerParticleParameterName(i);
variables[name] = "particleParams"+params->getParameterSuffix(i);
}
if (force.getNumGlobalParameters() > 0) {
globals.initialize(cc, force.getNumGlobalParameters(), "customExternalGlobals");
globals.upload(globalParamValues);
string argName = cc.getBondedUtilities().addArgument(globals, "float");
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = argName+"["+cc.intToString(i)+"]";
variables[name] = value;
}
}
stringstream compute;
for (int i = 0; i < (int) params->getParameterInfos().size(); i++) {
ComputeParameterInfo& parameter = params->getParameterInfos()[i];
string argName = cc.getBondedUtilities().addArgument(parameter.getArray(), parameter.getType());
compute< functions;
vector > functionNames;
compute << cc.getExpressionUtilities().createExpressions(expressions, variables, functions, functionNames, "temp");
map replacements;
replacements["COMPUTE_FORCE"] = compute.str();
cc.getBondedUtilities().addInteraction(atoms, cc.replaceStrings(CommonKernelSources::customExternalForce, replacements), force.getForceGroup());
}
double CommonCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
ContextSelector selector(cc);
if (globals.isInitialized()) {
bool changed = false;
for (int i = 0; i < (int) globalParamNames.size(); i++) {
float value = (float) context.getParameter(globalParamNames[i]);
if (value != globalParamValues[i])
changed = true;
globalParamValues[i] = value;
}
if (changed)
globals.upload(globalParamValues);
}
return 0.0;
}
void CommonCalcCustomExternalForceKernel::copyParametersToContext(ContextImpl& context, const CustomExternalForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumParticles()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumParticles()/numContexts;
if (numParticles != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of particles has changed");
if (numParticles == 0)
return;
// Record the per-particle parameters.
vector > paramVector(numParticles);
int particle;
for (int i = 0; i < numParticles; i++)
force.getParticleParameters(startIndex+i, particle, paramVector[i]);
params->setParameterValues(paramVector, true);
// Mark that the current reordering may be invalid.
cc.invalidateMolecules(info);
}
class CommonCalcCustomCompoundBondForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const CustomCompoundBondForce& force) : force(force) {
}
int getNumParticleGroups() {
return force.getNumBonds();
}
void getParticlesInGroup(int index, vector& particles) {
thread_local static vector parameters;
force.getBondParameters(index, particles, parameters);
}
bool areGroupsIdentical(int group1, int group2) {
thread_local static vector particles;
thread_local static vector parameters1, parameters2;
force.getBondParameters(group1, particles, parameters1);
force.getBondParameters(group2, particles, parameters2);
for (int i = 0; i < (int) parameters1.size(); i++)
if (parameters1[i] != parameters2[i])
return false;
return true;
}
private:
const CustomCompoundBondForce& force;
};
CommonCalcCustomCompoundBondForceKernel::~CommonCalcCustomCompoundBondForceKernel() {
ContextSelector selector(cc);
if (params != NULL)
delete params;
}
void CommonCalcCustomCompoundBondForceKernel::initialize(const System& system, const CustomCompoundBondForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumBonds()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumBonds()/numContexts;
numBonds = endIndex-startIndex;
if (numBonds == 0)
return;
int particlesPerBond = force.getNumParticlesPerBond();
vector > atoms(numBonds, vector(particlesPerBond));
params = new ComputeParameterSet(cc, force.getNumPerBondParameters(), numBonds, "customCompoundBondParams", false, cc.getUseDoublePrecision());
vector > paramVector(numBonds);
for (int i = 0; i < numBonds; i++)
force.getBondParameters(startIndex+i, atoms[i], paramVector[i]);
params->setParameterValues(paramVector, true);
info = new ForceInfo(force);
cc.addForce(info);
// Record the tabulated functions.
map functions;
vector > functionDefinitions;
vector functionList;
tabulatedFunctionArrays.resize(force.getNumTabulatedFunctions());
for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
functionList.push_back(&force.getTabulatedFunction(i));
string name = force.getTabulatedFunctionName(i);
tabulatedFunctionUpdateCount[name] = force.getTabulatedFunction(i).getUpdateCount();
functions[name] = cc.getExpressionUtilities().getFunctionPlaceholder(force.getTabulatedFunction(i));
int width;
vector f = cc.getExpressionUtilities().computeFunctionCoefficients(force.getTabulatedFunction(i), width);
tabulatedFunctionArrays[i].initialize(cc, f.size(), "TabulatedFunction");
tabulatedFunctionArrays[i].upload(f);
string arrayName = cc.getBondedUtilities().addArgument(tabulatedFunctionArrays[i], width == 1 ? "float" : "float"+cc.intToString(width));
functionDefinitions.push_back(make_pair(name, arrayName));
}
// Record information about parameters.
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
globalParamNames[i] = force.getGlobalParameterName(i);
globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
}
map variables;
for (int i = 0; i < particlesPerBond; i++) {
string index = cc.intToString(i+1);
variables["x"+index] = "pos"+index+".x";
variables["y"+index] = "pos"+index+".y";
variables["z"+index] = "pos"+index+".z";
}
for (int i = 0; i < force.getNumPerBondParameters(); i++) {
const string& name = force.getPerBondParameterName(i);
variables[name] = "bondParams"+params->getParameterSuffix(i);
}
if (force.getNumGlobalParameters() > 0) {
globals.initialize(cc, force.getNumGlobalParameters(), "customCompoundBondGlobals");
globals.upload(globalParamValues);
string argName = cc.getBondedUtilities().addArgument(globals, "float");
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = argName+"["+cc.intToString(i)+"]";
variables[name] = value;
}
}
// Generate the kernel.
Lepton::ParsedExpression energyExpression = CustomCompoundBondForceImpl::prepareExpression(force, functions);
map forceExpressions;
stringstream compute;
for (int i = 0; i < (int) params->getParameterInfos().size(); i++) {
ComputeParameterInfo& parameter = params->getParameterInfos()[i];
string argName = cc.getBondedUtilities().addArgument(parameter.getArray(), parameter.getType());
compute< forceNames;
for (int i = 0; i < particlesPerBond; i++) {
string istr = cc.intToString(i+1);
string forceName = "force"+istr;
forceNames.push_back(forceName);
compute<<"real3 "< replacements;
replacements["M_PI"] = cc.doubleToString(M_PI);
cc.getBondedUtilities().addPrefixCode(cc.replaceStrings(CommonKernelSources::pointFunctions, replacements));
}
double CommonCalcCustomCompoundBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
ContextSelector selector(cc);
if (globals.isInitialized()) {
bool changed = false;
for (int i = 0; i < (int) globalParamNames.size(); i++) {
float value = (float) context.getParameter(globalParamNames[i]);
if (value != globalParamValues[i])
changed = true;
globalParamValues[i] = value;
}
if (changed)
globals.upload(globalParamValues);
}
return 0.0;
}
void CommonCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force) {
ContextSelector selector(cc);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*force.getNumBonds()/numContexts;
int endIndex = (cc.getContextIndex()+1)*force.getNumBonds()/numContexts;
if (numBonds != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of bonds has changed");
if (numBonds == 0)
return;
// Record the per-bond parameters.
vector > paramVector(numBonds);
vector particles;
for (int i = 0; i < numBonds; i++)
force.getBondParameters(startIndex+i, particles, paramVector[i]);
params->setParameterValues(paramVector, true);
// See if any tabulated functions have changed.
for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
string name = force.getTabulatedFunctionName(i);
if (force.getTabulatedFunction(i).getUpdateCount() != tabulatedFunctionUpdateCount[name]) {
tabulatedFunctionUpdateCount[name] = force.getTabulatedFunction(i).getUpdateCount();
int width;
vector f = cc.getExpressionUtilities().computeFunctionCoefficients(force.getTabulatedFunction(i), width);
tabulatedFunctionArrays[i].upload(f);
}
}
// Mark that the current reordering may be invalid.
cc.invalidateMolecules(info);
}
class CommonCalcCustomCentroidBondForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const CustomCentroidBondForce& force) : force(force) {
}
int getNumParticleGroups() {
return force.getNumBonds();
}
void getParticlesInGroup(int index, vector& particles) {
thread_local static vector parameters;
thread_local static vector groups;
force.getBondParameters(index, groups, parameters);
for (int group : groups) {
vector groupParticles;
vector weights;
force.getGroupParameters(group, groupParticles, weights);
particles.insert(particles.end(), groupParticles.begin(), groupParticles.end());
}
}
bool areGroupsIdentical(int group1, int group2) {
thread_local static vector groups1, groups2;
thread_local static vector parameters1, parameters2;
force.getBondParameters(group1, groups1, parameters1);
force.getBondParameters(group2, groups2, parameters2);
for (int i = 0; i < (int) parameters1.size(); i++)
if (parameters1[i] != parameters2[i])
return false;
for (int i = 0; i < groups1.size(); i++) {
vector groupParticles;
vector weights1, weights2;
force.getGroupParameters(groups1[i], groupParticles, weights1);
force.getGroupParameters(groups2[i], groupParticles, weights2);
if (weights1.size() != weights2.size())
return false;
for (int j = 0; j < weights1.size(); j++)
if (weights1[j] != weights2[j])
return false;
}
return true;
}
private:
const CustomCentroidBondForce& force;
};
CommonCalcCustomCentroidBondForceKernel::~CommonCalcCustomCentroidBondForceKernel() {
ContextSelector selector(cc);
if (params != NULL)
delete params;
}
void CommonCalcCustomCentroidBondForceKernel::initialize(const System& system, const CustomCentroidBondForce& force) {
ContextSelector selector(cc);
numBonds = force.getNumBonds();
if (numBonds == 0)
return;
info = new ForceInfo(force);
cc.addForce(info);
// Record the groups.
numGroups = force.getNumGroups();
vector groupParticleVec;
vector groupWeightVec;
vector groupOffsetVec;
groupOffsetVec.push_back(0);
for (int i = 0; i < numGroups; i++) {
vector particles;
vector weights;
force.getGroupParameters(i, particles, weights);
groupParticleVec.insert(groupParticleVec.end(), particles.begin(), particles.end());
groupOffsetVec.push_back(groupParticleVec.size());
}
vector > normalizedWeights;
CustomCentroidBondForceImpl::computeNormalizedWeights(force, system, normalizedWeights);
for (int i = 0; i < numGroups; i++)
groupWeightVec.insert(groupWeightVec.end(), normalizedWeights[i].begin(), normalizedWeights[i].end());
groupParticles.initialize(cc, groupParticleVec.size(), "groupParticles");
groupParticles.upload(groupParticleVec);
if (cc.getUseDoublePrecision()) {
groupWeights.initialize(cc, groupParticleVec.size(), "groupWeights");
centerPositions.initialize(cc, numGroups, "centerPositions");
}
else {
groupWeights.initialize(cc, groupParticleVec.size(), "groupWeights");
centerPositions.initialize(cc, numGroups, "centerPositions");
}
groupWeights.upload(groupWeightVec, true);
groupOffsets.initialize(cc, groupOffsetVec.size(), "groupOffsets");
groupOffsets.upload(groupOffsetVec);
groupForces.initialize(cc, numGroups*3, "groupForces");
cc.addAutoclearBuffer(groupForces);
// Record the bonds.
int groupsPerBond = force.getNumGroupsPerBond();
vector bondGroupVec(numBonds*groupsPerBond);
params = new ComputeParameterSet(cc, force.getNumPerBondParameters(), numBonds, "customCentroidBondParams", false, cc.getUseDoublePrecision());
vector > paramVector(numBonds);
for (int i = 0; i < numBonds; i++) {
vector groups;
force.getBondParameters(i, groups, paramVector[i]);
for (int j = 0; j < groups.size(); j++)
bondGroupVec[i+j*numBonds] = groups[j];
}
params->setParameterValues(paramVector, true);
bondGroups.initialize(cc, bondGroupVec.size(), "bondGroups");
bondGroups.upload(bondGroupVec);
// Record the tabulated functions.
map functions;
vector > functionDefinitions;
vector functionList;
stringstream extraArgs;
tabulatedFunctionArrays.resize(force.getNumTabulatedFunctions());
for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
functionList.push_back(&force.getTabulatedFunction(i));
string name = force.getTabulatedFunctionName(i);
tabulatedFunctionUpdateCount[name] = force.getTabulatedFunction(i).getUpdateCount();
string arrayName = "table"+cc.intToString(i);
functionDefinitions.push_back(make_pair(name, arrayName));
functions[name] = cc.getExpressionUtilities().getFunctionPlaceholder(force.getTabulatedFunction(i));
int width;
vector f = cc.getExpressionUtilities().computeFunctionCoefficients(force.getTabulatedFunction(i), width);
tabulatedFunctionArrays[i].initialize(cc, f.size(), "TabulatedFunction");
tabulatedFunctionArrays[i].upload(f);
extraArgs << ", GLOBAL const float";
if (width > 1)
extraArgs << width;
extraArgs << "* RESTRICT " << arrayName;
}
// Record information about parameters.
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
globalParamNames[i] = force.getGlobalParameterName(i);
globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
}
map variables;
for (int i = 0; i < groupsPerBond; i++) {
string index = cc.intToString(i+1);
variables["x"+index] = "pos"+index+".x";
variables["y"+index] = "pos"+index+".y";
variables["z"+index] = "pos"+index+".z";
}
for (int i = 0; i < force.getNumPerBondParameters(); i++) {
const string& name = force.getPerBondParameterName(i);
variables[name] = "bondParams"+params->getParameterSuffix(i);
}
needEnergyParamDerivs = (force.getNumEnergyParameterDerivatives() > 0);
if (needEnergyParamDerivs)
extraArgs << ", GLOBAL mixed* RESTRICT energyParamDerivs";
if (force.getNumGlobalParameters() > 0) {
globals.initialize(cc, force.getNumGlobalParameters(), "customCentroidBondGlobals");
globals.upload(globalParamValues);
extraArgs << ", GLOBAL const float* RESTRICT globals";
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = "globals["+cc.intToString(i)+"]";
variables[name] = value;
}
}
// Generate the kernel.
Lepton::ParsedExpression energyExpression = CustomCentroidBondForceImpl::prepareExpression(force, functions);
map forceExpressions;
stringstream compute, initParamDerivs, saveParamDerivs;
for (int i = 0; i < groupsPerBond; i++) {
compute<<"int group"<<(i+1)<<" = bondGroups[index+"<<(i*numBonds)<<"];\n";
compute<<"real4 pos"<<(i+1)<<" = centerPositions[group"<<(i+1)<<"];\n";
}
for (int i = 0; i < (int) params->getParameterInfos().size(); i++) {
ComputeParameterInfo& parameter = params->getParameterInfos()[i];
extraArgs<<", GLOBAL const "<& allParamDerivNames = cc.getEnergyParamDerivNames();
int numDerivs = allParamDerivNames.size();
for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++)
for (int index = 0; index < numDerivs; index++)
if (allParamDerivNames[index] == force.getEnergyParameterDerivativeName(i))
saveParamDerivs << "energyParamDerivs[GLOBAL_ID*" << numDerivs << "+" << index << "] += energyParamDeriv" << i << ";\n";
}
vector forceNames;
for (int i = 0; i < groupsPerBond; i++) {
string istr = cc.intToString(i+1);
string forceName = "force"+istr;
forceNames.push_back(forceName);
compute<<"real3 "< replacements;
replacements["M_PI"] = cc.doubleToString(M_PI);
replacements["NUM_BONDS"] = cc.intToString(numBonds);
replacements["PADDED_NUM_ATOMS"] = cc.intToString(cc.getPaddedNumAtoms());
replacements["EXTRA_ARGS"] = extraArgs.str();
replacements["COMPUTE_FORCE"] = compute.str();
replacements["INIT_PARAM_DERIVS"] = initParamDerivs.str();
replacements["SAVE_PARAM_DERIVS"] = saveParamDerivs.str();
ComputeProgram program = cc.compileProgram(cc.replaceStrings(CommonKernelSources::pointFunctions+CommonKernelSources::customCentroidBond, replacements));
computeCentersKernel = program->createKernel("computeGroupCenters");
computeCentersKernel->addArg(numGroups);
computeCentersKernel->addArg(cc.getPosq());
computeCentersKernel->addArg(groupParticles);
computeCentersKernel->addArg(groupWeights);
computeCentersKernel->addArg(groupOffsets);
computeCentersKernel->addArg(centerPositions);
groupForcesKernel = program->createKernel("computeGroupForces");
groupForcesKernel->addArg(numGroups);
groupForcesKernel->addArg(groupForces);
groupForcesKernel->addArg(); // Energy buffer hasn't been created yet
groupForcesKernel->addArg(centerPositions);
groupForcesKernel->addArg(bondGroups);
for (int i = 0; i < 5; i++)
groupForcesKernel->addArg(); // Periodic box information will be set just before it is executed.
if (needEnergyParamDerivs)
groupForcesKernel->addArg(); // Deriv buffer hasn't been created yet.
for (auto& function : tabulatedFunctionArrays)
groupForcesKernel->addArg(function);
if (globals.isInitialized())
groupForcesKernel->addArg(globals);
for (auto& parameter : params->getParameterInfos())
groupForcesKernel->addArg(parameter.getArray());
applyForcesKernel = program->createKernel("applyForcesToAtoms");
applyForcesKernel->addArg(numGroups);
applyForcesKernel->addArg(groupParticles);
applyForcesKernel->addArg(groupWeights);
applyForcesKernel->addArg(groupOffsets);
applyForcesKernel->addArg(groupForces);
applyForcesKernel->addArg();
}
double CommonCalcCustomCentroidBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
if (numBonds == 0)
return 0.0;
ContextSelector selector(cc);
if (globals.isInitialized()) {
bool changed = false;
for (int i = 0; i < (int) globalParamNames.size(); i++) {
float value = (float) context.getParameter(globalParamNames[i]);
if (value != globalParamValues[i])
changed = true;
globalParamValues[i] = value;
}
if (changed)
globals.upload(globalParamValues);
}
computeCentersKernel->execute(32*numGroups);
groupForcesKernel->setArg(2, cc.getEnergyBuffer());
setPeriodicBoxArgs(cc, groupForcesKernel, 5);
if (needEnergyParamDerivs)
groupForcesKernel->setArg(10, cc.getEnergyParamDerivBuffer());
groupForcesKernel->execute(numBonds);
applyForcesKernel->setArg(5, cc.getLongForceBuffer());
applyForcesKernel->execute(32*numGroups);
return 0.0;
}
void CommonCalcCustomCentroidBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force) {
ContextSelector selector(cc);
if (numBonds != force.getNumBonds())
throw OpenMMException("updateParametersInContext: The number of bonds has changed");
if (numBonds == 0)
return;
// Record the per-bond parameters.
vector > paramVector(numBonds);
vector particles;
for (int i = 0; i < numBonds; i++)
force.getBondParameters(i, particles, paramVector[i]);
params->setParameterValues(paramVector, true);
// See if any tabulated functions have changed.
for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
string name = force.getTabulatedFunctionName(i);
if (force.getTabulatedFunction(i).getUpdateCount() != tabulatedFunctionUpdateCount[name]) {
tabulatedFunctionUpdateCount[name] = force.getTabulatedFunction(i).getUpdateCount();
int width;
vector f = cc.getExpressionUtilities().computeFunctionCoefficients(force.getTabulatedFunction(i), width);
tabulatedFunctionArrays[i].upload(f);
}
}
// Mark that the current reordering may be invalid.
cc.invalidateMolecules(info);
}
class CommonCalcCustomNonbondedForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const CustomNonbondedForce& force) : force(force) {
if (force.getNumInteractionGroups() > 0) {
groupsForParticle.resize(force.getNumParticles());
for (int i = 0; i < force.getNumInteractionGroups(); i++) {
set set1, set2;
force.getInteractionGroupParameters(i, set1, set2);
for (int p : set1)
groupsForParticle[p].insert(2*i);
for (int p : set2)
groupsForParticle[p].insert(2*i+1);
}
}
}
bool areParticlesIdentical(int particle1, int particle2) {
thread_local static vector params1, params2;
force.getParticleParameters(particle1, params1);
force.getParticleParameters(particle2, params2);
for (int i = 0; i < (int) params1.size(); i++)
if (params1[i] != params2[i])
return false;
if (groupsForParticle.size() > 0 && groupsForParticle[particle1] != groupsForParticle[particle2])
return false;
return true;
}
int getNumParticleGroups() {
return force.getNumExclusions();
}
void getParticlesInGroup(int index, vector& particles) {
int particle1, particle2;
force.getExclusionParticles(index, particle1, particle2);
particles.resize(2);
particles[0] = particle1;
particles[1] = particle2;
}
bool areGroupsIdentical(int group1, int group2) {
return true;
}
private:
const CustomNonbondedForce& force;
vector > groupsForParticle;
};
class CommonCalcCustomNonbondedForceKernel::LongRangePostComputation : public ComputeContext::ForcePostComputation {
public:
LongRangePostComputation(ComputeContext& cc, double& longRangeCoefficient, vector& longRangeCoefficientDerivs, CustomNonbondedForce* force) :
cc(cc), longRangeCoefficient(longRangeCoefficient), longRangeCoefficientDerivs(longRangeCoefficientDerivs), force(force) {
}
double computeForceAndEnergy(bool includeForces, bool includeEnergy, int groups) {
if ((groups&(1<getForceGroup())) == 0)
return 0;
if (!cc.getWorkThread().isCurrentThread())
cc.getWorkThread().flush();
Vec3 a, b, c;
cc.getPeriodicBoxVectors(a, b, c);
double volume = a[0]*b[1]*c[2];
map& derivs = cc.getEnergyParamDerivWorkspace();
for (int i = 0; i < longRangeCoefficientDerivs.size(); i++)
derivs[force->getEnergyParameterDerivativeName(i)] += longRangeCoefficientDerivs[i]/volume;
return longRangeCoefficient/volume;
}
private:
ComputeContext& cc;
double& longRangeCoefficient;
vector& longRangeCoefficientDerivs;
CustomNonbondedForce* force;
};
class CommonCalcCustomNonbondedForceKernel::LongRangeTask : public ComputeContext::WorkTask {
public:
LongRangeTask(ComputeContext& cc, Context& context, CustomNonbondedForceImpl::LongRangeCorrectionData& data,
double& longRangeCoefficient, vector& longRangeCoefficientDerivs, CustomNonbondedForce* force) :
cc(cc), context(context), data(data), longRangeCoefficient(longRangeCoefficient),
longRangeCoefficientDerivs(longRangeCoefficientDerivs), force(force) {
}
void execute() {
CustomNonbondedForceImpl::calcLongRangeCorrection(*force, data, context, longRangeCoefficient, longRangeCoefficientDerivs, cc.getThreadPool());
}
private:
ComputeContext& cc;
Context& context;
CustomNonbondedForceImpl::LongRangeCorrectionData& data;
double& longRangeCoefficient;
vector& longRangeCoefficientDerivs;
CustomNonbondedForce* force;
};
CommonCalcCustomNonbondedForceKernel::~CommonCalcCustomNonbondedForceKernel() {
ContextSelector selector(cc);
if (params != NULL)
delete params;
if (computedValues != NULL)
delete computedValues;
if (forceCopy != NULL)
delete forceCopy;
}
void CommonCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
ContextSelector selector(cc);
int forceIndex;
for (forceIndex = 0; forceIndex < system.getNumForces() && &system.getForce(forceIndex) != &force; ++forceIndex)
;
string prefix = (force.getNumInteractionGroups() == 0 ? "custom"+cc.intToString(forceIndex)+"_" : "");
// Record parameters and exclusions.
int numParticles = force.getNumParticles();
int paddedNumParticles = cc.getPaddedNumAtoms();
int numParams = force.getNumPerParticleParameters();
params = new ComputeParameterSet(cc, numParams, paddedNumParticles, "customNonbondedParameters", true);
if (force.getNumGlobalParameters() > 0)
globals.initialize(cc, force.getNumGlobalParameters(), "customNonbondedGlobals");
vector > paramVector(paddedNumParticles, vector(numParams, 0));
vector > exclusionList(numParticles);
for (int i = 0; i < numParticles; i++) {
vector parameters;
force.getParticleParameters(i, parameters);
paramVector[i].resize(parameters.size());
for (int j = 0; j < (int) parameters.size(); j++)
paramVector[i][j] = (float) parameters[j];
exclusionList[i].push_back(i);
}
for (int i = 0; i < force.getNumExclusions(); i++) {
int particle1, particle2;
force.getExclusionParticles(i, particle1, particle2);
exclusionList[particle1].push_back(particle2);
exclusionList[particle2].push_back(particle1);
}
params->setParameterValues(paramVector);
// Record the tabulated functions.
map functions;
vector > functionDefinitions;
vector functionList;
vector tableTypes;
stringstream tableArgs;
tabulatedFunctionArrays.resize(force.getNumTabulatedFunctions());
for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
functionList.push_back(&force.getTabulatedFunction(i));
string name = force.getTabulatedFunctionName(i);
tabulatedFunctionUpdateCount[name] = force.getTabulatedFunction(i).getUpdateCount();
string arrayName = prefix+"table"+cc.intToString(i);
functionDefinitions.push_back(make_pair(name, arrayName));
functions[name] = cc.getExpressionUtilities().getFunctionPlaceholder(force.getTabulatedFunction(i));
int width;
vector f = cc.getExpressionUtilities().computeFunctionCoefficients(force.getTabulatedFunction(i), width);
tabulatedFunctionArrays[i].initialize(cc, f.size(), "TabulatedFunction");
tabulatedFunctionArrays[i].upload(f);
if (force.getNumInteractionGroups() == 0)
cc.getNonbondedUtilities().addArgument(ComputeParameterInfo(tabulatedFunctionArrays[i], arrayName, "float", width));
if (width == 1)
tableTypes.push_back("float");
else
tableTypes.push_back("float"+cc.intToString(width));
tableArgs << ", GLOBAL const float";
if (width > 1)
tableArgs << width;
tableArgs << "* RESTRICT " << arrayName;
}
// Record information for the expressions.
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
globalParamNames[i] = force.getGlobalParameterName(i);
globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
}
if (globals.isInitialized())
globals.upload(globalParamValues);
bool useCutoff = (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff);
bool usePeriodic = (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff && force.getNonbondedMethod() != CustomNonbondedForce::CutoffNonPeriodic);
Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
Lepton::ParsedExpression forceExpression = energyExpression.differentiate("r").optimize();
map forceExpressions;
forceExpressions["real customEnergy = "] = energyExpression;
forceExpressions["tempForce -= "] = forceExpression;
// Record which per-particle parameters and computed values appear in the energy expression.
if (force.getNumComputedValues() > 0)
computedValues = new ComputeParameterSet(cc, force.getNumComputedValues(), paddedNumParticles, "customNonbondedComputedValues", true);
for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
string name = force.getPerParticleParameterName(i);
if (usesVariable(energyExpression, name+"1") || usesVariable(energyExpression, name+"2")) {
paramNames.push_back(name);
paramBuffers.push_back(params->getParameterInfos()[i]);
}
}
for (int i = 0; i < force.getNumComputedValues(); i++) {
string name, expression;
force.getComputedValueParameters(i, name, expression);
if (usesVariable(energyExpression, name+"1") || usesVariable(energyExpression, name+"2")) {
computedValueNames.push_back(name);
computedValueBuffers.push_back(computedValues->getParameterInfos()[i]);
}
}
// Create the kernels.
vector > variables;
ExpressionTreeNode rnode(new Operation::Variable("r"));
variables.push_back(make_pair(rnode, "r"));
variables.push_back(make_pair(ExpressionTreeNode(new Operation::Square(), rnode), "r2"));
variables.push_back(make_pair(ExpressionTreeNode(new Operation::Reciprocal(), rnode), "invR"));
for (int i = 0; i < paramNames.size(); i++) {
variables.push_back(makeVariable(paramNames[i]+"1", "((real) "+prefix+"params"+cc.intToString(i+1)+"1)"));
variables.push_back(makeVariable(paramNames[i]+"2", "((real) "+prefix+"params"+cc.intToString(i+1)+"2)"));
}
for (int i = 0; i < computedValueNames.size(); i++) {
variables.push_back(makeVariable(computedValueNames[i]+"1", prefix+"values"+cc.intToString(i+1)+"1"));
variables.push_back(makeVariable(computedValueNames[i]+"2", prefix+"values"+cc.intToString(i+1)+"2"));
}
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = "globals["+cc.intToString(i)+"]";
variables.push_back(makeVariable(name, prefix+value));
}
for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
string paramName = force.getEnergyParameterDerivativeName(i);
string derivVariable = cc.getNonbondedUtilities().addEnergyParameterDerivative(paramName);
Lepton::ParsedExpression derivExpression = energyExpression.differentiate(paramName).optimize();
forceExpressions[derivVariable+" += interactionScale*switchValue*"] = derivExpression;
}
stringstream compute;
compute << cc.getExpressionUtilities().createExpressions(forceExpressions, variables, functionList, functionDefinitions, prefix+"temp");
map replacements;
replacements["COMPUTE_FORCE"] = compute.str();
replacements["USE_SWITCH"] = (useCutoff && force.getUseSwitchingFunction() ? "1" : "0");
if (force.getUseSwitchingFunction()) {
// Compute the switching coefficients.
replacements["SWITCH_CUTOFF"] = cc.doubleToString(force.getSwitchingDistance());
replacements["SWITCH_C3"] = cc.doubleToString(10/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 3.0));
replacements["SWITCH_C4"] = cc.doubleToString(15/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 4.0));
replacements["SWITCH_C5"] = cc.doubleToString(6/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 5.0));
}
string source = cc.replaceStrings(CommonKernelSources::customNonbonded, replacements);
if (force.getNumInteractionGroups() > 0)
initInteractionGroups(force, source, tableTypes);
else {
cc.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup(), numParticles > 2000);
for (int i = 0; i < paramBuffers.size(); i++)
cc.getNonbondedUtilities().addParameter(ComputeParameterInfo(paramBuffers[i].getArray(), prefix+"params"+cc.intToString(i+1),
paramBuffers[i].getComponentType(), paramBuffers[i].getNumComponents()));
for (int i = 0; i < computedValueBuffers.size(); i++)
cc.getNonbondedUtilities().addParameter(ComputeParameterInfo(computedValueBuffers[i].getArray(), prefix+"values"+cc.intToString(i+1),
computedValueBuffers[i].getComponentType(), computedValueBuffers[i].getNumComponents()));
if (globals.isInitialized()) {
globals.upload(globalParamValues);
cc.getNonbondedUtilities().addArgument(ComputeParameterInfo(globals, prefix+"globals", "float", 1));
}
}
if (force.getNumComputedValues() > 0) {
// Create the kernel to calculate computed values.
stringstream valuesSource, args;
for (int i = 0; i < computedValues->getParameterInfos().size(); i++) {
ComputeParameterInfo& buffer = computedValues->getParameterInfos()[i];
string valueName = "values"+cc.intToString(i+1);
if (i > 0)
args << ", ";
args << "GLOBAL " << buffer.getType() << "* RESTRICT global_" << valueName;
valuesSource << buffer.getType() << " local_" << valueName << ";\n";
}
if (force.getNumGlobalParameters() > 0)
args << ", GLOBAL const float* globals";
for (int i = 0; i < params->getParameterInfos().size(); i++) {
ComputeParameterInfo& buffer = params->getParameterInfos()[i];
string paramName = "params"+cc.intToString(i+1);
args << ", GLOBAL const " << buffer.getType() << "* RESTRICT " << paramName;
}
map variables;
for (int i = 0; i < force.getNumPerParticleParameters(); i++)
variables[force.getPerParticleParameterName(i)] = "params"+params->getParameterSuffix(i, "[index]");
for (int i = 0; i < force.getNumGlobalParameters(); i++)
variables[force.getGlobalParameterName(i)] = "globals["+cc.intToString(i)+"]";
for (int i = 0; i < force.getNumComputedValues(); i++) {
string name, expression;
force.getComputedValueParameters(i, name, expression);
variables[name] = "local_values"+computedValues->getParameterSuffix(i);
map valueExpressions;
valueExpressions["local_values"+computedValues->getParameterSuffix(i)+" = "] = Lepton::Parser::parse(expression, functions).optimize();
valuesSource << cc.getExpressionUtilities().createExpressions(valueExpressions, variables, functionList, functionDefinitions, "value"+cc.intToString(i)+"_temp");
}
for (int i = 0; i < (int) computedValues->getParameterInfos().size(); i++) {
string valueName = "values"+cc.intToString(i+1);
valuesSource << "global_" << valueName << "[index] = local_" << valueName << ";\n";
}
map replacements;
replacements["PARAMETER_ARGUMENTS"] = args.str()+tableArgs.str();
replacements["COMPUTE_VALUES"] = valuesSource.str();
map defines;
defines["NUM_ATOMS"] = cc.intToString(cc.getNumAtoms());
ComputeProgram program = cc.compileProgram(cc.replaceStrings(CommonKernelSources::customNonbondedComputedValues, replacements), defines);
computedValuesKernel = program->createKernel("computePerParticleValues");
for (auto& value : computedValues->getParameterInfos())
computedValuesKernel->addArg(value.getArray());
if (globals.isInitialized())
computedValuesKernel->addArg(globals);
for (auto& parameter : params->getParameterInfos())
computedValuesKernel->addArg(parameter.getArray());
for (auto& function : tabulatedFunctionArrays)
computedValuesKernel->addArg(function);
}
info = new ForceInfo(force);
cc.addForce(info);
// Record information for the long range correction.
if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic && force.getUseLongRangeCorrection() && cc.getContextIndex() == 0) {
forceCopy = new CustomNonbondedForce(force);
longRangeCorrectionData = CustomNonbondedForceImpl::prepareLongRangeCorrection(force, cc.getThreadPool().getNumThreads());
cc.addPostComputation(new LongRangePostComputation(cc, longRangeCoefficient, longRangeCoefficientDerivs, forceCopy));
hasInitializedLongRangeCorrection = false;
}
else {
longRangeCoefficient = 0.0;
hasInitializedLongRangeCorrection = true;
}
}
void CommonCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbondedForce& force, const string& interactionSource, const vector& tableTypes) {
// Process groups to form tiles.
vector > atomLists;
vector > tiles;
vector tileGroup;
vector > duplicateAtomsForGroup;
for (int group = 0; group < force.getNumInteractionGroups(); group++) {
// Get the list of atoms in this group and sort them.
set set1, set2;
force.getInteractionGroupParameters(group, set1, set2);
vector atoms1, atoms2;
atoms1.insert(atoms1.begin(), set1.begin(), set1.end());
atoms2.insert(atoms2.begin(), set2.begin(), set2.end());
sort(atoms1.begin(), atoms1.end());
sort(atoms2.begin(), atoms2.end());
duplicateAtomsForGroup.push_back(vector());
set_intersection(set1.begin(), set1.end(), set2.begin(), set2.end(),
inserter(duplicateAtomsForGroup[group], duplicateAtomsForGroup[group].begin()));
sort(duplicateAtomsForGroup[group].begin(), duplicateAtomsForGroup[group].end());
// Find how many tiles we will create for this group.
int tileWidth = min(min(32, (int) atoms1.size()), (int) atoms2.size());
if (tileWidth == 0)
continue;
int numBlocks1 = (atoms1.size()+tileWidth-1)/tileWidth;
int numBlocks2 = (atoms2.size()+tileWidth-1)/tileWidth;
// Add the tiles.
int firstTile = tiles.size();
for (int i = 0; i < numBlocks1; i++)
for (int j = 0; j < numBlocks2; j++) {
tiles.push_back(make_pair(atomLists.size()+i, atomLists.size()+numBlocks1+j));
tileGroup.push_back(group);
}
// Add the atom lists.
for (int i = 0; i < numBlocks1; i++) {
vector atoms;
int first = i*tileWidth;
int last = min((i+1)*tileWidth, (int) atoms1.size());
for (int j = first; j < last; j++)
atoms.push_back(atoms1[j]);
atomLists.push_back(atoms);
}
for (int i = 0; i < numBlocks2; i++) {
vector atoms;
int first = i*tileWidth;
int last = min((i+1)*tileWidth, (int) atoms2.size());
for (int j = first; j < last; j++)
atoms.push_back(atoms2[j]);
atomLists.push_back(atoms);
}
}
// Build a lookup table for quickly identifying excluded interactions.
vector > exclusions(force.getNumParticles());
for (int i = 0; i < force.getNumExclusions(); i++) {
int p1, p2;
force.getExclusionParticles(i, p1, p2);
exclusions[p1].insert(p2);
exclusions[p2].insert(p1);
}
// Build the exclusion flags for each tile. While we're at it, filter out tiles
// where all interactions are excluded, and sort the tiles by size.
vector > exclusionFlags(tiles.size());
vector > tileOrder;
for (int tile = 0; tile < tiles.size(); tile++) {
bool swapped = false;
if (atomLists[tiles[tile].first].size() < atomLists[tiles[tile].second].size()) {
// For efficiency, we want the first axis to be the larger one.
int swap = tiles[tile].first;
tiles[tile].first = tiles[tile].second;
tiles[tile].second = swap;
swapped = true;
}
vector& atoms1 = atomLists[tiles[tile].first];
vector& atoms2 = atomLists[tiles[tile].second];
vector& duplicateAtoms = duplicateAtomsForGroup[tileGroup[tile]];
vector& flags = exclusionFlags[tile];
flags.resize(atoms1.size(), (int) (1LL< a2) == swapped && a1IsDuplicate && binary_search(duplicateAtoms.begin(), duplicateAtoms.end(), a2))
isExcluded = true; // Both atoms are in both sets, so skip duplicate interactions.
if (isExcluded) {
flags[i] &= -1-(1< tileSetStart;
tileSetStart.push_back(0);
int tileSetSize = 0;
for (int i = 0; i < tileOrder.size(); i++) {
int tile = tileOrder[i].second;
int size = atomLists[tiles[tile].first].size();
if (tileSetSize+size > 32) {
tileSetStart.push_back(i);
tileSetSize = 0;
}
tileSetSize += size;
}
tileSetStart.push_back(tileOrder.size());
// Build the data structures.
int numTileSets = tileSetStart.size()-1;
vector groupData;
for (int tileSet = 0; tileSet < numTileSets; tileSet++) {
int indexInTileSet = 0;
int minSize = 0;
if (cc.getSIMDWidth() < 32) {
// We need to include a barrier inside the inner loop, so ensure that all
// threads will loop the same number of times.
for (int i = tileSetStart[tileSet]; i < tileSetStart[tileSet+1]; i++)
minSize = max(minSize, (int) atomLists[tiles[tileOrder[i].second].first].size());
}
for (int i = tileSetStart[tileSet]; i < tileSetStart[tileSet+1]; i++) {
int tile = tileOrder[i].second;
vector& atoms1 = atomLists[tiles[tile].first];
vector& atoms2 = atomLists[tiles[tile].second];
int range = indexInTileSet + ((indexInTileSet+max(minSize, (int) atoms1.size()))<<16);
int allFlags = (1< 0 ? exclusionFlags[tile][j] : allFlags);
groupData.push_back(mm_int4(a1, a2, range, flags<(cc, groupData.size(), "interactionGroupData");
interactionGroupData.upload(groupData);
numGroupTiles.initialize(cc, 1, "numGroupTiles");
// Allocate space for a neighbor list, if necessary.
if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff && groupData.size() > cc.getNumThreadBlocks()) {
filteredGroupData.initialize(cc, groupData.size(), "filteredGroupData");
interactionGroupData.copyTo(filteredGroupData);
int numTiles = groupData.size()/32;
numGroupTiles.upload(&numTiles);
}
// Create the kernel.
hasParamDerivs = (force.getNumEnergyParameterDerivatives() > 0);
map replacements;
replacements["COMPUTE_INTERACTION"] = interactionSource;
const string suffixes[] = {"x", "y", "z", "w"};
stringstream localData;
int localDataSize = 0;
for (int i = 0; i < paramBuffers.size(); i++) {
localData<& allParamDerivNames = cc.getEnergyParamDerivNames();
int numDerivs = allParamDerivNames.size();
for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
string paramName = force.getEnergyParameterDerivativeName(i);
string derivVariable = cc.getNonbondedUtilities().addEnergyParameterDerivative(paramName);
initDerivs<<"mixed "< defines;
if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
defines["USE_CUTOFF"] = "1";
if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic)
defines["USE_PERIODIC"] = "1";
int localMemorySize = max(32, cc.getNonbondedUtilities().getForceThreadBlockSize());
defines["LOCAL_MEMORY_SIZE"] = cc.intToString(localMemorySize);
defines["WARPS_IN_BLOCK"] = cc.intToString(localMemorySize/32);
double cutoff = force.getCutoffDistance();
defines["CUTOFF_SQUARED"] = cc.doubleToString(cutoff*cutoff);
double paddedCutoff = cc.getNonbondedUtilities().padCutoff(cutoff);
defines["PADDED_CUTOFF_SQUARED"] = cc.doubleToString(paddedCutoff*paddedCutoff);
defines["PADDED_NUM_ATOMS"] = cc.intToString(cc.getPaddedNumAtoms());
defines["TILE_SIZE"] = "32";
defines["NUM_TILES"] = cc.intToString(numTileSets);
int numContexts = cc.getNumContexts();
int startIndex = cc.getContextIndex()*numTileSets/numContexts;
int endIndex = (cc.getContextIndex()+1)*numTileSets/numContexts;
defines["FIRST_TILE"] = cc.intToString(startIndex);
defines["LAST_TILE"] = cc.intToString(endIndex);
if ((localDataSize/4)%2 == 0 && !cc.getUseDoublePrecision())
defines["PARAMETER_SIZE_IS_EVEN"] = "1";
ComputeProgram program = cc.compileProgram(cc.replaceStrings(CommonKernelSources::customNonbondedGroups, replacements), defines);
interactionGroupKernel = program->createKernel("computeInteractionGroups");
prepareNeighborListKernel = program->createKernel("prepareToBuildNeighborList");
buildNeighborListKernel = program->createKernel("buildNeighborList");
numGroupThreadBlocks = cc.getNonbondedUtilities().getNumForceThreadBlocks();
}
double CommonCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
useNeighborList = (filteredGroupData.isInitialized() && cc.getNonbondedUtilities().getUseCutoff());
if (useNeighborList && cc.getContextIndex() > 0) {
// When using a neighbor list, run the whole calculation on a single device.
return 0.0;
}
ContextSelector selector(cc);
bool recomputeLongRangeCorrection = !hasInitializedLongRangeCorrection;
if (globals.isInitialized()) {
bool changed = false;
for (int i = 0; i < (int) globalParamNames.size(); i++) {
float value = (float) context.getParameter(globalParamNames[i]);
if (value != globalParamValues[i])
changed = true;
globalParamValues[i] = value;
}
if (changed) {
globals.upload(globalParamValues);
if (forceCopy != NULL)
recomputeLongRangeCorrection = true;
}
}
if (recomputeLongRangeCorrection) {
if (includeEnergy || forceCopy->getNumEnergyParameterDerivatives() > 0) {
cc.getWorkThread().addTask(new LongRangeTask(cc, context.getOwner(), longRangeCorrectionData, longRangeCoefficient, longRangeCoefficientDerivs, forceCopy));
hasInitializedLongRangeCorrection = true;
}
else
hasInitializedLongRangeCorrection = false;
}
if (computedValues != NULL)
computedValuesKernel->execute(cc.getNumAtoms());
if (interactionGroupData.isInitialized()) {
if (!hasInitializedKernel) {
hasInitializedKernel = true;
interactionGroupKernel->addArg(cc.getLongForceBuffer());
interactionGroupKernel->addArg(cc.getEnergyBuffer());
interactionGroupKernel->addArg(cc.getPosq());
interactionGroupKernel->addArg((useNeighborList ? filteredGroupData : interactionGroupData));
interactionGroupKernel->addArg(numGroupTiles);
interactionGroupKernel->addArg((int) useNeighborList);
for (int i = 0; i < 5; i++)
interactionGroupKernel->addArg(); // Periodic box information will be set just before it is executed.
interactionGroupKernel->addArg((int) cc.getEnergyParamDerivNames().size());
for (auto& buffer : paramBuffers)
interactionGroupKernel->addArg(buffer.getArray());
for (auto& buffer : computedValueBuffers)
interactionGroupKernel->addArg(buffer.getArray());
for (auto& function : tabulatedFunctionArrays)
interactionGroupKernel->addArg(function);
if (globals.isInitialized())
interactionGroupKernel->addArg(globals);
if (hasParamDerivs)
interactionGroupKernel->addArg(cc.getEnergyParamDerivBuffer());
if (useNeighborList) {
// Initialize kernels for building the interaction group neighbor list.
prepareNeighborListKernel->addArg(cc.getNonbondedUtilities().getRebuildNeighborList());
prepareNeighborListKernel->addArg(numGroupTiles);
buildNeighborListKernel->addArg(cc.getNonbondedUtilities().getRebuildNeighborList());
buildNeighborListKernel->addArg(numGroupTiles);
buildNeighborListKernel->addArg(cc.getPosq());
buildNeighborListKernel->addArg(interactionGroupData);
buildNeighborListKernel->addArg(filteredGroupData);
for (int i = 0; i < 5; i++)
buildNeighborListKernel->addArg(); // Periodic box information will be set just before it is executed.
}
}
int forceThreadBlockSize = max(32, cc.getNonbondedUtilities().getForceThreadBlockSize());
if (useNeighborList) {
// Rebuild the neighbor list, if necessary.
setPeriodicBoxArgs(cc, buildNeighborListKernel, 5);
prepareNeighborListKernel->execute(1, 1);
buildNeighborListKernel->execute(numGroupThreadBlocks*forceThreadBlockSize, forceThreadBlockSize);
}
setPeriodicBoxArgs(cc, interactionGroupKernel, 6);
interactionGroupKernel->execute(numGroupThreadBlocks*forceThreadBlockSize, forceThreadBlockSize);
}
return 0;
}
void CommonCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force) {
ContextSelector selector(cc);
int numParticles = force.getNumParticles();
if (numParticles != cc.getNumAtoms())
throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the per-particle parameters.
int paddedNumParticles = cc.getPaddedNumAtoms();
int numParams = force.getNumPerParticleParameters();
vector > paramVector(paddedNumParticles, vector(numParams, 0));
vector parameters;
for (int i = 0; i < numParticles; i++) {
force.getParticleParameters(i, parameters);
paramVector[i].resize(parameters.size());
for (int j = 0; j < (int) parameters.size(); j++)
paramVector[i][j] = (float) parameters[j];
}
params->setParameterValues(paramVector);
// If necessary, recompute the long range correction.
if (forceCopy != NULL) {
longRangeCorrectionData = CustomNonbondedForceImpl::prepareLongRangeCorrection(force, cc.getThreadPool().getNumThreads());
CustomNonbondedForceImpl::calcLongRangeCorrection(force, longRangeCorrectionData, context.getOwner(), longRangeCoefficient, longRangeCoefficientDerivs, cc.getThreadPool());
hasInitializedLongRangeCorrection = false;
*forceCopy = force;
}
// See if any tabulated functions have changed.
for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
string name = force.getTabulatedFunctionName(i);
if (force.getTabulatedFunction(i).getUpdateCount() != tabulatedFunctionUpdateCount[name]) {
tabulatedFunctionUpdateCount[name] = force.getTabulatedFunction(i).getUpdateCount();
int width;
vector f = cc.getExpressionUtilities().computeFunctionCoefficients(force.getTabulatedFunction(i), width);
tabulatedFunctionArrays[i].upload(f);
}
}
// Mark that the current reordering may be invalid.
cc.invalidateMolecules(info);
}
class CommonCalcGBSAOBCForceKernel::ForceInfo : public ComputeForceInfo {
public:
ForceInfo(const GBSAOBCForce& force) : force(force) {
}
bool areParticlesIdentical(int particle1, int particle2) {
double charge1, charge2, radius1, radius2, scale1, scale2;
force.getParticleParameters(particle1, charge1, radius1, scale1);
force.getParticleParameters(particle2, charge2, radius2, scale2);
return (charge1 == charge2 && radius1 == radius2 && scale1 == scale2);
}
private:
const GBSAOBCForce& force;
};
void CommonCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
ContextSelector selector(cc);
if (cc.getNumContexts() > 1)
throw OpenMMException("GBSAOBCForce does not support using multiple devices");
int forceIndex;
for (forceIndex = 0; forceIndex < system.getNumForces() && &system.getForce(forceIndex) != &force; ++forceIndex)
;
string prefix = "obc"+cc.intToString(forceIndex)+"_";
NonbondedUtilities& nb = cc.getNonbondedUtilities();
params.initialize(cc, cc.getPaddedNumAtoms(), "gbsaObcParams");
int elementSize = (cc.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
charges.initialize(cc, cc.getPaddedNumAtoms(), elementSize, "gbsaObcCharges");
bornRadii.initialize(cc, cc.getPaddedNumAtoms(), elementSize, "bornRadii");
obcChain.initialize(cc, cc.getPaddedNumAtoms(), elementSize, "obcChain");
bornSum.initialize(cc, cc.getPaddedNumAtoms(), "bornSum");
bornForce.initialize(cc, cc.getPaddedNumAtoms(), "bornForce");
cc.addAutoclearBuffer(bornSum);
cc.addAutoclearBuffer(bornForce);
vector chargeVec(cc.getPaddedNumAtoms());
vector paramsVector(cc.getPaddedNumAtoms(), mm_float2(1,1));
const double dielectricOffset = 0.009;
for (int i = 0; i < force.getNumParticles(); i++) {
double charge, radius, scalingFactor;
force.getParticleParameters(i, charge, radius, scalingFactor);
radius -= dielectricOffset;
chargeVec[i] = charge;
paramsVector[i] = mm_float2((float) radius, (float) (scalingFactor*radius));
}
charges.upload(chargeVec, true);
params.upload(paramsVector);
prefactor = -ONE_4PI_EPS0*((1.0/force.getSoluteDielectric())-(1.0/force.getSolventDielectric()));
surfaceAreaFactor = -6.0*4*M_PI*force.getSurfaceAreaEnergy();
bool useCutoff = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff);
bool usePeriodic = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff && force.getNonbondedMethod() != GBSAOBCForce::CutoffNonPeriodic);
cutoff = force.getCutoffDistance();
string source = CommonKernelSources::gbsaObc2;
map replacements;
replacements["CHARGE1"] = prefix+"charge1";
replacements["CHARGE2"] = prefix+"charge2";
replacements["OBC_PARAMS1"] = prefix+"obcParams1";
replacements["OBC_PARAMS2"] = prefix+"obcParams2";
replacements["BORN_FORCE1"] = prefix+"bornForce1";
replacements["BORN_FORCE2"] = prefix+"bornForce2";
source = cc.replaceStrings(source, replacements);
nb.addInteraction(useCutoff, usePeriodic, false, cutoff, vector >(), source, force.getForceGroup());
nb.addParameter(ComputeParameterInfo(charges, prefix+"charge", "float", 1));
nb.addParameter(ComputeParameterInfo(params, prefix+"obcParams", "float", 2));
nb.addParameter(ComputeParameterInfo(bornForce, prefix+"bornForce", "mm_long", 1));
info = new ForceInfo(force);
cc.addForce(info);
}
double CommonCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
ContextSelector selector(cc);
NonbondedUtilities& nb = cc.getNonbondedUtilities();
bool deviceIsCpu = cc.getIsCPU();
if (!hasCreatedKernels) {
// These Kernels cannot be created in initialize(), because the NonbondedUtilities has not been initialized yet then.
hasCreatedKernels = true;
maxTiles = (nb.getUseCutoff() ? nb.getInteractingTiles().getSize() : 0);
int numAtomBlocks = cc.getPaddedNumAtoms()/32;
map defines;
if (nb.getUseCutoff())
defines["USE_CUTOFF"] = "1";
if (nb.getUsePeriodic())
defines["USE_PERIODIC"] = "1";
defines["CUTOFF_SQUARED"] = cc.doubleToString(cutoff*cutoff);
defines["CUTOFF"] = cc.doubleToString(cutoff);
defines["PREFACTOR"] = cc.doubleToString(prefactor);
defines["SURFACE_AREA_FACTOR"] = cc.doubleToString(surfaceAreaFactor);
defines["NUM_ATOMS"] = cc.intToString(cc.getNumAtoms());
defines["PADDED_NUM_ATOMS"] = cc.intToString(cc.getPaddedNumAtoms());
defines["NUM_BLOCKS"] = cc.intToString(numAtomBlocks);
defines["FORCE_WORK_GROUP_SIZE"] = cc.intToString(nb.getForceThreadBlockSize());
defines["TILE_SIZE"] = "32";
int numExclusionTiles = nb.getExclusionTiles().getSize();
defines["NUM_TILES_WITH_EXCLUSIONS"] = cc.intToString(numExclusionTiles);
defines["FIRST_EXCLUSION_TILE"] = "0";
defines["LAST_EXCLUSION_TILE"] = cc.intToString(numExclusionTiles);
string file;
if (deviceIsCpu)
file = CommonKernelSources::gbsaObc_cpu;
else
file = CommonKernelSources::gbsaObc;
ComputeProgram program = cc.compileProgram(file, defines);
computeBornSumKernel = program->createKernel("computeBornSum");
computeBornSumKernel->addArg(bornSum);
computeBornSumKernel->addArg(cc.getPosq());
computeBornSumKernel->addArg(charges);
computeBornSumKernel->addArg(params);
if (nb.getUseCutoff()) {
computeBornSumKernel->addArg(nb.getInteractingTiles());
computeBornSumKernel->addArg(nb.getInteractionCount());
for (int i = 0; i < 5; i++)
computeBornSumKernel->addArg(); // The periodic box size arguments are set when the kernel is executed.
computeBornSumKernel->addArg(maxTiles);
computeBornSumKernel->addArg(nb.getBlockCenters());
computeBornSumKernel->addArg(nb.getBlockBoundingBoxes());
computeBornSumKernel->addArg(nb.getInteractingAtoms());
}
else
computeBornSumKernel->addArg(numAtomBlocks*(numAtomBlocks+1)/2);
computeBornSumKernel->addArg(nb.getExclusionTiles());
force1Kernel = program->createKernel("computeGBSAForce1");
force1Kernel->addArg(cc.getLongForceBuffer());
force1Kernel->addArg(bornForce);
force1Kernel->addArg(cc.getEnergyBuffer());
force1Kernel->addArg(cc.getPosq());
force1Kernel->addArg(charges);
force1Kernel->addArg(bornRadii);
force1Kernel->addArg(); // Whether to include energy.
if (nb.getUseCutoff()) {
force1Kernel->addArg(nb.getInteractingTiles());
force1Kernel->addArg(nb.getInteractionCount());
for (int i = 0; i < 5; i++)
force1Kernel->addArg(); // The periodic box size arguments are set when the kernel is executed.
force1Kernel->addArg(maxTiles);
force1Kernel->addArg(nb.getBlockCenters());
force1Kernel->addArg(nb.getBlockBoundingBoxes());
force1Kernel->addArg(nb.getInteractingAtoms());
}
else
force1Kernel->addArg(numAtomBlocks*(numAtomBlocks+1)/2);
force1Kernel->addArg(nb.getExclusionTiles());
program = cc.compileProgram(CommonKernelSources::gbsaObcReductions, defines);
reduceBornSumKernel = program->createKernel("reduceBornSum");
reduceBornSumKernel->addArg(1.0f);
reduceBornSumKernel->addArg(0.8f);
reduceBornSumKernel->addArg(4.85f);
reduceBornSumKernel->addArg(bornSum);
reduceBornSumKernel->addArg(params);
reduceBornSumKernel->addArg(bornRadii);
reduceBornSumKernel->addArg(obcChain);
reduceBornForceKernel = program->createKernel("reduceBornForce");
reduceBornForceKernel->addArg(bornForce);
reduceBornForceKernel->addArg(cc.getEnergyBuffer());
reduceBornForceKernel->addArg(params);
reduceBornForceKernel->addArg(bornRadii);
reduceBornForceKernel->addArg(obcChain);
}
force1Kernel->setArg(6, (int) includeEnergy);
if (nb.getUseCutoff()) {
setPeriodicBoxArgs(cc, computeBornSumKernel, 6);
setPeriodicBoxArgs(cc, force1Kernel, 9);
if (maxTiles < nb.getInteractingTiles().getSize()) {
maxTiles = nb.getInteractingTiles().getSize();
computeBornSumKernel->setArg(11, maxTiles);
force1Kernel->setArg(14, maxTiles);
}
}
computeBornSumKernel->execute(nb.getNumForceThreadBlocks()*nb.getForceThreadBlockSize(), nb.getForceThreadBlockSize());
reduceBornSumKernel->execute(cc.getPaddedNumAtoms());
force1Kernel->execute(nb.getNumForceThreadBlocks()*nb.getForceThreadBlockSize(), nb.getForceThreadBlockSize());
reduceBornForceKernel->execute(cc.getPaddedNumAtoms());
return 0.0;
}
void CommonCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force) {
// Make sure the new parameters are acceptable.
ContextSelector selector(cc);
int numParticles = force.getNumParticles();
if (numParticles != cc.getNumAtoms())
throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the per-particle parameters.
vector chargeVector(cc.getPaddedNumAtoms(), 0.0);
vector