Skip to content

GitLab

Commit 3279c73a authored by peastman's avatar peastman
Browse files

Created OpenCL and CUDA implementations of Continuous3DFunction

parent c9c218ac
Hide whitespace changes
Inline Side-by-side
Showing with 680 additions and 43 deletions
openmmapi/src/CustomCompoundBondForceImpl.cpp
View file @ 3279c73a
...@@ -147,7 +147,7 @@ ParsedExpression CustomCompoundBondForceImpl::prepareExpression(const CustomComp ...@@ -147,7 +147,7 @@ ParsedExpression CustomCompoundBondForceImpl::prepareExpression(const CustomComp
ExpressionTreeNode CustomCompoundBondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms, ExpressionTreeNode CustomCompoundBondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms,
map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) { map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
const Operation& op = node.getOperation(); const Operation& op = node.getOperation();
if (op.getId() != Operation::CUSTOM || op.getNumArguments() < 2) if (op.getId() != Operation::CUSTOM || (op.getName() != "distance" && op.getName() != "angle" && op.getName() != "dihedral"))
{ {
// This is not an angle or dihedral, so process its children. // This is not an angle or dihedral, so process its children.
......
platforms/cuda/src/CudaExpressionUtilities.cpp
View file @ 3279c73a
...@@ -133,36 +133,100 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express ...@@ -133,36 +133,100 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express
out << "x = (x-" << paramsFloat[2] << ")*" << paramsFloat[6] << ";\n"; out << "x = (x-" << paramsFloat[2] << ")*" << paramsFloat[6] << ";\n";
out << "y = (y-" << paramsFloat[4] << ")*" << paramsFloat[7] << ";\n"; out << "y = (y-" << paramsFloat[4] << ")*" << paramsFloat[7] << ";\n";
out << "int s = min((int) floor(x), " << paramsInt[0] << ");\n"; out << "int s = min((int) floor(x), " << paramsInt[0] << ");\n";
out << "int t = min((int) floor(y), " << paramsInt[0] << ");\n"; out << "int t = min((int) floor(y), " << paramsInt[1] << ");\n";
out << "int coeffIndex = 4*(s+" << paramsInt[0] << "*t);\n"; out << "int coeffIndex = 4*(s+" << paramsInt[0] << "*t);\n";
out << "float4 c[4];\n"; out << "float4 c[4];\n";
out << "c[0] = " << functionNames[i].second << "[coeffIndex];\n"; for (int j = 0; j < 4; j++)
out << "c[1] = " << functionNames[i].second << "[coeffIndex+1];\n"; out << "c[" << j << "] = " << functionNames[i].second << "[coeffIndex+" << j << "];\n";
out << "c[2] = " << functionNames[i].second << "[coeffIndex+2];\n"; out << "real da = x-s;\n";
out << "c[3] = " << functionNames[i].second << "[coeffIndex+3];\n"; out << "real db = y-t;\n";
out << "real da = x-s;";
out << "real db = y-t;";
for (int j = 0; j < nodes.size(); j++) { for (int j = 0; j < nodes.size(); j++) {
const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder(); const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder();
if (derivOrder[0] == 0 && derivOrder[1] == 0) { if (derivOrder[0] == 0 && derivOrder[1] == 0) {
out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[3].w*db + c[3].z)*db + c[3].y)*db + c[3].x;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[3].w*db + c[3].z)*db + c[3].y)*db + c[3].x;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[2].w*db + c[2].z)*db + c[2].y)*db + c[2].x;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[2].w*db + c[2].z)*db + c[2].y)*db + c[2].x;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[1].w*db + c[1].z)*db + c[1].y)*db + c[1].x;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[1].w*db + c[1].z)*db + c[1].y)*db + c[1].x;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[0].w*db + c[0].z)*db + c[0].y)*db + c[0].x;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[0].w*db + c[0].z)*db + c[0].y)*db + c[0].x;\n";
} }
else if (derivOrder[0] == 1 && derivOrder[1] == 0) { else if (derivOrder[0] == 1 && derivOrder[1] == 0) {
out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].w*da + 2.0f*c[2].w)*da + c[1].w;"; out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].w*da + 2.0f*c[2].w)*da + c[1].w;\n";
out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].z*da + 2.0f*c[2].z)*da + c[1].z;"; out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].z*da + 2.0f*c[2].z)*da + c[1].z;\n";
out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].y*da + 2.0f*c[2].y)*da + c[1].y;"; out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].y*da + 2.0f*c[2].y)*da + c[1].y;\n";
out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].x*da + 2.0f*c[2].x)*da + c[1].x;"; out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].x*da + 2.0f*c[2].x)*da + c[1].x;\n";
out << nodeNames[j] << " *= " << paramsFloat[6] << ";"; out << nodeNames[j] << " *= " << paramsFloat[6] << ";\n";
} }
else if (derivOrder[0] == 0 && derivOrder[1] == 1) { else if (derivOrder[0] == 0 && derivOrder[1] == 1) {
out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[3].w*db + 2.0f*c[3].z)*db + c[3].y;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[3].w*db + 2.0f*c[3].z)*db + c[3].y;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[2].w*db + 2.0f*c[2].z)*db + c[2].y;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[2].w*db + 2.0f*c[2].z)*db + c[2].y;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[1].w*db + 2.0f*c[1].z)*db + c[1].y;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[1].w*db + 2.0f*c[1].z)*db + c[1].y;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[0].w*db + 2.0f*c[0].z)*db + c[0].y;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[0].w*db + 2.0f*c[0].z)*db + c[0].y;\n";
out << nodeNames[j] << " *= " << paramsFloat[7] << ";"; out << nodeNames[j] << " *= " << paramsFloat[7] << ";\n";
}
else
throw OpenMMException("Unsupported derivative order for Continuous2DFunction");
}
out << "}\n";
}
else if (dynamic_cast<const Continuous3DFunction*>(functions[i]) != NULL) {
out << "real x = " << getTempName(node.getChildren()[0], temps) << ";\n";
out << "real y = " << getTempName(node.getChildren()[1], temps) << ";\n";
out << "real z = " << getTempName(node.getChildren()[2], temps) << ";\n";
out << "if (x >= " << paramsFloat[3] << " && x <= " << paramsFloat[4] << " && y >= " << paramsFloat[5] << " && y <= " << paramsFloat[6] << " && z >= " << paramsFloat[7] << " && z <= " << paramsFloat[8] << ") {\n";
out << "x = (x-" << paramsFloat[3] << ")*" << paramsFloat[9] << ";\n";
out << "y = (y-" << paramsFloat[5] << ")*" << paramsFloat[10] << ";\n";
out << "z = (z-" << paramsFloat[7] << ")*" << paramsFloat[11] << ";\n";
out << "int s = min((int) floor(x), " << paramsInt[0] << ");\n";
out << "int t = min((int) floor(y), " << paramsInt[1] << ");\n";
out << "int u = min((int) floor(z), " << paramsInt[2] << ");\n";
out << "int coeffIndex = 16*(s+" << paramsInt[0] << "*(t+" << paramsInt[1] << "*u));\n";
out << "float4 c[16];\n";
for (int j = 0; j < 16; j++)
out << "c[" << j << "] = " << functionNames[i].second << "[coeffIndex+" << j << "];\n";
out << "real da = x-s;\n";
out << "real db = y-t;\n";
out << "real dc = z-u;\n";
for (int j = 0; j < nodes.size(); j++) {
const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder();
if (derivOrder[0] == 0 && derivOrder[1] == 0 && derivOrder[2] == 0) {
out << "real value[4] = {0, 0, 0, 0};\n";
for (int k = 3; k >= 0; k--)
for (int m = 0; m < 4; m++) {
int base = k + 4*m;
out << "value[" << m << "] = db*value[" << m << "] + ((c[" << base << "].w*da + c[" << base << "].z)*da + c[" << base << "].y)*da + c[" << base << "].x;\n";
}
out << nodeNames[j] << " = value[0] + dc*(value[1] + dc*(value[2] + dc*value[3]));\n";
}
else if (derivOrder[0] == 1 && derivOrder[1] == 0 && derivOrder[2] == 0) {
out << "real derivx[4] = {0, 0, 0, 0};\n";
for (int k = 3; k >= 0; k--)
for (int m = 0; m < 4; m++) {
int base = k + 4*m;
out << "derivx[" << m << "] = db*derivx[" << m << "] + (3*c[" << base << "].w*da + 2*c[" << base << "].z)*da + c[" << base << "].y;\n";
}
out << nodeNames[j] << " = derivx[0] + dc*(derivx[1] + dc*(derivx[2] + dc*derivx[3]));\n";
out << nodeNames[j] << " *= " << paramsFloat[9] << ";\n";
}
else if (derivOrder[0] == 0 && derivOrder[1] == 1 && derivOrder[2] == 0) {
const string suffixes[] = {".x", ".y", ".z", ".w"};
out << "real derivy[4] = {0, 0, 0, 0};\n";
for (int k = 3; k >= 0; k--)
for (int m = 0; m < 4; m++) {
int base = 4*m;
string suffix = suffixes[m];
out << "derivy[" << m << "] = da*derivy[" << m << "] + (3*c[" << (base+3) << "]" << suffix << "*db + 2*c[" << (base+2) << "]" << suffix << ")*db + c[" << (base+1) << "]" << suffix << ";\n";
}
out << nodeNames[j] << " = derivy[0] + dc*(derivy[1] + dc*(derivy[2] + dc*derivy[3]));\n";
out << nodeNames[j] << " *= " << paramsFloat[10] << ";\n";
}
else if (derivOrder[0] == 0 && derivOrder[1] == 0 && derivOrder[2] == 1) {
out << "real derivz[4] = {0, 0, 0, 0};\n";
for (int k = 3; k >= 0; k--)
for (int m = 0; m < 4; m++) {
int base = k + 4*m;
out << "derivz[" << m << "] = db*derivz[" << m << "] + ((c[" << base << "].w*da + c[" << base << "].z)*da + c[" << base << "].y)*da + c[" << base << "].x;\n";
}
out << nodeNames[j] << " = derivz[1] + dc*(2*derivz[2] + dc*3*derivz[3]);\n";
out << nodeNames[j] << " *= " << paramsFloat[11] << ";\n";
} }
else else
throw OpenMMException("Unsupported derivative order for Continuous2DFunction"); throw OpenMMException("Unsupported derivative order for Continuous2DFunction");
...@@ -474,6 +538,31 @@ vector<float> CudaExpressionUtilities::computeFunctionCoefficients(const Tabulat ...@@ -474,6 +538,31 @@ vector<float> CudaExpressionUtilities::computeFunctionCoefficients(const Tabulat
width = 4; width = 4;
return f; return f;
} }
if (dynamic_cast<const Continuous3DFunction*>(&function) != NULL) {
// Compute the spline coefficients.
const Continuous3DFunction& fn = dynamic_cast<const Continuous3DFunction&>(function);
vector<double> values;
int xsize, ysize, zsize;
double xmin, xmax, ymin, ymax, zmin, zmax;
fn.getFunctionParameters(xsize, ysize, zsize, values, xmin, xmax, ymin, ymax, zmin, zmax);
vector<double> x(xsize), y(ysize), z(zsize);
for (int i = 0; i < xsize; i++)
x[i] = xmin+i*(xmax-xmin)/(xsize-1);
for (int i = 0; i < ysize; i++)
y[i] = ymin+i*(ymax-ymin)/(ysize-1);
for (int i = 0; i < zsize; i++)
z[i] = zmin+i*(zmax-zmin)/(zsize-1);
vector<vector<double> > c;
SplineFitter::create3DNaturalSpline(x, y, z, values, c);
vector<float> f(64*c.size());
for (int i = 0; i < (int) c.size(); i++) {
for (int j = 0; j < 64; j++)
f[64*i+j] = (float) c[i][j];
}
width = 4;
return f;
}
if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL) { if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL) {
// Record the tabulated values. // Record the tabulated values.
...@@ -546,6 +635,25 @@ vector<vector<double> > CudaExpressionUtilities::computeFunctionParameters(const ...@@ -546,6 +635,25 @@ vector<vector<double> > CudaExpressionUtilities::computeFunctionParameters(const
params[i].push_back((xsize-1)/(xmax-xmin)); params[i].push_back((xsize-1)/(xmax-xmin));
params[i].push_back((ysize-1)/(ymax-ymin)); params[i].push_back((ysize-1)/(ymax-ymin));
} }
else if (dynamic_cast<const Continuous3DFunction*>(functions[i]) != NULL) {
const Continuous3DFunction& fn = dynamic_cast<const Continuous3DFunction&>(*functions[i]);
vector<double> values;
int xsize, ysize, zsize;
double xmin, xmax, ymin, ymax, zmin, zmax;
fn.getFunctionParameters(xsize, ysize, zsize, values, xmin, xmax, ymin, ymax, zmin, zmax);
params[i].push_back(xsize-1);
params[i].push_back(ysize-1);
params[i].push_back(zsize-1);
params[i].push_back(xmin);
params[i].push_back(xmax);
params[i].push_back(ymin);
params[i].push_back(ymax);
params[i].push_back(zmin);
params[i].push_back(zmax);
params[i].push_back((xsize-1)/(xmax-xmin));
params[i].push_back((ysize-1)/(ymax-ymin));
params[i].push_back((zsize-1)/(zmax-zmin));
}
else if (dynamic_cast<const Discrete1DFunction*>(functions[i]) != NULL) { else if (dynamic_cast<const Discrete1DFunction*>(functions[i]) != NULL) {
const Discrete1DFunction& fn = dynamic_cast<const Discrete1DFunction&>(*functions[i]); const Discrete1DFunction& fn = dynamic_cast<const Discrete1DFunction&>(*functions[i]);
vector<double> values; vector<double> values;
...@@ -580,6 +688,8 @@ Lepton::CustomFunction* CudaExpressionUtilities::getFunctionPlaceholder(const Ta ...@@ -580,6 +688,8 @@ Lepton::CustomFunction* CudaExpressionUtilities::getFunctionPlaceholder(const Ta
return &fp1; return &fp1;
if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL) if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL)
return &fp2; return &fp2;
if (dynamic_cast<const Continuous3DFunction*>(&function) != NULL)
return &fp3;
if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL) if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL)
return &fp1; return &fp1;
if (dynamic_cast<const Discrete2DFunction*>(&function) != NULL) if (dynamic_cast<const Discrete2DFunction*>(&function) != NULL)
......
platforms/cuda/tests/TestCudaCustomCompoundBondForce.cpp
View file @ 3279c73a
...@@ -199,6 +199,103 @@ void testParallelComputation() { ...@@ -199,6 +199,103 @@ void testParallelComputation() {
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5); ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
} }
void testContinuous2DFunction() {
const int xsize = 10;
const int ysize = 11;
const double xmin = 0.4;
const double xmax = 1.1;
const double ymin = 0.0;
const double ymax = 0.9;
System system;
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomCompoundBondForce* forceField = new CustomCompoundBondForce(1, "fn(x1,y1)+1");
vector<int> particles(1, 0);
forceField->addBond(particles, vector<double>());
vector<double> table(xsize*ysize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
double x = xmin + i*(xmax-xmin)/xsize;
double y = ymin + j*(ymax-ymin)/ysize;
table[i+xsize*j] = sin(0.25*x)*cos(0.33*y);
}
}
forceField->addFunction("fn", new Continuous2DFunction(xsize, ysize, table, xmin, xmax, ymin, ymax));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
positions[0] = Vec3(x, y, 1.5);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double energy = 1;
Vec3 force(0, 0, 0);
if (x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
energy = sin(0.25*x)*cos(0.33*y)+1;
force[0] = -0.25*cos(0.25*x)*cos(0.33*y);
force[1] = 0.3*sin(0.25*x)*sin(0.33*y);
}
ASSERT_EQUAL_VEC(force, forces[0], 0.1);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.05);
}
}
}
void testContinuous3DFunction() {
const int xsize = 10;
const int ysize = 11;
const int zsize = 12;
const double xmin = 0.4;
const double xmax = 1.1;
const double ymin = 0.0;
const double ymax = 0.9;
const double zmin = 0.2;
const double zmax = 1.3;
System system;
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomCompoundBondForce* forceField = new CustomCompoundBondForce(1, "fn(x1,y1,z1)+1");
vector<int> particles(1, 0);
forceField->addBond(particles, vector<double>());
vector<double> table(xsize*ysize*zsize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
for (int k = 0; k < zsize; k++) {
double x = xmin + i*(xmax-xmin)/xsize;
double y = ymin + j*(ymax-ymin)/ysize;
double z = zmin + k*(zmax-zmin)/zsize;
table[i+xsize*j+xsize*ysize*k] = sin(0.25*x)*cos(0.33*y)*(1+z);
}
}
}
forceField->addFunction("fn", new Continuous3DFunction(xsize, ysize, zsize, table, xmin, xmax, ymin, ymax, zmin, zmax));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
for (double z = zmin-0.15; z < zmax+0.2; z += 0.1) {
positions[0] = Vec3(x, y, z);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double energy = 1;
Vec3 force(0, 0, 0);
if (x >= xmin && x <= xmax && y >= ymin && y <= ymax && z >= zmin && z <= zmax) {
energy = sin(0.25*x)*cos(0.33*y)*(1.0+z)+1;
force[0] = -0.25*cos(0.25*x)*cos(0.33*y)*(1.0+z);
force[1] = 0.3*sin(0.25*x)*sin(0.33*y)*(1.0+z);
force[2] = -sin(0.25*x)*cos(0.33*y);
}
ASSERT_EQUAL_VEC(force, forces[0], 0.1);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.05);
}
}
}
}
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
try { try {
if (argc > 1) if (argc > 1)
...@@ -206,6 +303,8 @@ int main(int argc, char* argv[]) { ...@@ -206,6 +303,8 @@ int main(int argc, char* argv[]) {
testBond(); testBond();
testPositionDependence(); testPositionDependence();
testParallelComputation(); testParallelComputation();
testContinuous2DFunction();
testContinuous3DFunction();
} }
catch(const exception& e) { catch(const exception& e) {
cout << "exception: " << e.what() << endl; cout << "exception: " << e.what() << endl;
......
platforms/cuda/tests/TestCudaCustomNonbondedForce.cpp
View file @ 3279c73a
...@@ -347,6 +347,64 @@ void testContinuous2DFunction() { ...@@ -347,6 +347,64 @@ void testContinuous2DFunction() {
} }
} }
void testContinuous3DFunction() {
const int xsize = 10;
const int ysize = 11;
const int zsize = 12;
const double xmin = 0.4;
const double xmax = 1.1;
const double ymin = 0.0;
const double ymax = 0.9;
const double zmin = 0.2;
const double zmax = 1.3;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r,a,b)+1");
forceField->addGlobalParameter("a", 0.0);
forceField->addGlobalParameter("b", 0.0);
forceField->addParticle(vector<double>());
forceField->addParticle(vector<double>());
vector<double> table(xsize*ysize*zsize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
for (int k = 0; k < zsize; k++) {
double x = xmin + i*(xmax-xmin)/xsize;
double y = ymin + j*(ymax-ymin)/ysize;
double z = zmin + k*(zmax-zmin)/zsize;
table[i+xsize*j+xsize*ysize*k] = sin(0.25*x)*cos(0.33*y)*(1+z);
}
}
}
forceField->addFunction("fn", new Continuous3DFunction(xsize, ysize, zsize, table, xmin, xmax, ymin, ymax, zmin, zmax));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
for (double z = zmin-0.15; z < zmax+0.2; z += 0.1) {
positions[1] = Vec3(x, 0, 0);
context.setParameter("a", y);
context.setParameter("b", z);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double energy = 1;
double force = 0;
if (x >= xmin && x <= xmax && y >= ymin && y <= ymax && z >= zmin && z <= zmax) {
energy = sin(0.25*x)*cos(0.33*y)*(1.0+z)+1.0;
force = -0.25*cos(0.25*x)*cos(0.33*y)*(1.0+z);
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], 0.1);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], 0.1);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.05);
}
}
}
}
void testDiscrete1DFunction() { void testDiscrete1DFunction() {
System system; System system;
system.addParticle(1.0); system.addParticle(1.0);
...@@ -869,6 +927,7 @@ int main(int argc, char* argv[]) { ...@@ -869,6 +927,7 @@ int main(int argc, char* argv[]) {
testPeriodic(); testPeriodic();
testContinuous1DFunction(); testContinuous1DFunction();
testContinuous2DFunction(); testContinuous2DFunction();
testContinuous3DFunction();
testDiscrete1DFunction(); testDiscrete1DFunction();
testDiscrete2DFunction(); testDiscrete2DFunction();
testDiscrete3DFunction(); testDiscrete3DFunction();
......
platforms/opencl/src/OpenCLExpressionUtilities.cpp
View file @ 3279c73a
...@@ -133,36 +133,100 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre ...@@ -133,36 +133,100 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre
out << "x = (x-" << paramsFloat[2] << ")*" << paramsFloat[6] << ";\n"; out << "x = (x-" << paramsFloat[2] << ")*" << paramsFloat[6] << ";\n";
out << "y = (y-" << paramsFloat[4] << ")*" << paramsFloat[7] << ";\n"; out << "y = (y-" << paramsFloat[4] << ")*" << paramsFloat[7] << ";\n";
out << "int s = min((int) floor(x), " << paramsInt[0] << ");\n"; out << "int s = min((int) floor(x), " << paramsInt[0] << ");\n";
out << "int t = min((int) floor(y), " << paramsInt[0] << ");\n"; out << "int t = min((int) floor(y), " << paramsInt[1] << ");\n";
out << "int coeffIndex = 4*(s+" << paramsInt[0] << "*t);\n"; out << "int coeffIndex = 4*(s+" << paramsInt[0] << "*t);\n";
out << "float4 c[4];\n"; out << "float4 c[4];\n";
out << "c[0] = " << functionNames[i].second << "[coeffIndex];\n"; for (int j = 0; j < 4; j++)
out << "c[1] = " << functionNames[i].second << "[coeffIndex+1];\n"; out << "c[" << j << "] = " << functionNames[i].second << "[coeffIndex+" << j << "];\n";
out << "c[2] = " << functionNames[i].second << "[coeffIndex+2];\n"; out << "real da = x-s;\n";
out << "c[3] = " << functionNames[i].second << "[coeffIndex+3];\n"; out << "real db = y-t;\n";
out << "real da = x-s;";
out << "real db = y-t;";
for (int j = 0; j < nodes.size(); j++) { for (int j = 0; j < nodes.size(); j++) {
const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder(); const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder();
if (derivOrder[0] == 0 && derivOrder[1] == 0) { if (derivOrder[0] == 0 && derivOrder[1] == 0) {
out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[3].w*db + c[3].z)*db + c[3].y)*db + c[3].x;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[3].w*db + c[3].z)*db + c[3].y)*db + c[3].x;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[2].w*db + c[2].z)*db + c[2].y)*db + c[2].x;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[2].w*db + c[2].z)*db + c[2].y)*db + c[2].x;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[1].w*db + c[1].z)*db + c[1].y)*db + c[1].x;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[1].w*db + c[1].z)*db + c[1].y)*db + c[1].x;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[0].w*db + c[0].z)*db + c[0].y)*db + c[0].x;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + ((c[0].w*db + c[0].z)*db + c[0].y)*db + c[0].x;\n";
} }
else if (derivOrder[0] == 1 && derivOrder[1] == 0) { else if (derivOrder[0] == 1 && derivOrder[1] == 0) {
out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].w*da + 2.0f*c[2].w)*da + c[1].w;"; out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].w*da + 2.0f*c[2].w)*da + c[1].w;\n";
out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].z*da + 2.0f*c[2].z)*da + c[1].z;"; out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].z*da + 2.0f*c[2].z)*da + c[1].z;\n";
out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].y*da + 2.0f*c[2].y)*da + c[1].y;"; out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].y*da + 2.0f*c[2].y)*da + c[1].y;\n";
out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].x*da + 2.0f*c[2].x)*da + c[1].x;"; out << nodeNames[j] << " = db*" << nodeNames[j] << " + (3.0f*c[3].x*da + 2.0f*c[2].x)*da + c[1].x;\n";
out << nodeNames[j] << " *= " << paramsFloat[6] << ";"; out << nodeNames[j] << " *= " << paramsFloat[6] << ";\n";
} }
else if (derivOrder[0] == 0 && derivOrder[1] == 1) { else if (derivOrder[0] == 0 && derivOrder[1] == 1) {
out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[3].w*db + 2.0f*c[3].z)*db + c[3].y;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[3].w*db + 2.0f*c[3].z)*db + c[3].y;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[2].w*db + 2.0f*c[2].z)*db + c[2].y;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[2].w*db + 2.0f*c[2].z)*db + c[2].y;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[1].w*db + 2.0f*c[1].z)*db + c[1].y;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[1].w*db + 2.0f*c[1].z)*db + c[1].y;\n";
out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[0].w*db + 2.0f*c[0].z)*db + c[0].y;"; out << nodeNames[j] << " = da*" << nodeNames[j] << " + (3.0f*c[0].w*db + 2.0f*c[0].z)*db + c[0].y;\n";
out << nodeNames[j] << " *= " << paramsFloat[7] << ";"; out << nodeNames[j] << " *= " << paramsFloat[7] << ";\n";
}
else
throw OpenMMException("Unsupported derivative order for Continuous2DFunction");
}
out << "}\n";
}
else if (dynamic_cast<const Continuous3DFunction*>(functions[i]) != NULL) {
out << "real x = " << getTempName(node.getChildren()[0], temps) << ";\n";
out << "real y = " << getTempName(node.getChildren()[1], temps) << ";\n";
out << "real z = " << getTempName(node.getChildren()[2], temps) << ";\n";
out << "if (x >= " << paramsFloat[3] << " && x <= " << paramsFloat[4] << " && y >= " << paramsFloat[5] << " && y <= " << paramsFloat[6] << " && z >= " << paramsFloat[7] << " && z <= " << paramsFloat[8] << ") {\n";
out << "x = (x-" << paramsFloat[3] << ")*" << paramsFloat[9] << ";\n";
out << "y = (y-" << paramsFloat[5] << ")*" << paramsFloat[10] << ";\n";
out << "z = (z-" << paramsFloat[7] << ")*" << paramsFloat[11] << ";\n";
out << "int s = min((int) floor(x), " << paramsInt[0] << ");\n";
out << "int t = min((int) floor(y), " << paramsInt[1] << ");\n";
out << "int u = min((int) floor(z), " << paramsInt[2] << ");\n";
out << "int coeffIndex = 16*(s+" << paramsInt[0] << "*(t+" << paramsInt[1] << "*u));\n";
out << "float4 c[16];\n";
for (int j = 0; j < 16; j++)
out << "c[" << j << "] = " << functionNames[i].second << "[coeffIndex+" << j << "];\n";
out << "real da = x-s;\n";
out << "real db = y-t;\n";
out << "real dc = z-u;\n";
for (int j = 0; j < nodes.size(); j++) {
const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder();
if (derivOrder[0] == 0 && derivOrder[1] == 0 && derivOrder[2] == 0) {
out << "real value[4] = {0, 0, 0, 0};\n";
for (int k = 3; k >= 0; k--)
for (int m = 0; m < 4; m++) {
int base = k + 4*m;
out << "value[" << m << "] = db*value[" << m << "] + ((c[" << base << "].w*da + c[" << base << "].z)*da + c[" << base << "].y)*da + c[" << base << "].x;\n";
}
out << nodeNames[j] << " = value[0] + dc*(value[1] + dc*(value[2] + dc*value[3]));\n";
}
else if (derivOrder[0] == 1 && derivOrder[1] == 0 && derivOrder[2] == 0) {
out << "real derivx[4] = {0, 0, 0, 0};\n";
for (int k = 3; k >= 0; k--)
for (int m = 0; m < 4; m++) {
int base = k + 4*m;
out << "derivx[" << m << "] = db*derivx[" << m << "] + (3*c[" << base << "].w*da + 2*c[" << base << "].z)*da + c[" << base << "].y;\n";
}
out << nodeNames[j] << " = derivx[0] + dc*(derivx[1] + dc*(derivx[2] + dc*derivx[3]));\n";
out << nodeNames[j] << " *= " << paramsFloat[9] << ";\n";
}
else if (derivOrder[0] == 0 && derivOrder[1] == 1 && derivOrder[2] == 0) {
const string suffixes[] = {".x", ".y", ".z", ".w"};
out << "real derivy[4] = {0, 0, 0, 0};\n";
for (int k = 3; k >= 0; k--)
for (int m = 0; m < 4; m++) {
int base = 4*m;
string suffix = suffixes[m];
out << "derivy[" << m << "] = da*derivy[" << m << "] + (3*c[" << (base+3) << "]" << suffix << "*db + 2*c[" << (base+2) << "]" << suffix << ")*db + c[" << (base+1) << "]" << suffix << ";\n";
}
out << nodeNames[j] << " = derivy[0] + dc*(derivy[1] + dc*(derivy[2] + dc*derivy[3]));\n";
out << nodeNames[j] << " *= " << paramsFloat[10] << ";\n";
}
else if (derivOrder[0] == 0 && derivOrder[1] == 0 && derivOrder[2] == 1) {
out << "real derivz[4] = {0, 0, 0, 0};\n";
for (int k = 3; k >= 0; k--)
for (int m = 0; m < 4; m++) {
int base = k + 4*m;
out << "derivz[" << m << "] = db*derivz[" << m << "] + ((c[" << base << "].w*da + c[" << base << "].z)*da + c[" << base << "].y)*da + c[" << base << "].x;\n";
}
out << nodeNames[j] << " = derivz[1] + dc*(2*derivz[2] + dc*3*derivz[3]);\n";
out << nodeNames[j] << " *= " << paramsFloat[11] << ";\n";
} }
else else
throw OpenMMException("Unsupported derivative order for Continuous2DFunction"); throw OpenMMException("Unsupported derivative order for Continuous2DFunction");
...@@ -474,6 +538,31 @@ vector<float> OpenCLExpressionUtilities::computeFunctionCoefficients(const Tabul ...@@ -474,6 +538,31 @@ vector<float> OpenCLExpressionUtilities::computeFunctionCoefficients(const Tabul
width = 4; width = 4;
return f; return f;
} }
if (dynamic_cast<const Continuous3DFunction*>(&function) != NULL) {
// Compute the spline coefficients.
const Continuous3DFunction& fn = dynamic_cast<const Continuous3DFunction&>(function);
vector<double> values;
int xsize, ysize, zsize;
double xmin, xmax, ymin, ymax, zmin, zmax;
fn.getFunctionParameters(xsize, ysize, zsize, values, xmin, xmax, ymin, ymax, zmin, zmax);
vector<double> x(xsize), y(ysize), z(zsize);
for (int i = 0; i < xsize; i++)
x[i] = xmin+i*(xmax-xmin)/(xsize-1);
for (int i = 0; i < ysize; i++)
y[i] = ymin+i*(ymax-ymin)/(ysize-1);
for (int i = 0; i < zsize; i++)
z[i] = zmin+i*(zmax-zmin)/(zsize-1);
vector<vector<double> > c;
SplineFitter::create3DNaturalSpline(x, y, z, values, c);
vector<float> f(64*c.size());
for (int i = 0; i < (int) c.size(); i++) {
for (int j = 0; j < 64; j++)
f[64*i+j] = (float) c[i][j];
}
width = 4;
return f;
}
if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL) { if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL) {
// Record the tabulated values. // Record the tabulated values.
...@@ -546,6 +635,25 @@ vector<vector<double> > OpenCLExpressionUtilities::computeFunctionParameters(con ...@@ -546,6 +635,25 @@ vector<vector<double> > OpenCLExpressionUtilities::computeFunctionParameters(con
params[i].push_back((xsize-1)/(xmax-xmin)); params[i].push_back((xsize-1)/(xmax-xmin));
params[i].push_back((ysize-1)/(ymax-ymin)); params[i].push_back((ysize-1)/(ymax-ymin));
} }
else if (dynamic_cast<const Continuous3DFunction*>(functions[i]) != NULL) {
const Continuous3DFunction& fn = dynamic_cast<const Continuous3DFunction&>(*functions[i]);
vector<double> values;
int xsize, ysize, zsize;
double xmin, xmax, ymin, ymax, zmin, zmax;
fn.getFunctionParameters(xsize, ysize, zsize, values, xmin, xmax, ymin, ymax, zmin, zmax);
params[i].push_back(xsize-1);
params[i].push_back(ysize-1);
params[i].push_back(zsize-1);
params[i].push_back(xmin);
params[i].push_back(xmax);
params[i].push_back(ymin);
params[i].push_back(ymax);
params[i].push_back(zmin);
params[i].push_back(zmax);
params[i].push_back((xsize-1)/(xmax-xmin));
params[i].push_back((ysize-1)/(ymax-ymin));
params[i].push_back((zsize-1)/(zmax-zmin));
}
else if (dynamic_cast<const Discrete1DFunction*>(functions[i]) != NULL) { else if (dynamic_cast<const Discrete1DFunction*>(functions[i]) != NULL) {
const Discrete1DFunction& fn = dynamic_cast<const Discrete1DFunction&>(*functions[i]); const Discrete1DFunction& fn = dynamic_cast<const Discrete1DFunction&>(*functions[i]);
vector<double> values; vector<double> values;
...@@ -580,6 +688,8 @@ Lepton::CustomFunction* OpenCLExpressionUtilities::getFunctionPlaceholder(const ...@@ -580,6 +688,8 @@ Lepton::CustomFunction* OpenCLExpressionUtilities::getFunctionPlaceholder(const
return &fp1; return &fp1;
if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL) if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL)
return &fp2; return &fp2;
if (dynamic_cast<const Continuous3DFunction*>(&function) != NULL)
return &fp3;
if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL) if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL)
return &fp1; return &fp1;
if (dynamic_cast<const Discrete2DFunction*>(&function) != NULL) if (dynamic_cast<const Discrete2DFunction*>(&function) != NULL)
......
platforms/opencl/tests/TestOpenCLCustomCompoundBondForce.cpp
View file @ 3279c73a
...@@ -199,6 +199,103 @@ void testParallelComputation() { ...@@ -199,6 +199,103 @@ void testParallelComputation() {
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5); ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
} }
void testContinuous2DFunction() {
const int xsize = 10;
const int ysize = 11;
const double xmin = 0.4;
const double xmax = 1.1;
const double ymin = 0.0;
const double ymax = 0.9;
System system;
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomCompoundBondForce* forceField = new CustomCompoundBondForce(1, "fn(x1,y1)+1");
vector<int> particles(1, 0);
forceField->addBond(particles, vector<double>());
vector<double> table(xsize*ysize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
double x = xmin + i*(xmax-xmin)/xsize;
double y = ymin + j*(ymax-ymin)/ysize;
table[i+xsize*j] = sin(0.25*x)*cos(0.33*y);
}
}
forceField->addFunction("fn", new Continuous2DFunction(xsize, ysize, table, xmin, xmax, ymin, ymax));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
positions[0] = Vec3(x, y, 1.5);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double energy = 1;
Vec3 force(0, 0, 0);
if (x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
energy = sin(0.25*x)*cos(0.33*y)+1;
force[0] = -0.25*cos(0.25*x)*cos(0.33*y);
force[1] = 0.3*sin(0.25*x)*sin(0.33*y);
}
ASSERT_EQUAL_VEC(force, forces[0], 0.1);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.05);
}
}
}
void testContinuous3DFunction() {
const int xsize = 10;
const int ysize = 11;
const int zsize = 12;
const double xmin = 0.4;
const double xmax = 1.1;
const double ymin = 0.0;
const double ymax = 0.9;
const double zmin = 0.2;
const double zmax = 1.3;
System system;
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomCompoundBondForce* forceField = new CustomCompoundBondForce(1, "fn(x1,y1,z1)+1");
vector<int> particles(1, 0);
forceField->addBond(particles, vector<double>());
vector<double> table(xsize*ysize*zsize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
for (int k = 0; k < zsize; k++) {
double x = xmin + i*(xmax-xmin)/xsize;
double y = ymin + j*(ymax-ymin)/ysize;
double z = zmin + k*(zmax-zmin)/zsize;
table[i+xsize*j+xsize*ysize*k] = sin(0.25*x)*cos(0.33*y)*(1+z);
}
}
}
forceField->addFunction("fn", new Continuous3DFunction(xsize, ysize, zsize, table, xmin, xmax, ymin, ymax, zmin, zmax));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
for (double z = zmin-0.15; z < zmax+0.2; z += 0.1) {
positions[0] = Vec3(x, y, z);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double energy = 1;
Vec3 force(0, 0, 0);
if (x >= xmin && x <= xmax && y >= ymin && y <= ymax && z >= zmin && z <= zmax) {
energy = sin(0.25*x)*cos(0.33*y)*(1.0+z)+1;
force[0] = -0.25*cos(0.25*x)*cos(0.33*y)*(1.0+z);
force[1] = 0.3*sin(0.25*x)*sin(0.33*y)*(1.0+z);
force[2] = -sin(0.25*x)*cos(0.33*y);
}
ASSERT_EQUAL_VEC(force, forces[0], 0.1);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.05);
}
}
}
}
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
try { try {
if (argc > 1) if (argc > 1)
...@@ -206,6 +303,8 @@ int main(int argc, char* argv[]) { ...@@ -206,6 +303,8 @@ int main(int argc, char* argv[]) {
testBond(); testBond();
testPositionDependence(); testPositionDependence();
testParallelComputation(); testParallelComputation();
testContinuous2DFunction();
testContinuous3DFunction();
} }
catch(const exception& e) { catch(const exception& e) {
cout << "exception: " << e.what() << endl; cout << "exception: " << e.what() << endl;
......
platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp
View file @ 3279c73a
...@@ -347,6 +347,64 @@ void testContinuous2DFunction() { ...@@ -347,6 +347,64 @@ void testContinuous2DFunction() {
} }
} }
void testContinuous3DFunction() {
const int xsize = 10;
const int ysize = 11;
const int zsize = 12;
const double xmin = 0.4;
const double xmax = 1.1;
const double ymin = 0.0;
const double ymax = 0.9;
const double zmin = 0.2;
const double zmax = 1.3;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r,a,b)+1");
forceField->addGlobalParameter("a", 0.0);
forceField->addGlobalParameter("b", 0.0);
forceField->addParticle(vector<double>());
forceField->addParticle(vector<double>());
vector<double> table(xsize*ysize*zsize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
for (int k = 0; k < zsize; k++) {
double x = xmin + i*(xmax-xmin)/xsize;
double y = ymin + j*(ymax-ymin)/ysize;
double z = zmin + k*(zmax-zmin)/zsize;
table[i+xsize*j+xsize*ysize*k] = sin(0.25*x)*cos(0.33*y)*(1+z);
}
}
}
forceField->addFunction("fn", new Continuous3DFunction(xsize, ysize, zsize, table, xmin, xmax, ymin, ymax, zmin, zmax));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
for (double z = zmin-0.15; z < zmax+0.2; z += 0.1) {
positions[1] = Vec3(x, 0, 0);
context.setParameter("a", y);
context.setParameter("b", z);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double energy = 1;
double force = 0;
if (x >= xmin && x <= xmax && y >= ymin && y <= ymax && z >= zmin && z <= zmax) {
energy = sin(0.25*x)*cos(0.33*y)*(1.0+z)+1.0;
force = -0.25*cos(0.25*x)*cos(0.33*y)*(1.0+z);
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], 0.1);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], 0.1);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.05);
}
}
}
}
void testDiscrete1DFunction() { void testDiscrete1DFunction() {
System system; System system;
system.addParticle(1.0); system.addParticle(1.0);
...@@ -868,6 +926,7 @@ int main(int argc, char* argv[]) { ...@@ -868,6 +926,7 @@ int main(int argc, char* argv[]) {
testPeriodic(); testPeriodic();
testContinuous1DFunction(); testContinuous1DFunction();
testContinuous2DFunction(); testContinuous2DFunction();
testContinuous3DFunction();
testDiscrete1DFunction(); testDiscrete1DFunction();
testDiscrete2DFunction(); testDiscrete2DFunction();
testDiscrete3DFunction(); testDiscrete3DFunction();
......
platforms/reference/tests/TestReferenceCustomCompoundBondForce.cpp
View file @ 3279c73a
...@@ -167,10 +167,111 @@ void testPositionDependence() { ...@@ -167,10 +167,111 @@ void testPositionDependence() {
ASSERT_EQUAL_VEC(Vec3(-0.3, -2, 0), state.getForces()[1], 1e-5); ASSERT_EQUAL_VEC(Vec3(-0.3, -2, 0), state.getForces()[1], 1e-5);
} }
void testContinuous2DFunction() {
const int xsize = 10;
const int ysize = 11;
const double xmin = 0.4;
const double xmax = 1.1;
const double ymin = 0.0;
const double ymax = 0.9;
ReferencePlatform platform;
System system;
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomCompoundBondForce* forceField = new CustomCompoundBondForce(1, "fn(x1,y1)+1");
vector<int> particles(1, 0);
forceField->addBond(particles, vector<double>());
vector<double> table(xsize*ysize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
double x = xmin + i*(xmax-xmin)/xsize;
double y = ymin + j*(ymax-ymin)/ysize;
table[i+xsize*j] = sin(0.25*x)*cos(0.33*y);
}
}
forceField->addFunction("fn", new Continuous2DFunction(xsize, ysize, table, xmin, xmax, ymin, ymax));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
positions[0] = Vec3(x, y, 1.5);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double energy = 1;
Vec3 force(0, 0, 0);
if (x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
energy = sin(0.25*x)*cos(0.33*y)+1;
force[0] = -0.25*cos(0.25*x)*cos(0.33*y);
force[1] = 0.3*sin(0.25*x)*sin(0.33*y);
}
ASSERT_EQUAL_VEC(force, forces[0], 0.1);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.05);
}
}
}
void testContinuous3DFunction() {
const int xsize = 10;
const int ysize = 11;
const int zsize = 12;
const double xmin = 0.4;
const double xmax = 1.1;
const double ymin = 0.0;
const double ymax = 0.9;
const double zmin = 0.2;
const double zmax = 1.3;
ReferencePlatform platform;
System system;
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomCompoundBondForce* forceField = new CustomCompoundBondForce(1, "fn(x1,y1,z1)+1");
vector<int> particles(1, 0);
forceField->addBond(particles, vector<double>());
vector<double> table(xsize*ysize*zsize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
for (int k = 0; k < zsize; k++) {
double x = xmin + i*(xmax-xmin)/xsize;
double y = ymin + j*(ymax-ymin)/ysize;
double z = zmin + k*(zmax-zmin)/zsize;
table[i+xsize*j+xsize*ysize*k] = sin(0.25*x)*cos(0.33*y)*(1+z);
}
}
}
forceField->addFunction("fn", new Continuous3DFunction(xsize, ysize, zsize, table, xmin, xmax, ymin, ymax, zmin, zmax));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
for (double z = zmin-0.15; z < zmax+0.2; z += 0.1) {
positions[0] = Vec3(x, y, z);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double energy = 1;
Vec3 force(0, 0, 0);
if (x >= xmin && x <= xmax && y >= ymin && y <= ymax && z >= zmin && z <= zmax) {
energy = sin(0.25*x)*cos(0.33*y)*(1.0+z)+1;
force[0] = -0.25*cos(0.25*x)*cos(0.33*y)*(1.0+z);
force[1] = 0.3*sin(0.25*x)*sin(0.33*y)*(1.0+z);
force[2] = -sin(0.25*x)*cos(0.33*y);
}
ASSERT_EQUAL_VEC(force, forces[0], 0.1);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.05);
}
}
}
}
int main() { int main() {
try { try {
testBond(); testBond();
testPositionDependence(); testPositionDependence();
testContinuous2DFunction();
testContinuous3DFunction();
} }
catch(const exception& e) { catch(const exception& e) {
cout << "exception: " << e.what() << endl; cout << "exception: " << e.what() << endl;
......
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