Skip to content

GitLab

"plugins/vscode:/vscode.git/clone" did not exist on "8fc504ed5d5c764a1f3c8989133f4fda84bc2a76"
Commit cf8a03e8 authored by peastman's avatar peastman
Browse files

Merged changes from main branch

parents f7f70136 31d02cdc
Hide whitespace changes
Inline Side-by-side
Showing with 1592 additions and 65 deletions
openmmapi/src/CustomNonbondedForceImpl.cpp
View file @ cf8a03e8
......@@ -203,11 +203,11 @@ double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedFo
// Count the total number of particle pairs for each pair of classes.
map<pair<int, int>, int> interactionCount;
map<pair<int, int>, long long int> interactionCount;
if (force.getNumInteractionGroups() == 0) {
// Count the particles of each class.
vector<int> classCounts(numClasses, 0);
vector<long long int> classCounts(numClasses, 0);
for (int i = 0; i < numParticles; i++)
classCounts[atomClass[i]]++;
for (int i = 0; i < numClasses; i++) {
......@@ -246,9 +246,10 @@ double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedFo
for (int i = 0; i < numClasses; i++)
for (int j = i; j < numClasses; j++)
sum += interactionCount[make_pair(i, j)]*integrateInteraction(expression, classes[i], classes[j], force, context);
int numInteractions = (numParticles*(numParticles+1))/2;
double nPart = (double) numParticles;
double numInteractions = (nPart*(nPart+1))/2;
sum /= numInteractions;
return 2*M_PI*numParticles*numParticles*sum;
return 2*M_PI*nPart*nPart*sum;
}
double CustomNonbondedForceImpl::integrateInteraction(Lepton::CompiledExpression& expression, const vector<double>& params1, const vector<double>& params2,
......
openmmapi/src/NonbondedForce.cpp
View file @ cf8a03e8
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2013 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -48,7 +48,7 @@ using std::stringstream;
using std::vector;
NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), switchingDistance(-1.0), rfDielectric(78.3),
ewaldErrorTol(5e-4), useSwitchingFunction(false), useDispersionCorrection(true), recipForceGroup(-1) {
ewaldErrorTol(5e-4), alpha(0.0), useSwitchingFunction(false), useDispersionCorrection(true), recipForceGroup(-1), nx(0), ny(0), nz(0) {
}
NonbondedForce::NonbondedMethod NonbondedForce::getNonbondedMethod() const {
......@@ -95,11 +95,24 @@ double NonbondedForce::getEwaldErrorTolerance() const {
return ewaldErrorTol;
}
void NonbondedForce::setEwaldErrorTolerance(double tol)
{
void NonbondedForce::setEwaldErrorTolerance(double tol) {
ewaldErrorTol = tol;
}
void NonbondedForce::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
alpha = this->alpha;
nx = this->nx;
ny = this->ny;
nz = this->nz;
}
void NonbondedForce::setPMEParameters(double alpha, int nx, int ny, int nz) {
this->alpha = alpha;
this->nx = nx;
this->ny = ny;
this->nz = nz;
}
int NonbondedForce::addParticle(double charge, double sigma, double epsilon) {
particles.push_back(ParticleInfo(charge, sigma, epsilon));
return particles.size()-1;
......
openmmapi/src/NonbondedForceImpl.cpp
View file @ cf8a03e8
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2010 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -149,16 +149,19 @@ void NonbondedForceImpl::calcEwaldParameters(const System& system, const Nonbond
}
void NonbondedForceImpl::calcPMEParameters(const System& system, const NonbondedForce& force, double& alpha, int& xsize, int& ysize, int& zsize) {
Vec3 boxVectors[3];
system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
double tol = force.getEwaldErrorTolerance();
alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(2.0*tol));
xsize = (int) ceil(2*alpha*boxVectors[0][0]/(3*pow(tol, 0.2)));
ysize = (int) ceil(2*alpha*boxVectors[1][1]/(3*pow(tol, 0.2)));
zsize = (int) ceil(2*alpha*boxVectors[2][2]/(3*pow(tol, 0.2)));
xsize = max(xsize, 5);
ysize = max(ysize, 5);
zsize = max(zsize, 5);
force.getPMEParameters(alpha, xsize, ysize, zsize);
if (alpha == 0.0) {
Vec3 boxVectors[3];
system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
double tol = force.getEwaldErrorTolerance();
alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(2.0*tol));
xsize = (int) ceil(2*alpha*boxVectors[0][0]/(3*pow(tol, 0.2)));
ysize = (int) ceil(2*alpha*boxVectors[1][1]/(3*pow(tol, 0.2)));
zsize = (int) ceil(2*alpha*boxVectors[2][2]/(3*pow(tol, 0.2)));
xsize = max(xsize, 5);
ysize = max(ysize, 5);
zsize = max(zsize, 5);
}
}
int NonbondedForceImpl::findZero(const NonbondedForceImpl::ErrorFunction& f, int initialGuess) {
......@@ -239,7 +242,8 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const
for (map<pair<double, double>, int>::const_iterator entry = classCounts.begin(); entry != classCounts.end(); ++entry) {
double sigma = entry->first.first;
double epsilon = entry->first.second;
int count = (entry->second*(entry->second+1))/2;
double count = (double) entry->second;
count *= (count + 1) / 2;
double sigma2 = sigma*sigma;
double sigma6 = sigma2*sigma2*sigma2;
sum1 += count*epsilon*sigma6*sigma6;
......@@ -251,7 +255,8 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const
for (map<pair<double, double>, int>::const_iterator class2 = classCounts.begin(); class2 != class1; ++class2) {
double sigma = 0.5*(class1->first.first+class2->first.first);
double epsilon = sqrt(class1->first.second*class2->first.second);
int count = class1->second*class2->second;
double count = (double) class1->second;
count *= (double) class2->second;
double sigma2 = sigma*sigma;
double sigma6 = sigma2*sigma2*sigma2;
sum1 += count*epsilon*sigma6*sigma6;
......@@ -259,8 +264,8 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const
if (useSwitch)
sum3 += count*epsilon*(evalIntegral(cutoff, switchDist, cutoff, sigma)-evalIntegral(switchDist, switchDist, cutoff, sigma));
}
int numParticles = system.getNumParticles();
int numInteractions = (numParticles*(numParticles+1))/2;
double numParticles = (double) system.getNumParticles();
double numInteractions = (numParticles*(numParticles+1))/2;
sum1 /= numInteractions;
sum2 /= numInteractions;
sum3 /= numInteractions;
......
openmmapi/src/SplineFitter.cpp
View file @ cf8a03e8
......@@ -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) 2010 Stanford University and the Authors. *
* Portions copyright (c) 2010-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -190,3 +190,535 @@ void SplineFitter::solveTridiagonalMatrix(const vector<double>& a, const vector<
for (int i = n-2; i >= 0; i--)
sol[i] -= gamma[i+1]*sol[i+1];
}
void SplineFitter::create2DNaturalSpline(const vector<double>& x, const vector<double>& y, const vector<double>& values, vector<vector<double> >& c) {
int xsize = x.size(), ysize = y.size();
if (xsize < 2 || ysize < 2)
throw OpenMMException("create2DNaturalSpline: must have at least two points along each axis");
if (values.size() != xsize*ysize)
throw OpenMMException("create2DNaturalSpline: incorrect number of values");
vector<double> d1(xsize*ysize), d2(xsize*ysize), d12(xsize*ysize);
vector<double> t(xsize), deriv(xsize);
// Compute derivatives with respect to x.
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < xsize; j++)
t[j] = values[j+xsize*i];
SplineFitter::createNaturalSpline(x, t, deriv);
for (int j = 0; j < xsize; j++)
d1[j+xsize*i] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
}
// Compute derivatives with respect to y.
t.resize(ysize);
deriv.resize(ysize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++)
t[j] = values[i+xsize*j];
SplineFitter::createNaturalSpline(y, t, deriv);
for (int j = 0; j < ysize; j++)
d2[i+xsize*j] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[j]);
}
// Compute cross derivatives.
t.resize(xsize);
deriv.resize(xsize);
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < xsize; j++)
t[j] = d2[j+xsize*i];
SplineFitter::createNaturalSpline(x, t, deriv);
for (int j = 0; j < xsize; j++)
d12[j+xsize*i] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
}
// Now compute the coefficients.
const int wt[] = {
1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4,
0, 0, 0, 0, 0, 0, 0, 0, 3, 0, -9, 6, -2, 0, 6, -4,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -6, 0, 0, -6, 4,
0, 0, 3, -2, 0, 0, 0, 0, 0, 0, -9, 6, 0, 0, 6, -4,
0, 0, 0, 0, 1, 0, -3, 2, -2, 0, 6, -4, 1, 0, -3, 2,
0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 3, -2, 1, 0, -3, 2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2, 0, 0, 3, -2,
0, 0, 0, 0, 0, 0, 3, -2, 0, 0, -6, 4, 0, 0, 3, -2,
0, 1, -2, 1, 0, 0, 0, 0, 0, -3, 6, -3, 0, 2, -4, 2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -6, 3, 0, -2, 4, -2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 2, -2,
0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 3, -3, 0, 0, -2, 2,
0, 0, 0, 0, 0, 1, -2, 1, 0, -2, 4, -2, 0, 1, -2, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 1, -2, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 1,
0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 2, -2, 0, 0, -1, 1
};
vector<double> rhs(16);
c.resize((xsize-1)*(ysize-1));
for (int i = 0; i < xsize-1; i++) {
for (int j = 0; j < ysize-1; j++) {
// Compute the 16 coefficients for patch (i, j).
int nexti = i+1;
int nextj = j+1;
double deltax = x[nexti]-x[i];
double deltay = y[nextj]-y[j];
double e[] = {values[i+j*xsize], values[nexti+j*xsize], values[nexti+nextj*xsize], values[i+nextj*xsize]};
double e1[] = {d1[i+j*xsize], d1[nexti+j*xsize], d1[nexti+nextj*xsize], d1[i+nextj*xsize]};
double e2[] = {d2[i+j*xsize], d2[nexti+j*xsize], d2[nexti+nextj*xsize], d2[i+nextj*xsize]};
double e12[] = {d12[i+j*xsize], d12[nexti+j*xsize], d12[nexti+nextj*xsize], d12[i+nextj*xsize]};
for (int k = 0; k < 4; k++) {
rhs[k] = e[k];
rhs[k+4] = e1[k]*deltax;
rhs[k+8] = e2[k]*deltay;
rhs[k+12] = e12[k]*deltax*deltay;
}
vector<double>& coeff = c[i+j*(xsize-1)];
coeff.resize(16);
for (int k = 0; k < 16; k++) {
double sum = 0.0;
for (int m = 0; m < 16; m++)
sum += wt[k+16*m]*rhs[m];
coeff[k] = sum;
}
}
}
}
double SplineFitter::evaluate2DSpline(const vector<double>& x, const vector<double>& y, const vector<double>& values, const vector<vector<double> >& c, double u, double v) {
int xsize = x.size();
int ysize = y.size();
if (u < x[0] || u > x[xsize-1] || v < y[0] || v > y[ysize-1])
throw OpenMMException("evaluate2DSpline: specified point is outside the range defined by the spline");
// Perform a binary search to identify the interval containing the point to evaluate.
int lowerx = 0;
int upperx = xsize-1;
while (upperx-lowerx > 1) {
int middle = (upperx+lowerx)/2;
if (x[middle] > u)
upperx = middle;
else
lowerx = middle;
}
int lowery = 0;
int uppery = ysize-1;
while (uppery-lowery > 1) {
int middle = (uppery+lowery)/2;
if (y[middle] > v)
uppery = middle;
else
lowery = middle;
}
double deltax = x[upperx]-x[lowerx];
double deltay = y[uppery]-y[lowery];
double da = (u-x[lowerx])/deltax;
double db = (v-y[lowery])/deltay;
const vector<double>& coeff = c[lowerx+(xsize-1)*lowery];
// Evaluate the spline to determine the value.
double value = 0;
for (int i = 3; i >= 0; i--)
value = da*value + ((coeff[i*4+3]*db + coeff[i*4+2])*db + coeff[i*4+1])*db + coeff[i*4+0];
return value;
}
void SplineFitter::evaluate2DSplineDerivatives(const vector<double>& x, const vector<double>& y, const vector<double>& values, const vector<vector<double> >& c, double u, double v, double& dx, double &dy) {
int xsize = x.size();
int ysize = y.size();
if (u < x[0] || u > x[xsize-1] || v < y[0] || v > y[ysize-1])
throw OpenMMException("evaluate2DSplineDerivatives: specified point is outside the range defined by the spline");
// Perform a binary search to identify the interval containing the point to evaluate.
int lowerx = 0;
int upperx = xsize-1;
while (upperx-lowerx > 1) {
int middle = (upperx+lowerx)/2;
if (x[middle] > u)
upperx = middle;
else
lowerx = middle;
}
int lowery = 0;
int uppery = ysize-1;
while (uppery-lowery > 1) {
int middle = (uppery+lowery)/2;
if (y[middle] > v)
uppery = middle;
else
lowery = middle;
}
double deltax = x[upperx]-x[lowerx];
double deltay = y[uppery]-y[lowery];
double da = (u-x[lowerx])/deltax;
double db = (v-y[lowery])/deltay;
const vector<double>& coeff = c[lowerx+(xsize-1)*lowery];
// Evaluate the spline to determine the derivatives.
dx = 0;
dy = 0;
for (int i = 3; i >= 0; i--) {
dx = db*dx + (3.0*coeff[i+3*4]*da + 2.0*coeff[i+2*4])*da + coeff[i+1*4];
dy = da*dy + (3.0*coeff[i*4+3]*db + 2.0*coeff[i*4+2])*db + coeff[i*4+1];
}
dx /= deltax;
dy /= deltay;
}
void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<double>& y, const vector<double>& z, const vector<double>& values, vector<vector<double> >& c) {
int xsize = x.size(), ysize = y.size(), zsize = z.size();
int xysize = xsize*ysize;
if (xsize < 2 || ysize < 2 || zsize < 2)
throw OpenMMException("create2DNaturalSpline: must have at least two points along each axis");
if (values.size() != xsize*ysize*zsize)
throw OpenMMException("create2DNaturalSpline: incorrect number of values");
vector<double> d1(xsize*ysize*zsize), d2(xsize*ysize*zsize), d3(xsize*ysize*zsize);
vector<double> d12(xsize*ysize*zsize), d13(xsize*ysize*zsize), d23(xsize*ysize*zsize), d123(xsize*ysize*zsize);
vector<double> t(xsize), deriv(xsize);
// Compute derivatives with respect to x.
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < zsize; j++) {
for (int k = 0; k < xsize; k++)
t[k] = values[k+xsize*i+xysize*j];
SplineFitter::createNaturalSpline(x, t, deriv);
for (int k = 0; k < xsize; k++)
d1[k+xsize*i+xysize*j] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[k]);
}
}
// Compute derivatives with respect to y.
t.resize(ysize);
deriv.resize(ysize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < zsize; j++) {
for (int k = 0; k < ysize; k++)
t[k] = values[i+xsize*k+xysize*j];
SplineFitter::createNaturalSpline(y, t, deriv);
for (int k = 0; k < ysize; k++)
d2[i+xsize*k+xysize*j] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[k]);
}
}
// Compute derivatives with respect to z.
t.resize(zsize);
deriv.resize(zsize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
for (int k = 0; k < zsize; k++)
t[k] = values[i+xsize*j+xysize*k];
SplineFitter::createNaturalSpline(z, t, deriv);
for (int k = 0; k < zsize; k++)
d3[i+xsize*j+xysize*k] = SplineFitter::evaluateSplineDerivative(z, t, deriv, z[k]);
}
}
// Compute second derivatives with respect to x and y.
t.resize(xsize);
deriv.resize(xsize);
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < zsize; j++) {
for (int k = 0; k < xsize; k++)
t[k] = d2[k+xsize*i+xysize*j];
SplineFitter::createNaturalSpline(x, t, deriv);
for (int k = 0; k < xsize; k++)
d12[k+xsize*i+xysize*j] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[k]);
}
}
// Compute second derivatives with respect to y and z.
t.resize(ysize);
deriv.resize(ysize);
for (int i = 0; i < zsize; i++) {
for (int j = 0; j < xsize; j++) {
for (int k = 0; k < ysize; k++)
t[k] = d3[j+xsize*k+xysize*i];
SplineFitter::createNaturalSpline(y, t, deriv);
for (int k = 0; k < ysize; k++)
d23[j+xsize*k+xysize*i] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[k]);
}
}
// Compute second derivatives with respect to x and z.
t.resize(zsize);
deriv.resize(zsize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
for (int k = 0; k < zsize; k++)
t[k] = d1[i+xsize*j+xysize*k];
SplineFitter::createNaturalSpline(z, t, deriv);
for (int k = 0; k < zsize; k++)
d13[i+xsize*j+xysize*k] = SplineFitter::evaluateSplineDerivative(z, t, deriv, z[k]);
}
}
// Compute third derivatives with respect to x, y, and z.
t.resize(xsize);
deriv.resize(xsize);
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < zsize; j++) {
for (int k = 0; k < xsize; k++)
t[k] = d23[k+xsize*i+xysize*j];
SplineFitter::createNaturalSpline(x, t, deriv);
for (int k = 0; k < xsize; k++)
d123[k+xsize*i+xysize*j] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[k]);
}
}
// Now compute the coefficients. This involves multiplying by a sparse 64x64 matrix, given
// here in packed form.
const int wt[] = {
1,0,1,
1,8,1,
4,0,-3,1,3,8,-2,9,-1,
4,0,2,1,-2,8,1,9,1,
1,16,1,
1,32,1,
4,16,-3,17,3,32,-2,33,-1,
4,16,2,17,-2,32,1,33,1,
4,0,-3,2,3,16,-2,18,-1,
4,8,-3,10,3,32,-2,34,-1,
16,0,9,1,-9,2,-9,3,9,8,6,9,3,10,-6,11,-3,16,6,17,-6,18,3,19,-3,32,4,33,2,34,2,35,1,
16,0,-6,1,6,2,6,3,-6,8,-3,9,-3,10,3,11,3,16,-4,17,4,18,-2,19,2,32,-2,33,-2,34,-1,35,-1,
4,0,2,2,-2,16,1,18,1,
4,8,2,10,-2,32,1,34,1,
16,0,-6,1,6,2,6,3,-6,8,-4,9,-2,10,4,11,2,16,-3,17,3,18,-3,19,3,32,-2,33,-1,34,-2,35,-1,
16,0,4,1,-4,2,-4,3,4,8,2,9,2,10,-2,11,-2,16,2,17,-2,18,2,19,-2,32,1,33,1,34,1,35,1,
1,24,1,
1,40,1,
4,24,-3,25,3,40,-2,41,-1,
4,24,2,25,-2,40,1,41,1,
1,48,1,
1,56,1,
4,48,-3,49,3,56,-2,57,-1,
4,48,2,49,-2,56,1,57,1,
4,24,-3,26,3,48,-2,50,-1,
4,40,-3,42,3,56,-2,58,-1,
16,24,9,25,-9,26,-9,27,9,40,6,41,3,42,-6,43,-3,48,6,49,-6,50,3,51,-3,56,4,57,2,58,2,59,1,
16,24,-6,25,6,26,6,27,-6,40,-3,41,-3,42,3,43,3,48,-4,49,4,50,-2,51,2,56,-2,57,-2,58,-1,59,-1,
4,24,2,26,-2,48,1,50,1,
4,40,2,42,-2,56,1,58,1,
16,24,-6,25,6,26,6,27,-6,40,-4,41,-2,42,4,43,2,48,-3,49,3,50,-3,51,3,56,-2,57,-1,58,-2,59,-1,
16,24,4,25,-4,26,-4,27,4,40,2,41,2,42,-2,43,-2,48,2,49,-2,50,2,51,-2,56,1,57,1,58,1,59,1,
4,0,-3,4,3,24,-2,28,-1,
4,8,-3,12,3,40,-2,44,-1,
16,0,9,1,-9,4,-9,5,9,8,6,9,3,12,-6,13,-3,24,6,25,-6,28,3,29,-3,40,4,41,2,44,2,45,1,
16,0,-6,1,6,4,6,5,-6,8,-3,9,-3,12,3,13,3,24,-4,25,4,28,-2,29,2,40,-2,41,-2,44,-1,45,-1,
4,16,-3,20,3,48,-2,52,-1,
4,32,-3,36,3,56,-2,60,-1,
16,16,9,17,-9,20,-9,21,9,32,6,33,3,36,-6,37,-3,48,6,49,-6,52,3,53,-3,56,4,57,2,60,2,61,1,
16,16,-6,17,6,20,6,21,-6,32,-3,33,-3,36,3,37,3,48,-4,49,4,52,-2,53,2,56,-2,57,-2,60,-1,61,-1,
16,0,9,2,-9,4,-9,6,9,16,6,18,3,20,-6,22,-3,24,6,26,-6,28,3,30,-3,48,4,50,2,52,2,54,1,
16,8,9,10,-9,12,-9,14,9,32,6,34,3,36,-6,38,-3,40,6,42,-6,44,3,46,-3,56,4,58,2,60,2,62,1,
64,0,-27,1,27,2,27,3,-27,4,27,5,-27,6,-27,7,27,8,-18,9,-9,10,18,11,9,12,18,13,9,14,-18,15,-9,16,-18,17,18,18,-9,19,9,20,18,21,-18,22,9,23,-9,24,-18,25,18,26,18,27,-18,28,-9,29,9,30,9,31,-9,32,-12,33,-6,34,-6,35,-3,36,12,37,6,38,6,39,3,40,-12,41,-6,42,12,43,6,44,-6,45,-3,46,6,47,3,48,-12,49,12,50,-6,51,6,52,-6,53,6,54,-3,55,3,56,-8,57,-4,58,-4,59,-2,60,-4,61,-2,62,-2,63,-1,
64,0,18,1,-18,2,-18,3,18,4,-18,5,18,6,18,7,-18,8,9,9,9,10,-9,11,-9,12,-9,13,-9,14,9,15,9,16,12,17,-12,18,6,19,-6,20,-12,21,12,22,-6,23,6,24,12,25,-12,26,-12,27,12,28,6,29,-6,30,-6,31,6,32,6,33,6,34,3,35,3,36,-6,37,-6,38,-3,39,-3,40,6,41,6,42,-6,43,-6,44,3,45,3,46,-3,47,-3,48,8,49,-8,50,4,51,-4,52,4,53,-4,54,2,55,-2,56,4,57,4,58,2,59,2,60,2,61,2,62,1,63,1,
16,0,-6,2,6,4,6,6,-6,16,-3,18,-3,20,3,22,3,24,-4,26,4,28,-2,30,2,48,-2,50,-2,52,-1,54,-1,
16,8,-6,10,6,12,6,14,-6,32,-3,34,-3,36,3,38,3,40,-4,42,4,44,-2,46,2,56,-2,58,-2,60,-1,62,-1,
64,0,18,1,-18,2,-18,3,18,4,-18,5,18,6,18,7,-18,8,12,9,6,10,-12,11,-6,12,-12,13,-6,14,12,15,6,16,9,17,-9,18,9,19,-9,20,-9,21,9,22,-9,23,9,24,12,25,-12,26,-12,27,12,28,6,29,-6,30,-6,31,6,32,6,33,3,34,6,35,3,36,-6,37,-3,38,-6,39,-3,40,8,41,4,42,-8,43,-4,44,4,45,2,46,-4,47,-2,48,6,49,-6,50,6,51,-6,52,3,53,-3,54,3,55,-3,56,4,57,2,58,4,59,2,60,2,61,1,62,2,63,1,
64,0,-12,1,12,2,12,3,-12,4,12,5,-12,6,-12,7,12,8,-6,9,-6,10,6,11,6,12,6,13,6,14,-6,15,-6,16,-6,17,6,18,-6,19,6,20,6,21,-6,22,6,23,-6,24,-8,25,8,26,8,27,-8,28,-4,29,4,30,4,31,-4,32,-3,33,-3,34,-3,35,-3,36,3,37,3,38,3,39,3,40,-4,41,-4,42,4,43,4,44,-2,45,-2,46,2,47,2,48,-4,49,4,50,-4,51,4,52,-2,53,2,54,-2,55,2,56,-2,57,-2,58,-2,59,-2,60,-1,61,-1,62,-1,63,-1,
4,0,2,4,-2,24,1,28,1,
4,8,2,12,-2,40,1,44,1,
16,0,-6,1,6,4,6,5,-6,8,-4,9,-2,12,4,13,2,24,-3,25,3,28,-3,29,3,40,-2,41,-1,44,-2,45,-1,
16,0,4,1,-4,4,-4,5,4,8,2,9,2,12,-2,13,-2,24,2,25,-2,28,2,29,-2,40,1,41,1,44,1,45,1,
4,16,2,20,-2,48,1,52,1,
4,32,2,36,-2,56,1,60,1,
16,16,-6,17,6,20,6,21,-6,32,-4,33,-2,36,4,37,2,48,-3,49,3,52,-3,53,3,56,-2,57,-1,60,-2,61,-1,
16,16,4,17,-4,20,-4,21,4,32,2,33,2,36,-2,37,-2,48,2,49,-2,52,2,53,-2,56,1,57,1,60,1,61,1,
16,0,-6,2,6,4,6,6,-6,16,-4,18,-2,20,4,22,2,24,-3,26,3,28,-3,30,3,48,-2,50,-1,52,-2,54,-1,
16,8,-6,10,6,12,6,14,-6,32,-4,34,-2,36,4,38,2,40,-3,42,3,44,-3,46,3,56,-2,58,-1,60,-2,62,-1,
64,0,18,1,-18,2,-18,3,18,4,-18,5,18,6,18,7,-18,8,12,9,6,10,-12,11,-6,12,-12,13,-6,14,12,15,6,16,12,17,-12,18,6,19,-6,20,-12,21,12,22,-6,23,6,24,9,25,-9,26,-9,27,9,28,9,29,-9,30,-9,31,9,32,8,33,4,34,4,35,2,36,-8,37,-4,38,-4,39,-2,40,6,41,3,42,-6,43,-3,44,6,45,3,46,-6,47,-3,48,6,49,-6,50,3,51,-3,52,6,53,-6,54,3,55,-3,56,4,57,2,58,2,59,1,60,4,61,2,62,2,63,1,
64,0,-12,1,12,2,12,3,-12,4,12,5,-12,6,-12,7,12,8,-6,9,-6,10,6,11,6,12,6,13,6,14,-6,15,-6,16,-8,17,8,18,-4,19,4,20,8,21,-8,22,4,23,-4,24,-6,25,6,26,6,27,-6,28,-6,29,6,30,6,31,-6,32,-4,33,-4,34,-2,35,-2,36,4,37,4,38,2,39,2,40,-3,41,-3,42,3,43,3,44,-3,45,-3,46,3,47,3,48,-4,49,4,50,-2,51,2,52,-4,53,4,54,-2,55,2,56,-2,57,-2,58,-1,59,-1,60,-2,61,-2,62,-1,63,-1,
16,0,4,2,-4,4,-4,6,4,16,2,18,2,20,-2,22,-2,24,2,26,-2,28,2,30,-2,48,1,50,1,52,1,54,1,
16,8,4,10,-4,12,-4,14,4,32,2,34,2,36,-2,38,-2,40,2,42,-2,44,2,46,-2,56,1,58,1,60,1,62,1,
64,0,-12,1,12,2,12,3,-12,4,12,5,-12,6,-12,7,12,8,-8,9,-4,10,8,11,4,12,8,13,4,14,-8,15,-4,16,-6,17,6,18,-6,19,6,20,6,21,-6,22,6,23,-6,24,-6,25,6,26,6,27,-6,28,-6,29,6,30,6,31,-6,32,-4,33,-2,34,-4,35,-2,36,4,37,2,38,4,39,2,40,-4,41,-2,42,4,43,2,44,-4,45,-2,46,4,47,2,48,-3,49,3,50,-3,51,3,52,-3,53,3,54,-3,55,3,56,-2,57,-1,58,-2,59,-1,60,-2,61,-1,62,-2,63,-1,
64,0,8,1,-8,2,-8,3,8,4,-8,5,8,6,8,7,-8,8,4,9,4,10,-4,11,-4,12,-4,13,-4,14,4,15,4,16,4,17,-4,18,4,19,-4,20,-4,21,4,22,-4,23,4,24,4,25,-4,26,-4,27,4,28,4,29,-4,30,-4,31,4,32,2,33,2,34,2,35,2,36,-2,37,-2,38,-2,39,-2,40,2,41,2,42,-2,43,-2,44,2,45,2,46,-2,47,-2,48,2,49,-2,50,2,51,-2,52,2,53,-2,54,2,55,-2,56,1,57,1,58,1,59,1,60,1,61,1,62,1,63,1
};
vector<vector<int> > weight(64);
int index = 0;
for (int i = 0; i < 64; i++) {
int numElements = wt[index++];
for (int j = 0; j < numElements; j++) {
weight[i].push_back(wt[index++]);
weight[i].push_back(wt[index++]);
}
}
vector<double> rhs(64);
c.resize((xsize-1)*(ysize-1)*(zsize-1));
for (int i = 0; i < xsize-1; i++) {
for (int j = 0; j < ysize-1; j++) {
for (int k = 0; k < zsize-1; k++) {
// Compute the 64 coefficients for patch (i, j, k).
int nexti = i+1;
int nextj = j+1;
int nextk = k+1;
double deltax = x[nexti]-x[i];
double deltay = y[nextj]-y[j];
double deltaz = z[nextj]-z[j];
double e[] = {values[i+j*xsize+k*xysize], values[nexti+j*xsize+k*xysize], values[i+nextj*xsize+k*xysize], values[nexti+nextj*xsize+k*xysize], values[i+j*xsize+nextk*xysize], values[nexti+j*xsize+nextk*xysize], values[i+nextj*xsize+nextk*xysize], values[nexti+nextj*xsize+nextk*xysize]};
double e1[] = {d1[i+j*xsize+k*xysize], d1[nexti+j*xsize+k*xysize], d1[i+nextj*xsize+k*xysize], d1[nexti+nextj*xsize+k*xysize], d1[i+j*xsize+nextk*xysize], d1[nexti+j*xsize+nextk*xysize], d1[i+nextj*xsize+nextk*xysize], d1[nexti+nextj*xsize+nextk*xysize]};
double e2[] = {d2[i+j*xsize+k*xysize], d2[nexti+j*xsize+k*xysize], d2[i+nextj*xsize+k*xysize], d2[nexti+nextj*xsize+k*xysize], d2[i+j*xsize+nextk*xysize], d2[nexti+j*xsize+nextk*xysize], d2[i+nextj*xsize+nextk*xysize], d2[nexti+nextj*xsize+nextk*xysize]};
double e3[] = {d3[i+j*xsize+k*xysize], d3[nexti+j*xsize+k*xysize], d3[i+nextj*xsize+k*xysize], d3[nexti+nextj*xsize+k*xysize], d3[i+j*xsize+nextk*xysize], d3[nexti+j*xsize+nextk*xysize], d3[i+nextj*xsize+nextk*xysize], d3[nexti+nextj*xsize+nextk*xysize]};
double e12[] = {d12[i+j*xsize+k*xysize], d12[nexti+j*xsize+k*xysize], d12[i+nextj*xsize+k*xysize], d12[nexti+nextj*xsize+k*xysize], d12[i+j*xsize+nextk*xysize], d12[nexti+j*xsize+nextk*xysize], d12[i+nextj*xsize+nextk*xysize], d12[nexti+nextj*xsize+nextk*xysize]};
double e13[] = {d13[i+j*xsize+k*xysize], d13[nexti+j*xsize+k*xysize], d13[i+nextj*xsize+k*xysize], d13[nexti+nextj*xsize+k*xysize], d13[i+j*xsize+nextk*xysize], d13[nexti+j*xsize+nextk*xysize], d13[i+nextj*xsize+nextk*xysize], d13[nexti+nextj*xsize+nextk*xysize]};
double e23[] = {d23[i+j*xsize+k*xysize], d23[nexti+j*xsize+k*xysize], d23[i+nextj*xsize+k*xysize], d23[nexti+nextj*xsize+k*xysize], d23[i+j*xsize+nextk*xysize], d23[nexti+j*xsize+nextk*xysize], d23[i+nextj*xsize+nextk*xysize], d23[nexti+nextj*xsize+nextk*xysize]};
double e123[] = {d123[i+j*xsize+k*xysize], d123[nexti+j*xsize+k*xysize], d123[i+nextj*xsize+k*xysize], d123[nexti+nextj*xsize+k*xysize], d123[i+j*xsize+nextk*xysize], d123[nexti+j*xsize+nextk*xysize], d123[i+nextj*xsize+nextk*xysize], d123[nexti+nextj*xsize+nextk*xysize]};
for (int m = 0; m < 8; m++) {
rhs[m] = e[m];
rhs[m+8] = e1[m]*deltax;
rhs[m+16] = e2[m]*deltay;
rhs[m+24] = e3[m]*deltaz;
rhs[m+32] = e12[m]*deltax*deltay;
rhs[m+40] = e13[m]*deltax*deltaz;
rhs[m+48] = e23[m]*deltay*deltaz;
rhs[m+56] = e123[m]*deltax*deltay*deltaz;
}
vector<double>& coeff = c[i+j*(xsize-1)+k*(xsize-1)*(ysize-1)];
coeff.resize(64);
for (int m = 0; m < 64; m++) {
double sum = 0.0;
int numElements = weight[m].size();
for (int n = 0; n < numElements; n += 2)
sum += weight[m][n+1]*rhs[weight[m][n]];
coeff[m] = sum;
}
}
}
}
}
double SplineFitter::evaluate3DSpline(const vector<double>& x, const vector<double>& y, const vector<double>& z, const vector<double>& values, const vector<vector<double> >& c, double u, double v, double w) {
int xsize = x.size();
int ysize = y.size();
int zsize = z.size();
if (u < x[0] || u > x[xsize-1] || v < y[0] || v > y[ysize-1] || w < z[0] || w > z[zsize-1])
throw OpenMMException("evaluate3DSpline: specified point is outside the range defined by the spline");
// Perform a binary search to identify the interval containing the point to evaluate.
int lowerx = 0;
int upperx = xsize-1;
while (upperx-lowerx > 1) {
int middle = (upperx+lowerx)/2;
if (x[middle] > u)
upperx = middle;
else
lowerx = middle;
}
int lowery = 0;
int uppery = ysize-1;
while (uppery-lowery > 1) {
int middle = (uppery+lowery)/2;
if (y[middle] > v)
uppery = middle;
else
lowery = middle;
}
int lowerz = 0;
int upperz = zsize-1;
while (upperz-lowerz > 1) {
int middle = (upperz+lowerz)/2;
if (z[middle] > w)
upperz = middle;
else
lowerz = middle;
}
double deltax = x[upperx]-x[lowerx];
double deltay = y[uppery]-y[lowery];
double deltaz = z[upperz]-z[lowerz];
double da = (u-x[lowerx])/deltax;
double db = (v-y[lowery])/deltay;
double dc = (w-z[lowerz])/deltaz;
const vector<double>& coeff = c[lowerx+(xsize-1)*lowery+(xsize-1)*(ysize-1)*lowerz];
// Evaluate the spline to determine the value and gradients.
double value[] = {0, 0, 0, 0};
for (int i = 3; i >= 0; i--) {
for (int j = 0; j < 4; j++) {
int base = 4*i + 16*j;
value[j] = db*value[j] + ((coeff[base+3]*da + coeff[base+2])*da + coeff[base+1])*da + coeff[base];
}
}
return value[0] + dc*(value[1] + dc*(value[2] + dc*value[3]));
}
void SplineFitter::evaluate3DSplineDerivatives(const vector<double>& x, const vector<double>& y, const vector<double>& z, const vector<double>& values, const vector<vector<double> >& c, double u, double v, double w, double& dx, double& dy, double& dz) {
int xsize = x.size();
int ysize = y.size();
int zsize = z.size();
if (u < x[0] || u > x[xsize-1] || v < y[0] || v > y[ysize-1] || w < z[0] || w > z[zsize-1])
throw OpenMMException("evaluate3DSpline: specified point is outside the range defined by the spline");
// Perform a binary search to identify the interval containing the point to evaluate.
int lowerx = 0;
int upperx = xsize-1;
while (upperx-lowerx > 1) {
int middle = (upperx+lowerx)/2;
if (x[middle] > u)
upperx = middle;
else
lowerx = middle;
}
int lowery = 0;
int uppery = ysize-1;
while (uppery-lowery > 1) {
int middle = (uppery+lowery)/2;
if (y[middle] > v)
uppery = middle;
else
lowery = middle;
}
int lowerz = 0;
int upperz = zsize-1;
while (upperz-lowerz > 1) {
int middle = (upperz+lowerz)/2;
if (z[middle] > w)
upperz = middle;
else
lowerz = middle;
}
double deltax = x[upperx]-x[lowerx];
double deltay = y[uppery]-y[lowery];
double deltaz = z[upperz]-z[lowerz];
double da = (u-x[lowerx])/deltax;
double db = (v-y[lowery])/deltay;
double dc = (w-z[lowerz])/deltaz;
const vector<double>& coeff = c[lowerx+(xsize-1)*lowery+(xsize-1)*(ysize-1)*lowerz];
// Evaluate the spline to determine the derivatives.
double derivx[] = {0, 0, 0, 0};
double derivy[] = {0, 0, 0, 0};
double derivz[] = {0, 0, 0, 0};
for (int i = 3; i >= 0; i--) {
for (int j = 0; j < 4; j++) {
int base = 4*i + 16*j;
derivx[j] = db*derivx[j] + (3.0*coeff[base+3]*da + 2.0*coeff[base+2])*da + coeff[base+1];
derivz[j] = db*derivz[j] + ((coeff[base+3]*da + coeff[base+2])*da + coeff[base+1])*da + coeff[base];
base = i + 16*j;
derivy[j] = da*derivy[j] + (3.0*coeff[base+12]*db + 2.0*coeff[base+8])*db + coeff[base+4];
}
}
dx = derivx[0] + dc*(derivx[1] + dc*(derivx[2] + dc*derivx[3]));
dy = derivy[0] + dc*(derivy[1] + dc*(derivy[2] + dc*derivy[3]));
dz = derivz[1] + dc*(2.0*derivz[2] + 3.0*dc*derivz[3]);
dx /= deltax;
dy /= deltay;
dz /= deltaz;
}
openmmapi/src/TabulatedFunction.cpp 0 → 100644
View file @ cf8a03e8
/* -------------------------------------------------------------------------- *
* 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) 2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/TabulatedFunction.h"
#include "openmm/OpenMMException.h"
using namespace OpenMM;
using namespace std;
Continuous1DFunction::Continuous1DFunction(const vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("Continuous1DFunction: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("Continuous1DFunction: a tabulated function must have at least two points");
this->values = values;
this->min = min;
this->max = max;
}
void Continuous1DFunction::getFunctionParameters(vector<double>& values, double& min, double& max) const {
values = this->values;
min = this->min;
max = this->max;
}
void Continuous1DFunction::setFunctionParameters(const vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("Continuous1DFunction: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("Continuous1DFunction: a tabulated function must have at least two points");
this->values = values;
this->min = min;
this->max = max;
}
Continuous2DFunction::Continuous2DFunction(int xsize, int ysize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax) {
if (xsize < 2 || ysize < 2)
throw OpenMMException("Continuous2DFunction: must have at least two points along each axis");
if (values.size() != xsize*ysize)
throw OpenMMException("Continuous2DFunction: incorrect number of values");
if (xmax <= xmin)
throw OpenMMException("Continuous2DFunction: xmax <= xmin for a tabulated function.");
if (ymax <= ymin)
throw OpenMMException("Continuous2DFunction: ymax <= ymin for a tabulated function.");
this->values = values;
this->xsize = xsize;
this->ysize = ysize;
this->xmin = xmin;
this->xmax = xmax;
this->ymin = ymin;
this->ymax = ymax;
}
void Continuous2DFunction::getFunctionParameters(int& xsize, int& ysize, vector<double>& values, double& xmin, double& xmax, double& ymin, double& ymax) const {
values = this->values;
xsize = this->xsize;
ysize = this->ysize;
xmin = this->xmin;
xmax = this->xmax;
ymin = this->ymin;
ymax = this->ymax;
}
void Continuous2DFunction::setFunctionParameters(int xsize, int ysize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax) {
if (xsize < 2 || ysize < 2)
throw OpenMMException("Continuous2DFunction: must have at least two points along each axis");
if (values.size() != xsize*ysize)
throw OpenMMException("Continuous2DFunction: incorrect number of values");
if (xmax <= xmin)
throw OpenMMException("Continuous2DFunction: xmax <= xmin for a tabulated function.");
if (ymax <= ymin)
throw OpenMMException("Continuous2DFunction: ymax <= ymin for a tabulated function.");
this->values = values;
this->xsize = xsize;
this->ysize = ysize;
this->xmin = xmin;
this->xmax = xmax;
this->ymin = ymin;
this->ymax = ymax;
}
Continuous3DFunction::Continuous3DFunction(int xsize, int ysize, int zsize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax) {
if (xsize < 2 || ysize < 2 || zsize < 2)
throw OpenMMException("Continuous3DFunction: must have at least two points along each axis");
if (values.size() != xsize*ysize*zsize)
throw OpenMMException("Continuous3DFunction: incorrect number of values");
if (xmax <= xmin)
throw OpenMMException("Continuous3DFunction: xmax <= xmin for a tabulated function.");
if (ymax <= ymin)
throw OpenMMException("Continuous3DFunction: ymax <= ymin for a tabulated function.");
if (zmax <= zmin)
throw OpenMMException("Continuous3DFunction: zmax <= zmin for a tabulated function.");
this->values = values;
this->xsize = xsize;
this->ysize = ysize;
this->zsize = zsize;
this->xmin = xmin;
this->xmax = xmax;
this->ymin = ymin;
this->ymax = ymax;
this->zmin = zmin;
this->zmax = zmax;
}
void Continuous3DFunction::getFunctionParameters(int& xsize, int& ysize, int& zsize, vector<double>& values, double& xmin, double& xmax, double& ymin, double& ymax, double& zmin, double& zmax) const {
values = this->values;
xsize = this->xsize;
ysize = this->ysize;
zsize = this->zsize;
xmin = this->xmin;
xmax = this->xmax;
ymin = this->ymin;
ymax = this->ymax;
zmin = this->zmin;
zmax = this->zmax;
}
void Continuous3DFunction::setFunctionParameters(int xsize, int ysize, int zsize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax) {
if (xsize < 2 || ysize < 2 || zsize < 2)
throw OpenMMException("Continuous3DFunction: must have at least two points along each axis");
if (values.size() != xsize*ysize*zsize)
throw OpenMMException("Continuous3DFunction: incorrect number of values");
if (xmax <= xmin)
throw OpenMMException("Continuous3DFunction: xmax <= xmin for a tabulated function.");
if (ymax <= ymin)
throw OpenMMException("Continuous3DFunction: ymax <= ymin for a tabulated function.");
if (zmax <= zmin)
throw OpenMMException("Continuous3DFunction: zmax <= zmin for a tabulated function.");
this->values = values;
this->xsize = xsize;
this->ysize = ysize;
this->zsize = zsize;
this->xmin = xmin;
this->xmax = xmax;
this->ymin = ymin;
this->ymax = ymax;
this->zmin = zmin;
this->zmax = zmax;
}
Discrete1DFunction::Discrete1DFunction(const vector<double>& values) {
this->values = values;
}
void Discrete1DFunction::getFunctionParameters(vector<double>& values) const {
values = this->values;
}
void Discrete1DFunction::setFunctionParameters(const vector<double>& values) {
this->values = values;
}
Discrete2DFunction::Discrete2DFunction(int xsize, int ysize, const vector<double>& values) {
if (values.size() != xsize*ysize)
throw OpenMMException("Discrete2DFunction: incorrect number of values");
this->xsize = xsize;
this->ysize = ysize;
this->values = values;
}
void Discrete2DFunction::getFunctionParameters(int& xsize, int& ysize, vector<double>& values) const {
xsize = this->xsize;
ysize = this->ysize;
values = this->values;
}
void Discrete2DFunction::setFunctionParameters(int xsize, int ysize, const vector<double>& values) {
if (values.size() != xsize*ysize)
throw OpenMMException("Discrete2DFunction: incorrect number of values");
this->xsize = xsize;
this->ysize = ysize;
this->values = values;
}
Discrete3DFunction::Discrete3DFunction(int xsize, int ysize, int zsize, const vector<double>& values) {
if (values.size() != xsize*ysize*zsize)
throw OpenMMException("Discrete3DFunction: incorrect number of values");
this->xsize = xsize;
this->ysize = ysize;
this->zsize = zsize;
this->values = values;
}
void Discrete3DFunction::getFunctionParameters(int& xsize, int& ysize, int& zsize, vector<double>& values) const {
xsize = this->xsize;
ysize = this->ysize;
zsize = this->zsize;
values = this->values;
}
void Discrete3DFunction::setFunctionParameters(int xsize, int ysize, int zsize, const vector<double>& values) {
if (values.size() != xsize*ysize*zsize)
throw OpenMMException("Discrete3DFunction: incorrect number of values");
this->xsize = xsize;
this->ysize = ysize;
this->zsize = zsize;
this->values = values;
}
openmmapi/src/ThreadPool.cpp
View file @ cf8a03e8
......@@ -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-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -60,8 +60,10 @@ static void* threadBody(void* args) {
return 0;
}
ThreadPool::ThreadPool() {
numThreads = getNumProcessors();
ThreadPool::ThreadPool(int numThreads) {
if (numThreads <= 0)
numThreads = getNumProcessors();
this->numThreads = numThreads;
pthread_cond_init(&startCondition, NULL);
pthread_cond_init(&endCondition, NULL);
pthread_mutex_init(&lock, NULL);
......
platforms/cpu/CMakeLists.txt
View file @ cf8a03e8
......@@ -32,9 +32,9 @@ SET(STATIC_TARGET ${OPENMMCPU_LIBRARY_NAME}_static)
# Ensure that debug libraries have "_d" appended to their names.
# CMake gets this right on Windows automatically with this definition.
IF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
IF (MSVC)
SET(CMAKE_DEBUG_POSTFIX "_d" CACHE INTERNAL "" FORCE)
ENDIF (${CMAKE_GENERATOR} MATCHES "Visual Studio")
ENDIF (MSVC)
# But on Unix or Cygwin we have to add the suffix manually
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
......@@ -92,3 +92,7 @@ FILE(GLOB CORE_HEADERS include/*.h)
INSTALL_FILES(/include/openmm/cpu FILES ${CORE_HEADERS})
SUBDIRS (sharedTarget)
IF(OPENMM_BUILD_STATIC_LIB)
SUBDIRS (staticTarget)
ENDIF(OPENMM_BUILD_STATIC_LIB)
platforms/cpu/include/CpuBondForce.h 0 → 100644
View file @ cf8a03e8
#ifndef OPENMM_CPUBONDFORCE_H_
#define OPENMM_CPUBONDFORCE_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2014 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 "ReferenceBondIxn.h"
#include "windowsExportCpu.h"
#include "openmm/internal/ThreadPool.h"
#include <list>
#include <set>
#include <vector>
namespace OpenMM {
/**
* This class parallelizes the calculation of bonded forces.
*/
class OPENMM_EXPORT_CPU CpuBondForce {
public:
class ComputeForceTask;
CpuBondForce();
/**
* Analyze the set of bonds and decide which to compute with each thread.
*/
void initialize(int numAtoms, int numBonds, int numAtomsPerBond, int** bondAtoms, ThreadPool& threads);
/**
* Compute the forces from all bonds.
*/
void calculateForce(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** parameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn);
/**
* This routine contains the code executed by each thread.
*/
void threadComputeForce(ThreadPool& threads, int threadIndex, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** parameters,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn);
private:
bool canAssignBond(int bond, int thread, std::vector<int>& atomThread);
void assignBond(int bond, int thread, std::vector<int>& atomThread, std::vector<int>& bondThread, std::vector<std::set<int> >& atomBonds, std::list<int>& candidateBonds);
int numBonds, numAtomsPerBond;
int** bondAtoms;
ThreadPool* threads;
std::vector<std::vector<int> > threadBonds;
std::vector<int> extraBonds;
};
} // namespace OpenMM
#endif /*OPENMM_CPUBONDFORCE_H_*/
platforms/cpu/include/CpuKernels.h
View file @ cf8a03e8
......@@ -32,6 +32,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "CpuBondForce.h"
#include "CpuGBSAOBCForce.h"
#include "CpuLangevinDynamics.h"
#include "CpuNeighborList.h"
......@@ -49,6 +50,7 @@ namespace OpenMM {
*/
class CpuCalcForcesAndEnergyKernel : public CalcForcesAndEnergyKernel {
public:
class InitForceTask;
class SumForceTask;
CpuCalcForcesAndEnergyKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data, ContextImpl& context);
/**
......@@ -85,6 +87,86 @@ private:
Kernel referenceKernel;
};
/**
* This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
class CpuCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
public:
CpuCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) :
CalcPeriodicTorsionForceKernel(name, platform), data(data), torsionIndexArray(NULL), torsionParamArray(NULL) {
}
~CpuCalcPeriodicTorsionForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the PeriodicTorsionForce this kernel will be used for
*/
void initialize(const System& system, const PeriodicTorsionForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the PeriodicTorsionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force);
private:
CpuPlatform::PlatformData& data;
int numTorsions;
int **torsionIndexArray;
RealOpenMM **torsionParamArray;
CpuBondForce bondForce;
};
/**
* This kernel is invoked by RBTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
class CpuCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
public:
CpuCalcRBTorsionForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) :
CalcRBTorsionForceKernel(name, platform), data(data), torsionIndexArray(NULL), torsionParamArray(NULL) {
}
~CpuCalcRBTorsionForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the RBTorsionForce this kernel will be used for
*/
void initialize(const System& system, const RBTorsionForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the RBTorsionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const RBTorsionForce& force);
private:
CpuPlatform::PlatformData& data;
int numTorsions;
int **torsionIndexArray;
RealOpenMM **torsionParamArray;
CpuBondForce bondForce;
};
/**
* This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
*/
......
platforms/cpu/include/CpuPlatform.h
View file @ cf8a03e8
......@@ -55,27 +55,37 @@ public:
return name;
}
double getSpeed() const;
const std::string& getPropertyValue(const Context& context, const std::string& property) const;
bool supportsDoublePrecision() const;
static bool isProcessorSupported();
void contextCreated(ContextImpl& context, const std::map<std::string, std::string>& properties) const;
void contextDestroyed(ContextImpl& context) const;
/**
* This is the name of the parameter for selecting the number of threads to use.
*/
static const std::string& CpuThreads() {
static const std::string key = "CpuThreads";
return key;
}
/**
* We cannot use the standard mechanism for platform data, because that is already used by the superclass.
* Instead, we maintain a table of ContextImpls to PlatformDatas.
*/
static PlatformData& getPlatformData(ContextImpl& context);
static const PlatformData& getPlatformData(const ContextImpl& context);
private:
static std::map<ContextImpl*, PlatformData*> contextData;
static std::map<const ContextImpl*, PlatformData*> contextData;
};
class CpuPlatform::PlatformData {
public:
PlatformData(int numParticles);
PlatformData(int numParticles, int numThreads);
AlignedArray<float> posq;
std::vector<AlignedArray<float> > threadForce;
ThreadPool threads;
bool isPeriodic;
CpuRandom random;
std::map<std::string, std::string> propertyValues;
};
} // namespace OpenMM
......
platforms/cpu/src/CpuBondForce.cpp 0 → 100644
View file @ cf8a03e8
/* -------------------------------------------------------------------------- *
* 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) 2014 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 "CpuBondForce.h"
#include "openmm/OpenMMException.h"
using namespace OpenMM;
using namespace std;
class CpuBondForce::ComputeForceTask : public ThreadPool::Task {
public:
ComputeForceTask(CpuBondForce& owner, vector<RealVec>& atomCoordinates, RealOpenMM** parameters, vector<RealVec>& forces,
vector<RealOpenMM>& threadEnergy, RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn) : owner(owner), atomCoordinates(atomCoordinates),
parameters(parameters), forces(forces), threadEnergy(threadEnergy), totalEnergy(totalEnergy), referenceBondIxn(referenceBondIxn) {
}
void execute(ThreadPool& threads, int threadIndex) {
RealOpenMM* energy = (totalEnergy == NULL ? NULL : &threadEnergy[threadIndex]);
owner.threadComputeForce(threads, threadIndex, atomCoordinates, parameters, forces, energy, referenceBondIxn);
}
CpuBondForce& owner;
vector<RealVec>& atomCoordinates;
RealOpenMM** parameters;
vector<RealVec>& forces;
vector<RealOpenMM>& threadEnergy;
RealOpenMM* totalEnergy;
ReferenceBondIxn& referenceBondIxn;
};
CpuBondForce::CpuBondForce() {
}
void CpuBondForce::initialize(int numAtoms, int numBonds, int numAtomsPerBond, int** bondAtoms, ThreadPool& threads) {
this->numBonds = numBonds;
this->numAtomsPerBond = numAtomsPerBond;
this->bondAtoms = bondAtoms;
this->threads = &threads;
int numThreads = threads.getNumThreads();
int targetBondsPerThread = numBonds/numThreads;
// Record the bonds that include each atom.
vector<set<int> > atomBonds(numAtoms);
for (int bond = 0; bond < numBonds; bond++) {
for (int i = 0; i < numAtomsPerBond; i++)
atomBonds[bondAtoms[bond][i]].insert(bond);
}
// Divide bonds into groups.
vector<int> atomThread(numAtoms, -1);
vector<int> bondThread(numBonds, -1);
threadBonds.resize(numThreads);
int numProcessed = 0;
int thread = 0;
list<int> candidateBonds;
while (thread < numThreads) {
// Find the next unassigned bond.
while (numProcessed < numBonds && bondThread[numProcessed] != -1)
numProcessed++;
if (numProcessed == numBonds)
break; // We've gone through the whole list of bonds.
// See if this bond can be assigned to this thread.
if (!canAssignBond(numProcessed, thread, atomThread)) {
numProcessed++;
continue;
}
// Assign this bond to the thread.
assignBond(numProcessed++, thread, atomThread, bondThread, atomBonds, candidateBonds);
// Assign additional bonds that have been identified as involving atoms assigned to this thread.
while (!candidateBonds.empty() && threadBonds[thread].size() < targetBondsPerThread) {
int bond = *candidateBonds.begin();
if (bondThread[bond] == -1 && canAssignBond(bond, thread, atomThread))
assignBond(bond, thread, atomThread, bondThread, atomBonds, candidateBonds);
candidateBonds.pop_front();
}
// If we have assigned enough bonds to this thread, move on to the next one.
if (threadBonds[thread].size() >= targetBondsPerThread) {
candidateBonds.clear();
thread++;
}
}
// Look through the remaining bonds and see whether any of them can be assigned.
candidateBonds.clear();
for (int bond = 0; bond < numBonds; bond++) {
if (bondThread[bond] == -1) {
// See whether this bond can be assigned to a thread.
bool canAssign = true;
int assignment = -1;
for (int i = 0; i < numAtomsPerBond; i++) {
int thread = atomThread[bondAtoms[bond][i]];
if (thread == -1 || thread == assignment)
continue;
if (assignment == -1)
assignment = thread;
else {
canAssign = false;
break;
}
}
if (canAssign) {
// Assign this bond to a thread.
if (assignment == -1)
assignment = numThreads-1;
assignBond(bond, assignment, atomThread, bondThread, atomBonds, candidateBonds);
}
else {
// Add it to the list of "extra" bonds.
extraBonds.push_back(bond);
}
}
}
}
bool CpuBondForce::canAssignBond(int bond, int thread, vector<int>& atomThread) {
for (int i = 0; i < numAtomsPerBond; i++) {
int atom = bondAtoms[bond][i];
if (atomThread[atom] != -1 && atomThread[atom] != thread)
return false;
}
return true;
}
void CpuBondForce::assignBond(int bond, int thread, vector<int>& atomThread, vector<int>& bondThread, vector<set<int> >& atomBonds, list<int>& candidateBonds) {
// Assign the bond to a thread.
bondThread[bond] = thread;
threadBonds[thread].push_back(bond);
// Mark every atom in this bond as also belonging to the thread, and add all of their
// bonds to the list of candidates.
for (int i = 0; i < numAtomsPerBond; i++) {
int& atom = atomThread[bondAtoms[bond][i]];
if (atom == thread)
continue;
if (atom != -1)
throw OpenMMException("CpuBondForce: Internal error: atoms assigned to threads incorrectly");
atom = thread;
for (set<int>::const_iterator iter = atomBonds[atom].begin(); iter != atomBonds[atom].end(); ++iter)
candidateBonds.push_back(*iter);
}
}
void CpuBondForce::calculateForce(vector<RealVec>& atomCoordinates, RealOpenMM** parameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn) {
// Have the worker threads compute their forces.
vector<RealOpenMM> threadEnergy(threads->getNumThreads(), 0);
ComputeForceTask task(*this, atomCoordinates, parameters, forces, threadEnergy, totalEnergy, referenceBondIxn);
threads->execute(task);
threads->waitForThreads();
// Compute any "extra" bonds.
for (int i = 0; i < extraBonds.size(); i++) {
int bond = extraBonds[i];
referenceBondIxn.calculateBondIxn(bondAtoms[bond], atomCoordinates, parameters[bond], forces, totalEnergy);
}
// Compute the total energy.
if (totalEnergy != NULL)
for (int i = 0; i < threads->getNumThreads(); i++)
*totalEnergy += threadEnergy[i];
}
void CpuBondForce::threadComputeForce(ThreadPool& threads, int threadIndex, vector<RealVec>& atomCoordinates, RealOpenMM** parameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy, ReferenceBondIxn& referenceBondIxn) {
vector<int>& bonds = threadBonds[threadIndex];
int numBonds = bonds.size();
for (int i = 0; i < numBonds; i++) {
int bond = bonds[i];
referenceBondIxn.calculateBondIxn(bondAtoms[bond], atomCoordinates, parameters[bond], forces, totalEnergy);
}
}
\ No newline at end of file
platforms/cpu/src/CpuGBSAOBCForce.cpp
View file @ cf8a03e8
......@@ -404,10 +404,10 @@ void CpuGBSAOBCForce::getDeltaR(const fvec4& posI, const fvec4& x, const fvec4&
fvec4 CpuGBSAOBCForce::fastLog(const fvec4& x) {
// Evaluate log(x) using a lookup table for speed.
if (any((x < TABLE_MIN) | (x >= TABLE_MAX)))
return fvec4(logf(x[0]), logf(x[1]), logf(x[2]), logf(x[3]));
fvec4 x1 = (x-TABLE_MIN)*logDXInv;
ivec4 index = floor(x1);
if (any((index < 0) | (index >= NUM_TABLE_POINTS)))
return fvec4(logf(x[0]), logf(x[1]), logf(x[2]), logf(x[3]));
fvec4 coeff2 = x1-index;
fvec4 coeff1 = 1.0f-coeff2;
fvec4 t1(&logTable[index[0]]);
......
platforms/cpu/src/CpuKernelFactory.cpp
View file @ cf8a03e8
......@@ -41,6 +41,10 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform&
CpuPlatform::PlatformData& data = CpuPlatform::getPlatformData(context);
if (name == CalcForcesAndEnergyKernel::Name())
return new CpuCalcForcesAndEnergyKernel(name, platform, data, context);
if (name == CalcPeriodicTorsionForceKernel::Name())
return new CpuCalcPeriodicTorsionForceKernel(name, platform, data);
if (name == CalcRBTorsionForceKernel::Name())
return new CpuCalcRBTorsionForceKernel(name, platform, data);
if (name == CalcNonbondedForceKernel::Name())
return new CpuCalcNonbondedForceKernel(name, platform, data);
if (name == CalcGBSAOBCForceKernel::Name())
......
platforms/cpu/src/CpuKernels.cpp
View file @ cf8a03e8
......@@ -35,6 +35,8 @@
#include "ReferenceKernelFactory.h"
#include "ReferenceKernels.h"
#include "ReferenceLJCoulomb14.h"
#include "ReferenceProperDihedralBond.h"
#include "ReferenceRbDihedralBond.h"
#include "openmm/Context.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
......@@ -130,6 +132,46 @@ public:
CpuPlatform::PlatformData& data;
};
class CpuCalcForcesAndEnergyKernel::InitForceTask : public ThreadPool::Task {
public:
InitForceTask(int numParticles, ContextImpl& context, CpuPlatform::PlatformData& data) : numParticles(numParticles), context(context), data(data) {
}
void execute(ThreadPool& threads, int threadIndex) {
// Convert the positions to single precision and apply periodic boundary conditions
AlignedArray<float>& posq = data.posq;
vector<RealVec>& posData = extractPositions(context);
RealVec boxSize = extractBoxSize(context);
double invBoxSize[3] = {1/boxSize[0], 1/boxSize[1], 1/boxSize[2]};
int numParticles = context.getSystem().getNumParticles();
int numThreads = threads.getNumThreads();
int start = threadIndex*numParticles/numThreads;
int end = (threadIndex+1)*numParticles/numThreads;
if (data.isPeriodic)
for (int i = start; i < end; i++)
for (int j = 0; j < 3; j++) {
RealOpenMM x = posData[i][j];
double base = floor(x*invBoxSize[j])*boxSize[j];
posq[4*i+j] = (float) (x-base);
}
else
for (int i = start; i < end; i++) {
posq[4*i] = (float) posData[i][0];
posq[4*i+1] = (float) posData[i][1];
posq[4*i+2] = (float) posData[i][2];
}
// Clear the forces.
fvec4 zero(0.0f);
for (int j = 0; j < numParticles; j++)
zero.store(&data.threadForce[threadIndex][j*4]);
}
int numParticles;
ContextImpl& context;
CpuPlatform::PlatformData& data;
};
CpuCalcForcesAndEnergyKernel::CpuCalcForcesAndEnergyKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data, ContextImpl& context) :
CalcForcesAndEnergyKernel(name, platform), data(data) {
// Create a Reference platform version of this kernel.
......@@ -145,33 +187,11 @@ void CpuCalcForcesAndEnergyKernel::initialize(const System& system) {
void CpuCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups) {
referenceKernel.getAs<ReferenceCalcForcesAndEnergyKernel>().beginComputation(context, includeForce, includeEnergy, groups);
// Convert the positions to single precision and apply periodic boundary conditions
// Convert positions to single precision and clear the forces.
AlignedArray<float>& posq = data.posq;
vector<RealVec>& posData = extractPositions(context);
RealVec boxSize = extractBoxSize(context);
double invBoxSize[3] = {1/boxSize[0], 1/boxSize[1], 1/boxSize[2]};
int numParticles = context.getSystem().getNumParticles();
if (data.isPeriodic)
for (int i = 0; i < numParticles; i++)
for (int j = 0; j < 3; j++) {
RealOpenMM x = posData[i][j];
double base = floor(x*invBoxSize[j])*boxSize[j];
posq[4*i+j] = (float) (x-base);
}
else
for (int i = 0; i < numParticles; i++) {
posq[4*i] = (float) posData[i][0];
posq[4*i+1] = (float) posData[i][1];
posq[4*i+2] = (float) posData[i][2];
}
// Clear the forces.
fvec4 zero(0.0f);
for (int i = 0; i < (int) data.threadForce.size(); i++)
for (int j = 0; j < numParticles; j++)
zero.store(&data.threadForce[i][j*4]);
InitForceTask task(context.getSystem().getNumParticles(), context, data);
data.threads.execute(task);
data.threads.waitForThreads();
}
double CpuCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups) {
......@@ -183,6 +203,134 @@ double CpuCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, boo
return referenceKernel.getAs<ReferenceCalcForcesAndEnergyKernel>().finishComputation(context, includeForce, includeEnergy, groups);
}
CpuCalcPeriodicTorsionForceKernel::~CpuCalcPeriodicTorsionForceKernel() {
if (torsionIndexArray != NULL) {
for (int i = 0; i < numTorsions; i++) {
delete[] torsionIndexArray[i];
delete[] torsionParamArray[i];
}
delete[] torsionIndexArray;
delete[] torsionParamArray;
}
}
void CpuCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
numTorsions = force.getNumTorsions();
torsionIndexArray = new int*[numTorsions];
for (int i = 0; i < numTorsions; i++)
torsionIndexArray[i] = new int[4];
torsionParamArray = new RealOpenMM*[numTorsions];
for (int i = 0; i < numTorsions; i++)
torsionParamArray[i] = new RealOpenMM[3];
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4, periodicity;
double phase, k;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, k);
torsionIndexArray[i][0] = particle1;
torsionIndexArray[i][1] = particle2;
torsionIndexArray[i][2] = particle3;
torsionIndexArray[i][3] = particle4;
torsionParamArray[i][0] = (RealOpenMM) k;
torsionParamArray[i][1] = (RealOpenMM) phase;
torsionParamArray[i][2] = (RealOpenMM) periodicity;
}
bondForce.initialize(system.getNumParticles(), numTorsions, 4, torsionIndexArray, data.threads);
}
double CpuCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
RealOpenMM energy = 0;
ReferenceProperDihedralBond periodicTorsionBond;
bondForce.calculateForce(posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, periodicTorsionBond);
return energy;
}
void CpuCalcPeriodicTorsionForceKernel::copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force) {
if (numTorsions != force.getNumTorsions())
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
// Record the values.
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4, periodicity;
double phase, k;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, k);
if (particle1 != torsionIndexArray[i][0] || particle2 != torsionIndexArray[i][1] || particle3 != torsionIndexArray[i][2] || particle4 != torsionIndexArray[i][3])
throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
torsionParamArray[i][0] = (RealOpenMM) k;
torsionParamArray[i][1] = (RealOpenMM) phase;
torsionParamArray[i][2] = (RealOpenMM) periodicity;
}
}
CpuCalcRBTorsionForceKernel::~CpuCalcRBTorsionForceKernel() {
if (torsionIndexArray != NULL) {
for (int i = 0; i < numTorsions; i++) {
delete[] torsionIndexArray[i];
delete[] torsionParamArray[i];
}
delete[] torsionIndexArray;
delete[] torsionParamArray;
}
}
void CpuCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
numTorsions = force.getNumTorsions();
torsionIndexArray = new int*[numTorsions];
for (int i = 0; i < numTorsions; i++)
torsionIndexArray[i] = new int[4];
torsionParamArray = new RealOpenMM*[numTorsions];
for (int i = 0; i < numTorsions; i++)
torsionParamArray[i] = new RealOpenMM[6];
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4;
double c0, c1, c2, c3, c4, c5;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
torsionIndexArray[i][0] = particle1;
torsionIndexArray[i][1] = particle2;
torsionIndexArray[i][2] = particle3;
torsionIndexArray[i][3] = particle4;
torsionParamArray[i][0] = (RealOpenMM) c0;
torsionParamArray[i][1] = (RealOpenMM) c1;
torsionParamArray[i][2] = (RealOpenMM) c2;
torsionParamArray[i][3] = (RealOpenMM) c3;
torsionParamArray[i][4] = (RealOpenMM) c4;
torsionParamArray[i][5] = (RealOpenMM) c5;
}
bondForce.initialize(system.getNumParticles(), numTorsions, 4, torsionIndexArray, data.threads);
}
double CpuCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
RealOpenMM energy = 0;
ReferenceRbDihedralBond rbTorsionBond;
bondForce.calculateForce(posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, rbTorsionBond);
return energy;
}
void CpuCalcRBTorsionForceKernel::copyParametersToContext(ContextImpl& context, const RBTorsionForce& force) {
if (numTorsions != force.getNumTorsions())
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
// Record the values.
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4;
double c0, c1, c2, c3, c4, c5;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
if (particle1 != torsionIndexArray[i][0] || particle2 != torsionIndexArray[i][1] || particle3 != torsionIndexArray[i][2] || particle4 != torsionIndexArray[i][3])
throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
torsionParamArray[i][0] = (RealOpenMM) c0;
torsionParamArray[i][1] = (RealOpenMM) c1;
torsionParamArray[i][2] = (RealOpenMM) c2;
torsionParamArray[i][3] = (RealOpenMM) c3;
torsionParamArray[i][4] = (RealOpenMM) c4;
torsionParamArray[i][5] = (RealOpenMM) c5;
}
}
class CpuCalcNonbondedForceKernel::PmeIO : public CalcPmeReciprocalForceKernel::IO {
public:
PmeIO(float* posq, float* force, int numParticles) : posq(posq), force(force), numParticles(numParticles) {
......
platforms/cpu/src/CpuPlatform.cpp
View file @ cf8a03e8
......@@ -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-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -35,25 +35,52 @@
#include "CpuSETTLE.h"
#include "ReferenceConstraints.h"
#include "openmm/internal/hardware.h"
#include <sstream>
using namespace OpenMM;
using namespace std;
#ifdef OPENMM_CPU_BUILDING_STATIC_LIBRARY
extern "C" void registerCpuPlatform() {
if (CpuPlatform::isProcessorSupported())
Platform::registerPlatform(new CpuPlatform());
}
#else
extern "C" OPENMM_EXPORT_CPU void registerPlatforms() {
// Only register this platform if the CPU supports SSE 4.1.
if (CpuPlatform::isProcessorSupported())
Platform::registerPlatform(new CpuPlatform());
}
#endif
map<ContextImpl*, CpuPlatform::PlatformData*> CpuPlatform::contextData;
map<const ContextImpl*, CpuPlatform::PlatformData*> CpuPlatform::contextData;
CpuPlatform::CpuPlatform() {
CpuKernelFactory* factory = new CpuKernelFactory();
registerKernelFactory(CalcForcesAndEnergyKernel::Name(), factory);
registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
platformProperties.push_back(CpuThreads());
int threads = getNumProcessors();
char* threadsEnv = getenv("OPENMM_CPU_THREADS");
if (threadsEnv != NULL)
stringstream(threadsEnv) >> threads;
stringstream defaultThreads;
defaultThreads << threads;
setPropertyDefaultValue(CpuThreads(), defaultThreads.str());
}
const string& CpuPlatform::getPropertyValue(const Context& context, const string& property) const {
const ContextImpl& impl = getContextImpl(context);
const PlatformData& data = getPlatformData(impl);
map<string, string>::const_iterator value = data.propertyValues.find(property);
if (value != data.propertyValues.end())
return value->second;
return ReferencePlatform::getPropertyValue(context, property);
}
double CpuPlatform::getSpeed() const {
......@@ -78,7 +105,11 @@ bool CpuPlatform::isProcessorSupported() {
void CpuPlatform::contextCreated(ContextImpl& context, const map<string, string>& properties) const {
ReferencePlatform::contextCreated(context, properties);
PlatformData* data = new PlatformData(context.getSystem().getNumParticles());
const string& threadsPropValue = (properties.find(CpuThreads()) == properties.end() ?
getPropertyDefaultValue(CpuThreads()) : properties.find(CpuThreads())->second);
int numThreads;
stringstream(threadsPropValue) >> numThreads;
PlatformData* data = new PlatformData(context.getSystem().getNumParticles(), numThreads);
contextData[&context] = data;
ReferenceConstraints& constraints = *(ReferenceConstraints*) reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData())->constraints;
if (constraints.settle != NULL) {
......@@ -98,10 +129,17 @@ CpuPlatform::PlatformData& CpuPlatform::getPlatformData(ContextImpl& context) {
return *contextData[&context];
}
CpuPlatform::PlatformData::PlatformData(int numParticles) : posq(4*numParticles) {
int numThreads = threads.getNumThreads();
const CpuPlatform::PlatformData& CpuPlatform::getPlatformData(const ContextImpl& context) {
return *contextData[&context];
}
CpuPlatform::PlatformData::PlatformData(int numParticles, int numThreads) : posq(4*numParticles), threads(numThreads) {
numThreads = threads.getNumThreads();
threadForce.resize(numThreads);
for (int i = 0; i < numThreads; i++)
threadForce[i].resize(4*numParticles);
isPeriodic = false;
stringstream threadsProperty;
threadsProperty << numThreads;
propertyValues[CpuThreads()] = threadsProperty.str();
}
platforms/cpu/src/CpuSETTLE.cpp
View file @ cf8a03e8
......@@ -41,7 +41,9 @@ public:
inverseMasses(inverseMasses), tolerance(tolerance), threadSettle(threadSettle) {
}
void execute(ThreadPool& threads, int threadIndex) {
threadSettle[threadIndex]->apply(atomCoordinates, atomCoordinatesP, inverseMasses, tolerance);
if (threadIndex < threadSettle.size()) {
threadSettle[threadIndex]->apply(atomCoordinates, atomCoordinatesP, inverseMasses, tolerance);
}
}
vector<OpenMM::RealVec>& atomCoordinates;
vector<OpenMM::RealVec>& atomCoordinatesP;
......
platforms/cpu/staticTarget/CMakeLists.txt 0 → 100644
View file @ cf8a03e8
FOREACH(file ${SOURCE_FILES})
IF (file MATCHES ".*Vec8.*")
IF (MSVC)
SET_SOURCE_FILES_PROPERTIES(${file} PROPERTIES COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} /arch:AVX /D__AVX__")
ELSE (MSVC)
SET_SOURCE_FILES_PROPERTIES(${file} PROPERTIES COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -msse4.1 -mavx")
ENDIF (MSVC)
ELSE (file MATCHES ".*Vec8.*")
IF (NOT MSVC)
SET_SOURCE_FILES_PROPERTIES(${file} PROPERTIES COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -msse4.1")
ENDIF (NOT MSVC)
ENDIF (file MATCHES ".*Vec8.*")
ENDFOREACH(file)
ADD_LIBRARY(${STATIC_TARGET} STATIC ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(MAIN_OPENMM_LIB ${OPENMM_LIBRARY_NAME}_d)
ELSE (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(MAIN_OPENMM_LIB ${OPENMM_LIBRARY_NAME})
ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
TARGET_LINK_LIBRARIES(${STATIC_TARGET} ${MAIN_OPENMM_LIB} ${PTHREADS_LIB_STATIC})
#-DPTW32_STATIC_LIB only works for the windows pthreads.
SET_TARGET_PROPERTIES(${STATIC_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_CPU_BUILDING_STATIC_LIBRARY -DPTW32_STATIC_LIB")
INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${STATIC_TARGET})
platforms/cpu/tests/TestCpuGBSAOBCForce.cpp
View file @ cf8a03e8
......@@ -223,6 +223,10 @@ void testForce(int numParticles, NonbondedForce::NonbondedMethod method, GBSAOBC
int main() {
try {
if (!CpuPlatform::isProcessorSupported()) {
cout << "CPU is not supported. Exiting." << endl;
return 0;
}
testSingleParticle();
testCutoffAndPeriodic();
for (int i = 5; i < 11; i++) {
......
platforms/cpu/tests/TestCpuLangevinIntegrator.cpp
View file @ cf8a03e8
......@@ -264,6 +264,10 @@ void testRandomSeed() {
int main() {
try {
if (!CpuPlatform::isProcessorSupported()) {
cout << "CPU is not supported. Exiting." << endl;
return 0;
}
testSingleBond();
testTemperature();
testConstraints();
......
platforms/cpu/tests/TestCpuPeriodicTorsionForce.cpp 0 → 100644
View file @ cf8a03e8
/* -------------------------------------------------------------------------- *
* 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-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CPU implementation of PeriodicTorsionForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CpuPlatform.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CpuPlatform platform;
const double TOL = 1e-5;
void testPeriodicTorsions() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
PeriodicTorsionForce* forceField = new PeriodicTorsionForce();
forceField->addTorsion(0, 1, 2, 3, 2, PI_M/3, 1.1);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(1, 0, 2);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double torque = -2*1.1*std::sin(2*PI_M/3);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), TOL);
}
// Try changing the torsion parameters and make sure it's still correct.
forceField->setTorsionParameters(0, 0, 1, 2, 3, 3, PI_M/3.2, 1.3);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double dtheta = (3*PI_M/2)-(PI_M/3.2);
double torque = -3*1.3*std::sin(dtheta);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(1.3*(1+std::cos(dtheta)), state.getPotentialEnergy(), TOL);
}
}
void testParallelComputation() {
System system;
const int numParticles = 200;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
PeriodicTorsionForce* force = new PeriodicTorsionForce();
for (int i = 3; i < numParticles; i++)
force->addTorsion(i-3, i-2, i-1, i, 2, 1.1, i);
system.addForce(force);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(i, i%2, i%3);
VerletIntegrator integrator1(0.01);
ReferencePlatform reference;
Context context1(system, integrator1, reference);
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
VerletIntegrator integrator2(0.01);
Context context2(system, integrator2, platform);
context2.setPositions(positions);
State state2 = context2.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
}
int main(int argc, char* argv[]) {
try {
testPeriodicTorsions();
testParallelComputation();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment