Skip to content

GitLab

Commit 6c1d5eb1 authored by peastman's avatar peastman
Browse files

Created CPU implementation of GayBerneForce

parent e45f1f28
Hide whitespace changes
Inline Side-by-side
Showing with 675 additions and 15 deletions
platforms/cpu/include/CpuGayBerneForce.h 0 → 100644
View file @ 6c1d5eb1
/* -------------------------------------------------------------------------- *
* 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) 2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#ifndef OPENMM_CPU_GAYBERNEFORCE_H__
#define OPENMM_CPU_GAYBERNEFORCE_H__
#include "openmm/GayBerneForce.h"
#include "CpuNeighborList.h"
#include "CpuPlatform.h"
#include "RealVec.h"
#include <set>
#include <utility>
namespace OpenMM {
class OPENMM_EXPORT CpuGayBerneForce {
public:
struct Matrix;
/**
* Constructor.
*/
CpuGayBerneForce(const GayBerneForce& force);
/**
* Destructor.
*/
~CpuGayBerneForce();
/**
* Compute the interaction.
*
* @param positions the positions of the atoms
* @param forces forces will be added to this vector
* @param boxVectors the periodic box vectors
* @param data the platform data for the current context
* @return the energy of the interaction
*/
RealOpenMM calculateForce(const std::vector<RealVec>& positions, std::vector<RealVec>& forces, const RealVec* boxVectors, CpuPlatform::PlatformData& data);
private:
struct ParticleInfo;
struct ExceptionInfo;
std::vector<ParticleInfo> particles;
std::vector<ExceptionInfo> exceptions;
std::set<std::pair<int, int> > exclusions;
std::vector<std::set<int> > particleExclusions;
GayBerneForce::NonbondedMethod nonbondedMethod;
RealOpenMM cutoffDistance, switchingDistance;
bool useSwitchingFunction;
std::vector<RealOpenMM> s;
std::vector<Matrix> A, B, G;
CpuNeighborList* neighborList;
void computeEllipsoidFrames(const std::vector<RealVec>& positions);
void applyTorques(const std::vector<RealVec>& positions, std::vector<RealVec>& forces, const std::vector<RealVec>& torques);
RealOpenMM computeOneInteraction(int particle1, int particle2, RealOpenMM sigma, RealOpenMM epsilon, const std::vector<RealVec>& positions,
std::vector<RealVec>& forces, std::vector<RealVec>& torques, const RealVec* boxVectors);
};
struct CpuGayBerneForce::ParticleInfo {
int xparticle, yparticle;
RealOpenMM sigmaOver2, sqrtEpsilon, rx, ry, rz, ex, ey, ez;
bool isPointParticle;
};
struct CpuGayBerneForce::ExceptionInfo {
int particle1, particle2;
RealOpenMM sigma, epsilon;
};
struct CpuGayBerneForce::Matrix {
RealOpenMM v[3][3];
RealVec operator*(const RealVec& r) {
return RealVec(v[0][0]*r[0] + v[0][1]*r[1] + v[0][2]*r[2],
v[1][0]*r[0] + v[1][1]*r[1] + v[1][2]*r[2],
v[2][0]*r[0] + v[2][1]*r[1] + v[2][2]*r[2]);
}
Matrix operator+(const Matrix& m) {
Matrix result;
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
result.v[i][j] = v[i][j]+m.v[i][j];
return result;
}
RealOpenMM determinant() {
return (v[0][0]*v[1][1]*v[2][2] + v[0][1]*v[1][2]*v[2][0] + v[0][2]*v[1][0]*v[2][1] -
v[0][0]*v[1][2]*v[2][1] - v[0][1]*v[1][0]*v[2][2] - v[0][2]*v[1][1]*v[2][0]);
}
Matrix inverse() {
RealOpenMM invDet = 1/determinant();
Matrix result;
result.v[0][0] = invDet*(v[1][1]*v[2][2] - v[1][2]*v[2][1]);
result.v[1][0] = -invDet*(v[1][0]*v[2][2] - v[1][2]*v[2][0]);
result.v[2][0] = invDet*(v[1][0]*v[2][1] - v[1][1]*v[2][0]);
result.v[0][1] = -invDet*(v[0][1]*v[2][2] - v[0][2]*v[2][1]);
result.v[1][1] = invDet*(v[0][0]*v[2][2] - v[0][2]*v[2][0]);
result.v[2][1] = -invDet*(v[0][0]*v[2][1] - v[0][1]*v[2][0]);
result.v[0][2] = invDet*(v[0][1]*v[1][2] - v[0][2]*v[1][1]);
result.v[1][2] = -invDet*(v[0][0]*v[1][2] - v[0][2]*v[1][0]);
result.v[2][2] = invDet*(v[0][0]*v[1][1] - v[0][1]*v[1][0]);
return result;
}
};
static RealVec operator*(const RealVec& r, CpuGayBerneForce::Matrix& m) {
return RealVec(m.v[0][0]*r[0] + m.v[1][0]*r[1] + m.v[2][0]*r[2],
m.v[0][1]*r[0] + m.v[1][1]*r[1] + m.v[2][1]*r[2],
m.v[0][2]*r[0] + m.v[1][2]*r[1] + m.v[2][2]*r[2]);
}
} // namespace OpenMM
#endif // OPENMM_CPU_GAYBERNEFORCE_H__
platforms/cpu/include/CpuKernels.h
View file @ 6c1d5eb1
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013 Stanford University and the Authors. *
* Portions copyright (c) 2013-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -36,6 +36,7 @@
#include "CpuCustomGBForce.h"
#include "CpuCustomManyParticleForce.h"
#include "CpuCustomNonbondedForce.h"
#include "CpuGayBerneForce.h"
#include "CpuGBSAOBCForce.h"
#include "CpuLangevinDynamics.h"
#include "CpuNeighborList.h"
......@@ -445,6 +446,42 @@ private:
NonbondedMethod nonbondedMethod;
};
/**
* This kernel is invoked by GayBerneForce to calculate the forces acting on the system.
*/
class CpuCalcGayBerneForceKernel : public CalcGayBerneForceKernel {
public:
CpuCalcGayBerneForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) : CalcGayBerneForceKernel(name, platform),
data(data), ixn(NULL) {
}
~CpuCalcGayBerneForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the GayBerneForce this kernel will be used for
*/
void initialize(const System& system, const GayBerneForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the GayBerneForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const GayBerneForce& force);
private:
CpuPlatform::PlatformData& data;
CpuGayBerneForce* ixn;
};
/**
* This kernel is invoked by LangevinIntegrator to take one time step.
*/
......
platforms/cpu/src/CpuGayBerneForce.cpp 0 → 100644
View file @ 6c1d5eb1
/* -------------------------------------------------------------------------- *
* 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) 2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "CpuGayBerneForce.h"
#include "ReferenceForce.h"
#include "openmm/OpenMMException.h"
#include "openmm/GayBerneForce.h"
#include <cmath>
using namespace OpenMM;
using namespace std;
CpuGayBerneForce::CpuGayBerneForce(const GayBerneForce& force) : neighborList(NULL) {
// Record the force parameters.
int numParticles = force.getNumParticles();
particles.resize(numParticles);
for (int i = 0; i < numParticles; i++) {
ParticleInfo& p = particles[i];
double sigma, epsilon, sx, sy, sz, ex, ey, ez;
force.getParticleParameters(i, sigma, epsilon, p.xparticle, p.yparticle, sx, sy, sz, ex, ey, ez);
p.sigmaOver2 = 0.5*sigma;
p.sqrtEpsilon = sqrt(epsilon);
p.rx = 0.5*sx;
p.ry = 0.5*sy;
p.rz = 0.5*sz;
p.ex = ex;
p.ey = ey;
p.ez = ez;
p.isPointParticle = (sx == sigma && sy == sigma && sz == sigma && ex == 1.0 && ey == 1.0 && ez == 1.0);
}
int numExceptions = force.getNumExceptions();
exceptions.resize(numExceptions);
particleExclusions.resize(numParticles);
for (int i = 0; i < numExceptions; i++) {
ExceptionInfo& e = exceptions[i];
double sigma, epsilon;
force.getExceptionParameters(i, e.particle1, e.particle2, sigma, epsilon);
e.sigma = sigma;
e.epsilon = epsilon;
exclusions.insert(make_pair(min(e.particle1, e.particle2), max(e.particle1, e.particle2)));
particleExclusions[e.particle1].insert(e.particle2);
particleExclusions[e.particle2].insert(e.particle1);
}
nonbondedMethod = force.getNonbondedMethod();
cutoffDistance = force.getCutoffDistance();
switchingDistance = force.getSwitchingDistance();
useSwitchingFunction = force.getUseSwitchingFunction();
if (nonbondedMethod != GayBerneForce::NoCutoff)
neighborList = new CpuNeighborList(4);
// Allocate workspace for calculations.
s.resize(numParticles);
A.resize(numParticles);
B.resize(numParticles);
G.resize(numParticles);
// We can precompute the shape factors.
for (int i = 0; i < numParticles; i++) {
ParticleInfo& p = particles[i];
s[i] = (p.rx*p.ry + p.rz*p.rz)*sqrtf(p.rx*p.ry);
}
}
CpuGayBerneForce::~CpuGayBerneForce() {
if (neighborList != NULL)
delete neighborList;
}
RealOpenMM CpuGayBerneForce::calculateForce(const vector<RealVec>& positions, vector<RealVec>& forces, const RealVec* boxVectors, CpuPlatform::PlatformData& data) {
if (nonbondedMethod == GayBerneForce::CutoffPeriodic) {
double minAllowedSize = 1.999999*cutoffDistance;
if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
}
// Build the neighbor list.
int numParticles = particles.size();
if (nonbondedMethod != GayBerneForce::NoCutoff)
neighborList->computeNeighborList(numParticles, data.posq, particleExclusions, boxVectors, nonbondedMethod == GayBerneForce::CutoffPeriodic, cutoffDistance, data.threads);
// First find the orientations of the particles and compute the matrices we'll be needing.
computeEllipsoidFrames(positions);
// Compute standard interactions.
double energy = 0;
vector<RealVec> torques(numParticles, Vec3());
if (neighborList == NULL) {
for (int i = 1; i < numParticles; i++) {
if (particles[i].sqrtEpsilon == 0.0f)
continue;
for (int j = 0; j < i; j++) {
if (particles[j].sqrtEpsilon == 0.0f)
continue;
if (particleExclusions[i].find(j) != particleExclusions[i].end())
continue; // This interaction will be handled by an exception.
RealOpenMM sigma = particles[i].sigmaOver2+particles[j].sigmaOver2;
RealOpenMM epsilon = particles[i].sqrtEpsilon*particles[j].sqrtEpsilon;
energy += computeOneInteraction(i, j, sigma, epsilon, positions, forces, torques, boxVectors);
}
}
}
else {
for (int blockIndex = 0; blockIndex < neighborList->getNumBlocks(); blockIndex++) {
const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
const vector<char>& exclusions = neighborList->getBlockExclusions(blockIndex);
for (int i = 0; i < (int) neighbors.size(); i++) {
int first = neighbors[i];
if (particles[first].sqrtEpsilon == 0.0f)
continue;
for (int k = 0; k < 4; k++) {
if ((exclusions[i] & (1<<k)) == 0) {
int second = blockAtom[k];
if (particles[second].sqrtEpsilon == 0.0f)
continue;
RealOpenMM sigma = particles[first].sigmaOver2+particles[second].sigmaOver2;
RealOpenMM epsilon = particles[first].sqrtEpsilon*particles[second].sqrtEpsilon;
energy += computeOneInteraction(first, second, sigma, epsilon, positions, forces, torques, boxVectors);
}
}
}
}
}
// Compute exceptions.
int numExceptions = exceptions.size();
for (int i = 0; i < numExceptions; i++) {
ExceptionInfo& e = exceptions[i];
energy += computeOneInteraction(e.particle1, e.particle2, e.sigma, e.epsilon, positions, forces, torques, boxVectors);
}
// Apply torques.
applyTorques(positions, forces, torques);
return energy;
}
void CpuGayBerneForce::computeEllipsoidFrames(const vector<RealVec>& positions) {
int numParticles = particles.size();
for (int particle = 0; particle < numParticles; particle++) {
ParticleInfo& p = particles[particle];
// Compute the local coordinate system of the ellipsoid;
RealVec xdir, ydir, zdir;
if (p.xparticle == -1) {
xdir = RealVec(1, 0, 0);
ydir = RealVec(0, 1, 0);
}
else {
xdir = positions[particle]-positions[p.xparticle];
xdir /= SQRT(xdir.dot(xdir));
if (p.yparticle == -1) {
if (xdir[1] > -0.5 && xdir[1] < 0.5)
ydir = RealVec(0, 1, 0);
else
ydir = RealVec(1, 0, 0);
}
else
ydir = positions[particle]-positions[p.yparticle];
ydir -= xdir*(xdir.dot(ydir));
ydir /= SQRT(ydir.dot(ydir));
}
zdir = xdir.cross(ydir);
// Compute matrices we will need later.
RealOpenMM (&a)[3][3] = A[particle].v;
RealOpenMM (&b)[3][3] = B[particle].v;
RealOpenMM (&g)[3][3] = G[particle].v;
a[0][0] = xdir[0];
a[0][1] = xdir[1];
a[0][2] = xdir[2];
a[1][0] = ydir[0];
a[1][1] = ydir[1];
a[1][2] = ydir[2];
a[2][0] = zdir[0];
a[2][1] = zdir[1];
a[2][2] = zdir[2];
RealVec r2(p.rx*p.rx, p.ry*p.ry, p.rz*p.rz);
RealVec e2(1/sqrt(p.ex), 1/sqrt(p.ey), 1/sqrt(p.ez));
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) {
b[i][j] = 0;
g[i][j] = 0;
for (int k = 0; k < 3; k++) {
b[i][j] += a[k][i]*e2[k]*a[k][j];
g[i][j] += a[k][i]*r2[k]*a[k][j];
}
}
}
}
void CpuGayBerneForce::applyTorques(const vector<RealVec>& positions, vector<RealVec>& forces, const vector<RealVec>& torques) {
int numParticles = particles.size();
for (int particle = 0; particle < numParticles; particle++) {
ParticleInfo& p = particles[particle];
RealVec pos = positions[particle];
if (p.xparticle != -1) {
// Apply a force to the x particle.
RealVec dx = positions[p.xparticle]-pos;
double dx2 = dx.dot(dx);
RealVec f = torques[particle].cross(dx)/dx2;
forces[p.xparticle] += f;
forces[particle] -= f;
if (p.yparticle != -1) {
// Apply a force to the y particle. This is based on the component of the torque
// that was not already applied to the x particle.
RealVec dy = positions[p.yparticle]-pos;
double dy2 = dy.dot(dy);
RealVec torque = dx*(torques[particle].dot(dx)/dx2);
f = torque.cross(dy)/dy2;
forces[p.yparticle] += f;
forces[particle] -= f;
}
}
}
}
RealOpenMM CpuGayBerneForce::computeOneInteraction(int particle1, int particle2, RealOpenMM sigma, RealOpenMM epsilon, const vector<RealVec>& positions,
vector<RealVec>& forces, vector<RealVec>& torques, const RealVec* boxVectors) {
// Compute the displacement and check against the cutoff.
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (nonbondedMethod == GayBerneForce::CutoffPeriodic)
ReferenceForce::getDeltaRPeriodic(positions[particle2], positions[particle1], boxVectors, deltaR);
else
ReferenceForce::getDeltaR(positions[particle2], positions[particle1], deltaR);
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (nonbondedMethod != GayBerneForce::NoCutoff && r >= cutoffDistance)
return 0;
RealOpenMM rInv = 1/r;
RealVec dr(deltaR[ReferenceForce::XIndex], deltaR[ReferenceForce::YIndex], deltaR[ReferenceForce::ZIndex]);
RealVec drUnit = dr*rInv;
// Compute the switching function.
RealOpenMM switchValue = 1, switchDeriv = 0;
if (useSwitchingFunction && r > switchingDistance) {
RealOpenMM t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
switchValue = 1+t*t*t*(-10+t*(15-t*6));
switchDeriv = t*t*(-30+t*(60-t*30))/(cutoffDistance-switchingDistance);
}
// Interactions between two point particles can be computed more easily.
if (particles[particle1].isPointParticle && particles[particle2].isPointParticle) {
RealOpenMM sig2 = sigma*rInv;
sig2 *= sig2;
RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM energy = 4*epsilon*(sig6-1)*sig6;
RealVec force = drUnit*(switchValue*4*epsilon*(12*sig6 - 6)*sig6*rInv - energy*switchDeriv);
forces[particle1] += force;
forces[particle2] -= force;
return energy*switchValue;
}
// Compute vectors and matrices we'll be needing.
Matrix B12 = B[particle1]+B[particle2];
Matrix G12 = G[particle1]+G[particle2];
Matrix B12inv = B12.inverse();
Matrix G12inv = G12.inverse();
RealOpenMM detG12 = G12.determinant();
// Estimate the distance between the ellipsoids and compute the first terms needed for the energy.
RealOpenMM sigma12 = 1/SQRT(0.5*drUnit.dot(G12inv*drUnit));
RealOpenMM h12 = r - sigma12;
RealOpenMM rho = sigma/(h12+sigma);
RealOpenMM rho2 = rho*rho;
RealOpenMM rho6 = rho2*rho2*rho2;
RealOpenMM u = 4*epsilon*(rho6*rho6-rho6);
RealOpenMM eta = SQRT(2*s[particle1]*s[particle2]/detG12);
RealOpenMM chi = 2*drUnit.dot(B12inv*drUnit);
chi *= chi;
RealOpenMM energy = u*eta*chi;
// Compute the terms needed for the force.
RealVec kappa = G12inv*dr;
RealVec iota = B12inv*dr;
RealOpenMM rInv2 = rInv*rInv;
RealOpenMM dUSLJdr = 24*epsilon*(2*rho6-1)*rho6*rho/sigma;
RealOpenMM temp = 0.5*sigma12*sigma12*sigma12*rInv2;
RealVec dudr = (drUnit + (kappa-drUnit*kappa.dot(drUnit))*temp)*dUSLJdr;
RealVec dchidr = (iota-drUnit*iota.dot(drUnit))*(-8*rInv2*SQRT(chi));
RealVec force = (dchidr*u + dudr*chi)*(eta*switchValue) - drUnit*(energy*switchDeriv);
forces[particle1] += force;
forces[particle2] -= force;
// Compute the terms needed for the torque.
for (int j = 0; j < 2; j++) {
int particle = (j == 0 ? particle1 : particle2);
RealVec dudq = (kappa*G[particle]).cross(kappa*(temp*dUSLJdr));
RealVec dchidq = (iota*B[particle]).cross(iota)*(-4*rInv2);
RealOpenMM (&g12)[3][3] = G12.v;
RealOpenMM (&a)[3][3] = A[particle].v;
ParticleInfo& p = particles[particle];
RealVec scale = RealVec(p.rx*p.rx, p.ry*p.ry, p.rz*p.rz)*(-0.5*eta/detG12);
Matrix D;
RealOpenMM (&d)[3][3] = D.v;
d[0][0] = scale[0]*(g12[1][2]*g12[0][1]*a[0][2] + 2*g12[1][1]*g12[2][2]*a[0][0] -
g12[1][1]*a[0][2]*g12[0][2] - 2*g12[1][2]*a[0][0]*g12[2][1] +
a[0][1]*g12[0][2]*g12[2][1] - a[0][1]*g12[0][1]*g12[2][2] -
g12[1][0]*g12[2][2]*a[0][1] + g12[2][0]*g12[1][2]*a[0][1] +
g12[1][0]*a[0][2]*g12[2][1] - a[0][2]*g12[2][0]*g12[1][1]);
d[0][1] = scale[0]*( g12[0][2]*a[0][0]*g12[2][1] - g12[2][2]*a[0][0]*g12[0][1] +
2*g12[0][0]*g12[2][2]*a[0][1] - g12[0][0]*a[0][2]*g12[1][2] -
2*g12[2][0]*g12[0][2]*a[0][1] + a[0][2]*g12[1][0]*g12[0][2] -
g12[2][2]*g12[1][0]*a[0][0] + g12[2][0]*a[0][0]*g12[1][2] +
g12[2][0]*a[0][2]*g12[0][1] - a[0][2]*g12[0][0]*g12[2][1]);
d[0][2] = scale[0]*( g12[0][1]*g12[1][2]*a[0][0] - g12[0][2]*a[0][0]*g12[1][1] -
g12[0][0]*g12[1][2]*a[0][1] + g12[1][0]*g12[0][2]*a[0][1] -
a[0][1]*g12[0][0]*g12[2][1] - g12[2][0]*g12[1][1]*a[0][0] +
2*g12[1][1]*g12[0][0]*a[0][2] - 2*g12[1][0]*a[0][2]*g12[0][1] +
g12[1][0]*g12[2][1]*a[0][0] + g12[2][0]*a[0][1]*g12[0][1]);
d[1][0] = scale[1]*(-g12[1][1]*a[1][2]*g12[0][2] + 2*g12[1][1]*g12[2][2]*a[1][0] +
g12[1][2]*g12[0][1]*a[1][2] - 2*g12[1][2]*a[1][0]*g12[2][1] +
a[1][1]*g12[0][2]*g12[2][1] - a[1][1]*g12[0][1]*g12[2][2] -
g12[1][0]*g12[2][2]*a[1][1] + g12[2][0]*g12[1][2]*a[1][1] -
a[1][2]*g12[2][0]*g12[1][1] + g12[1][0]*a[1][2]*g12[2][1]);
d[1][1] = scale[1]*( g12[0][2]*a[1][0]*g12[2][1] - g12[0][1]*g12[2][2]*a[1][0] +
2*g12[2][2]*g12[0][0]*a[1][1] - a[1][2]*g12[0][0]*g12[1][2] -
2*g12[2][0]*a[1][1]*g12[0][2] - g12[1][0]*g12[2][2]*a[1][0] +
g12[2][0]*g12[1][2]*a[1][0] + g12[1][0]*a[1][2]*g12[0][2] -
g12[0][0]*a[1][2]*g12[2][1] + a[1][2]*g12[0][1]*g12[2][0]);
d[1][2] = scale[1]*( g12[0][1]*g12[1][2]*a[1][0] - g12[0][2]*a[1][0]*g12[1][1] -
g12[0][0]*g12[1][2]*a[1][1] + g12[1][0]*g12[0][2]*a[1][1] +
2*g12[1][1]*g12[0][0]*a[1][2] - g12[0][0]*a[1][1]*g12[2][1] +
g12[0][1]*g12[2][0]*a[1][1] - a[1][0]*g12[2][0]*g12[1][1] -
2*g12[1][0]*g12[0][1]*a[1][2] + g12[1][0]*a[1][0]*g12[2][1]);
d[2][0] = scale[2]*( -g12[1][1]*g12[0][2]*a[2][2] + g12[0][1]*g12[1][2]*a[2][2] +
2*g12[1][1]*a[2][0]*g12[2][2] - g12[0][1]*a[2][1]*g12[2][2] +
g12[0][2]*g12[2][1]*a[2][1] - 2*a[2][0]*g12[2][1]*g12[1][2] -
g12[1][0]*a[2][1]*g12[2][2] + g12[1][2]*g12[2][0]*a[2][1] -
g12[1][1]*g12[2][0]*a[2][2] + g12[2][1]*g12[1][0]*a[2][2]);
d[2][1] = scale[2]*( -g12[0][1]*g12[2][2]*a[2][0] + g12[0][2]*a[2][0]*g12[2][1] +
2*a[2][1]*g12[0][0]*g12[2][2] - g12[1][2]*a[2][2]*g12[0][0] -
2*a[2][1]*g12[0][2]*g12[2][0] - g12[1][0]*a[2][0]*g12[2][2] +
g12[1][0]*g12[0][2]*a[2][2] + g12[1][2]*g12[2][0]*a[2][0] -
g12[0][0]*a[2][2]*g12[2][1] + a[2][2]*g12[0][1]*g12[2][0]);
d[2][2] = scale[2]*( g12[0][1]*g12[1][2]*a[2][0] - g12[0][2]*a[2][0]*g12[1][1] -
g12[0][0]*g12[1][2]*a[2][1] + g12[1][0]*g12[0][2]*a[2][1] -
g12[1][1]*g12[2][0]*a[2][0] - g12[2][1]*a[2][1]*g12[0][0] +
2*g12[1][1]*a[2][2]*g12[0][0] + g12[2][1]*g12[1][0]*a[2][0] +
g12[2][0]*g12[0][1]*a[2][1] - 2*a[2][2]*g12[1][0]*g12[0][1]);
RealVec detadq;
for (int i = 0; i < 3; i++)
detadq += RealVec(a[i][0], a[i][1], a[i][2]).cross(RealVec(d[i][0], d[i][1], d[i][2]));
RealVec torque = (dchidq*(u*eta) + detadq*(u*chi) + dudq*(eta*chi))*switchValue;
torques[particle] -= torque;
}
return switchValue*energy;
}
platforms/cpu/src/CpuKernelFactory.cpp
View file @ 6c1d5eb1
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013 Stanford University and the Authors. *
* Portions copyright (c) 2013-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -57,6 +57,8 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform&
return new CpuCalcGBSAOBCForceKernel(name, platform, data);
if (name == CalcCustomGBForceKernel::Name())
return new CpuCalcCustomGBForceKernel(name, platform, data);
if (name == CalcGayBerneForceKernel::Name())
return new CpuCalcGayBerneForceKernel(name, platform, data);
if (name == IntegrateLangevinStepKernel::Name())
return new CpuIntegrateLangevinStepKernel(name, platform, data);
throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str());
......
platforms/cpu/src/CpuKernels.cpp
View file @ 6c1d5eb1
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013-2015 Stanford University and the Authors. *
* Portions copyright (c) 2013-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -1205,6 +1205,25 @@ void CpuCalcCustomManyParticleForceKernel::copyParametersToContext(ContextImpl&
}
}
CpuCalcGayBerneForceKernel::~CpuCalcGayBerneForceKernel() {
if (ixn != NULL)
delete ixn;
}
void CpuCalcGayBerneForceKernel::initialize(const System& system, const GayBerneForce& force) {
ixn = new CpuGayBerneForce(force);
}
double CpuCalcGayBerneForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
return ixn->calculateForce(extractPositions(context), extractForces(context), extractBoxVectors(context), data);
}
void CpuCalcGayBerneForceKernel::copyParametersToContext(ContextImpl& context, const GayBerneForce& force) {
delete ixn;
ixn = NULL;
ixn = new CpuGayBerneForce(force);
}
CpuIntegrateLangevinStepKernel::~CpuIntegrateLangevinStepKernel() {
if (dynamics)
delete dynamics;
......
platforms/cpu/src/CpuPlatform.cpp
View file @ 6c1d5eb1
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013-2014 Stanford University and the Authors. *
* Portions copyright (c) 2013-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -69,6 +69,7 @@ CpuPlatform::CpuPlatform() {
registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
platformProperties.push_back(CpuThreads());
int threads = getNumProcessors();
......
platforms/cpu/tests/TestCpuGayBerneForce.cpp 0 → 100644
View file @ 6c1d5eb1
/* -------------------------------------------------------------------------- *
* 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) 2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "CpuTests.h"
#include "TestGayBerneForce.h"
void runPlatformTests() {
}
platforms/reference/include/ReferenceGayBerneForce.h
View file @ 6c1d5eb1
......@@ -95,17 +95,6 @@ struct ReferenceGayBerneForce::Matrix {
v[2][0]*r[0] + v[2][1]*r[1] + v[2][2]*r[2]);
}
Matrix operator*(const Matrix& m) {
Matrix result;
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) {
result.v[i][j] = 0;
for (int k = 0; k < 3; k++)
result.v[i][j] += v[i][k]*m.v[k][j];
}
return result;
}
Matrix operator+(const Matrix& m) {
Matrix result;
for (int i = 0; i < 3; i++)
......
tests/TestGayBerneForce.h
View file @ 6c1d5eb1
......@@ -184,6 +184,42 @@ void testEnergyConservation() {
}
}
void testExceptions() {
// Create two Lennard-Jones particles for which the energy scale factors vary,
// then override their interaction with an exception.
const double sigma = 0.5;
const double epsilon = 1.5;
System system;
for (int i = 0; i < 6; i++)
system.addParticle(1.0);
GayBerneForce* gb = new GayBerneForce();
system.addForce(gb);
gb->addParticle(sigma, epsilon, 1, 2, sigma, sigma, sigma, 1.1, 1.5, 1.8);
gb->addParticle(1, 0, -1, -1, 1, 1, 1, 1, 1, 1);
gb->addParticle(1, 0, -1, -1, 1, 1, 1, 1, 1, 1);
gb->addParticle(sigma, epsilon, 4, 5, sigma, sigma, sigma, 1.2, 1.6, 1.7);
gb->addParticle(1, 0, -1, -1, 1, 1, 1, 1, 1, 1);
gb->addParticle(1, 0, -1, -1, 1, 1, 1, 1, 1, 1);
gb->addException(0, 3, sigma, 3.5*epsilon);
vector<Vec3> positions(6);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
positions[3] = Vec3(1, 0, 0);
positions[4] = Vec3(2, 0, 0);
positions[5] = Vec3(1, 1, 0);
VerletIntegrator integ(0.001);
Context context(system, integ, platform);
context.setPositions(positions);
double expectedEnergy = 3.5*4*epsilon*(pow(sigma, 12.0)-pow(sigma, 6.0));
double expectedForce = 3.5*4*epsilon*(12*pow(sigma, 12.0)-6*pow(sigma, 6.0));
double expectedScale = pow(2.0/(1/sqrt(1.1) + 1/sqrt(1.2)), 2.0);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(expectedEnergy*expectedScale, state.getPotentialEnergy(), 1e-5);
ASSERT_EQUAL_VEC(Vec3(expectedForce*expectedScale, 0, 0), state.getForces()[3], 1e-5);
}
void runPlatformTests();
int main(int argc, char* argv[]) {
......@@ -192,6 +228,7 @@ int main(int argc, char* argv[]) {
testPointParticles();
testEnergyScales();
testEnergyConservation();
testExceptions();
runPlatformTests();
}
catch(const exception& e) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment