Skip to content

GitLab

Commit 708c4246 authored by Lee-Ping's avatar Lee-Ping
Browse files

Merge branch 'master' of github.com:SimTk/openmm

parents d284e2b8 ec90b4bb
Hide whitespace changes
Inline Side-by-side
Showing with 1124 additions and 98 deletions
platforms/opencl/src/OpenCLPlatform.cpp
View file @ 708c4246
......@@ -40,9 +40,15 @@ using std::string;
using std::stringstream;
using std::vector;
#ifdef OPENMM_OPENCL_BUILDING_STATIC_LIBRARY
extern "C" void registerOpenCLPlatform() {
Platform::registerPlatform(new OpenCLPlatform());
}
#else
extern "C" OPENMM_EXPORT_OPENCL void registerPlatforms() {
Platform::registerPlatform(new OpenCLPlatform());
}
#endif
OpenCLPlatform::OpenCLPlatform() {
OpenCLKernelFactory* factory = new OpenCLKernelFactory();
......
platforms/opencl/staticTarget/CMakeLists.txt 0 → 100644
View file @ 708c4246
#
# Include OpenCL related files.
#
# INCLUDE(${CMAKE_CURRENT_SOURCE_DIR}/../FindOpenCL.cmake)
INCLUDE_DIRECTORIES(${OPENCL_INCLUDE_DIR})
FILE(GLOB OPENCL_KERNELS ${CL_SOURCE_DIR}/kernels/*.cl)
ADD_CUSTOM_COMMAND(OUTPUT ${CL_KERNELS_CPP} ${CL_KERNELS_H}
COMMAND ${CMAKE_COMMAND}
ARGS -D CL_SOURCE_DIR=${CL_SOURCE_DIR} -D CL_KERNELS_CPP=${CL_KERNELS_CPP} -D CL_KERNELS_H=${CL_KERNELS_H} -D CL_SOURCE_CLASS=${CL_SOURCE_CLASS} -P ${CMAKE_CURRENT_SOURCE_DIR}/../EncodeCLFiles.cmake
DEPENDS ${OPENCL_KERNELS}
)
SET_SOURCE_FILES_PROPERTIES(${CL_KERNELS_CPP} ${CL_KERNELS_H} PROPERTIES GENERATED TRUE)
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} ${OPENCL_LIBRARIES} ${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_OPENCL_BUILDING_STATIC_LIBRARY -DPTW32_STATIC_LIB")
INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${STATIC_TARGET})
platforms/opencl/tests/TestOpenCLCustomCompoundBondForce.cpp
View file @ 708c4246
......@@ -199,6 +199,103 @@ void testParallelComputation() {
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->addTabulatedFunction("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->addTabulatedFunction("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[]) {
try {
if (argc > 1)
......@@ -206,6 +303,8 @@ int main(int argc, char* argv[]) {
testBond();
testPositionDependence();
testParallelComputation();
testContinuous2DFunction();
testContinuous3DFunction();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/opencl/tests/TestOpenCLCustomGBForce.cpp
View file @ 708c4246
......@@ -277,7 +277,7 @@ void testTabulatedFunction() {
vector<double> table;
for (int i = 0; i < 21; i++)
table.push_back(std::sin(0.25*i));
force->addFunction("fn", table, 1.0, 6.0);
force->addTabulatedFunction("fn", new Continuous1DFunction(table, 1.0, 6.0));
system.addForce(force);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
......
platforms/opencl/tests/TestOpenCLCustomHbondForce.cpp
View file @ 708c4246
......@@ -214,7 +214,7 @@ void testCustomFunctions() {
vector<double> function(2);
function[0] = 0;
function[1] = 1;
custom->addFunction("foo", function, 0, 10);
custom->addTabulatedFunction("foo", new Continuous1DFunction(function, 0, 10));
system.addForce(custom);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
......
platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp
View file @ 708c4246
......@@ -7,7 +7,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: *
* *
......@@ -261,7 +261,7 @@ void testPeriodic() {
ASSERT_EQUAL_TOL(1.9+1+0.9, state.getPotentialEnergy(), TOL);
}
void testTabulatedFunction() {
void testContinuous1DFunction() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
......@@ -271,21 +271,20 @@ void testTabulatedFunction() {
forceField->addParticle(vector<double>());
vector<double> table;
for (int i = 0; i < 21; i++)
table.push_back(std::sin(0.25*i));
forceField->addFunction("fn", table, 1.0, 6.0);
table.push_back(sin(0.25*i));
forceField->addTabulatedFunction("fn", new Continuous1DFunction(table, 1.0, 6.0));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
double tol = 0.01;
for (int i = 1; i < 30; i++) {
double x = (7.0/30.0)*i;
positions[1] = Vec3(x, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double force = (x < 1.0 || x > 6.0 ? 0.0 : -std::cos(x-1.0));
double energy = (x < 1.0 || x > 6.0 ? 0.0 : std::sin(x-1.0))+1.0;
double force = (x < 1.0 || x > 6.0 ? 0.0 : -cos(x-1.0));
double energy = (x < 1.0 || x > 6.0 ? 0.0 : sin(x-1.0))+1.0;
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.02);
......@@ -295,11 +294,211 @@ void testTabulatedFunction() {
positions[1] = Vec3(x, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double energy = (x < 1.0 || x > 6.0 ? 0.0 : std::sin(x-1.0))+1.0;
double energy = (x < 1.0 || x > 6.0 ? 0.0 : sin(x-1.0))+1.0;
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 1e-4);
}
}
void testContinuous2DFunction() {
const int xsize = 20;
const int ysize = 21;
const double xmin = 0.4;
const double xmax = 1.5;
const double ymin = 0.0;
const double ymax = 2.1;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r,a)+1");
forceField->addGlobalParameter("a", 0.0);
forceField->addParticle(vector<double>());
forceField->addParticle(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->addTabulatedFunction("fn", new Continuous2DFunction(xsize, ysize, table, xmin, xmax, ymin, ymax));
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) {
positions[1] = Vec3(x, 0, 0);
context.setParameter("a", y);
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) {
energy = sin(0.25*x)*cos(0.33*y)+1.0;
force = -0.25*cos(0.25*x)*cos(0.33*y);
}
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.02);
}
}
}
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->addTabulatedFunction("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() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r-1)+1");
forceField->addParticle(vector<double>());
forceField->addParticle(vector<double>());
vector<double> table;
for (int i = 0; i < 21; i++)
table.push_back(sin(0.25*i));
forceField->addTabulatedFunction("fn", new Discrete1DFunction(table));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
for (int i = 0; i < (int) table.size(); i++) {
positions[1] = Vec3(i+1, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[0], 1e-6);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], 1e-6);
ASSERT_EQUAL_TOL(table[i]+1.0, state.getPotentialEnergy(), 1e-6);
}
}
void testDiscrete2DFunction() {
const int xsize = 10;
const int ysize = 5;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r-1,a)+1");
forceField->addGlobalParameter("a", 0.0);
forceField->addParticle(vector<double>());
forceField->addParticle(vector<double>());
vector<double> table;
for (int i = 0; i < xsize; i++)
for (int j = 0; j < ysize; j++)
table.push_back(sin(0.25*i)+cos(0.33*j));
forceField->addTabulatedFunction("fn", new Discrete2DFunction(xsize, ysize, table));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
for (int i = 0; i < (int) table.size(); i++) {
positions[1] = Vec3((i%xsize)+1, 0, 0);
context.setPositions(positions);
context.setParameter("a", i/xsize);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[0], 1e-6);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], 1e-6);
ASSERT_EQUAL_TOL(table[i]+1.0, state.getPotentialEnergy(), 1e-6);
}
}
void testDiscrete3DFunction() {
const int xsize = 8;
const int ysize = 5;
const int zsize = 6;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r-1,a,b)+1");
forceField->addGlobalParameter("a", 0.0);
forceField->addGlobalParameter("b", 0.0);
forceField->addParticle(vector<double>());
forceField->addParticle(vector<double>());
vector<double> table;
for (int i = 0; i < xsize; i++)
for (int j = 0; j < ysize; j++)
for (int k = 0; k < zsize; k++)
table.push_back(sin(0.25*i)+cos(0.33*j)+0.12345*k);
forceField->addTabulatedFunction("fn", new Discrete3DFunction(xsize, ysize, zsize, table));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
for (int i = 0; i < (int) table.size(); i++) {
positions[1] = Vec3((i%xsize)+1, 0, 0);
context.setPositions(positions);
context.setParameter("a", (i/xsize)%ysize);
context.setParameter("b", i/(xsize*ysize));
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[0], 1e-6);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], 1e-6);
ASSERT_EQUAL_TOL(table[i]+1.0, state.getPotentialEnergy(), 1e-6);
}
}
void testCoulombLennardJones() {
const int numMolecules = 300;
const int numParticles = numMolecules*2;
......@@ -725,7 +924,12 @@ int main(int argc, char* argv[]) {
testExclusions();
testCutoff();
testPeriodic();
testTabulatedFunction();
testContinuous1DFunction();
testContinuous2DFunction();
testContinuous3DFunction();
testDiscrete1DFunction();
testDiscrete2DFunction();
testDiscrete3DFunction();
testCoulombLennardJones();
testParallelComputation();
testSwitchingFunction();
......
platforms/reference/include/ReferenceProperDihedralBond.h
View file @ 708c4246
......@@ -29,7 +29,7 @@
// ---------------------------------------------------------------------------------------
class ReferenceProperDihedralBond : public ReferenceBondIxn {
class OPENMM_EXPORT ReferenceProperDihedralBond : public ReferenceBondIxn {
private:
......
platforms/reference/include/ReferenceRbDihedralBond.h
View file @ 708c4246
......@@ -29,7 +29,7 @@
// ---------------------------------------------------------------------------------------
class ReferenceRbDihedralBond : public ReferenceBondIxn {
class OPENMM_EXPORT ReferenceRbDihedralBond : public ReferenceBondIxn {
private:
......
platforms/reference/include/ReferenceTabulatedFunction.h 0 → 100644
View file @ 708c4246
#ifndef OPENMM_REFERENCETABULATEDFUNCTION_H_
#define OPENMM_REFERENCETABULATEDFUNCTION_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 "openmm/TabulatedFunction.h"
#include "openmm/internal/windowsExport.h"
#include "lepton/CustomFunction.h"
#include <vector>
namespace OpenMM {
/**
* Given a TabulatedFunction, wrap it in an appropriate subclass of Lepton::CustomFunction.
*/
extern "C" Lepton::CustomFunction* createReferenceTabulatedFunction(const TabulatedFunction& function);
/**
* This class adapts a Continuous1DFunction into a Lepton::CustomFunction.
*/
class OPENMM_EXPORT ReferenceContinuous1DFunction : public Lepton::CustomFunction {
public:
ReferenceContinuous1DFunction(const Continuous1DFunction& function);
int getNumArguments() const;
double evaluate(const double* arguments) const;
double evaluateDerivative(const double* arguments, const int* derivOrder) const;
CustomFunction* clone() const;
private:
const Continuous1DFunction& function;
double min, max;
std::vector<double> x, values, derivs;
};
/**
* This class adapts a Continuous2DFunction into a Lepton::CustomFunction.
*/
class OPENMM_EXPORT ReferenceContinuous2DFunction : public Lepton::CustomFunction {
public:
ReferenceContinuous2DFunction(const Continuous2DFunction& function);
int getNumArguments() const;
double evaluate(const double* arguments) const;
double evaluateDerivative(const double* arguments, const int* derivOrder) const;
CustomFunction* clone() const;
private:
const Continuous2DFunction& function;
int xsize, ysize;
double xmin, xmax, ymin, ymax;
std::vector<double> x, y, values;
std::vector<std::vector<double> > c;
};
/**
* This class adapts a Continuous3DFunction into a Lepton::CustomFunction.
*/
class OPENMM_EXPORT ReferenceContinuous3DFunction : public Lepton::CustomFunction {
public:
ReferenceContinuous3DFunction(const Continuous3DFunction& function);
int getNumArguments() const;
double evaluate(const double* arguments) const;
double evaluateDerivative(const double* arguments, const int* derivOrder) const;
CustomFunction* clone() const;
private:
const Continuous3DFunction& function;
int xsize, ysize, zsize;
double xmin, xmax, ymin, ymax, zmin, zmax;
std::vector<double> x, y, z, values;
std::vector<std::vector<double> > c;
};
/**
* This class adapts a Discrete1DFunction into a Lepton::CustomFunction.
*/
class OPENMM_EXPORT ReferenceDiscrete1DFunction : public Lepton::CustomFunction {
public:
ReferenceDiscrete1DFunction(const Discrete1DFunction& function);
int getNumArguments() const;
double evaluate(const double* arguments) const;
double evaluateDerivative(const double* arguments, const int* derivOrder) const;
CustomFunction* clone() const;
private:
const Discrete1DFunction& function;
std::vector<double> values;
};
/**
* This class adapts a Discrete2DFunction into a Lepton::CustomFunction.
*/
class OPENMM_EXPORT ReferenceDiscrete2DFunction : public Lepton::CustomFunction {
public:
ReferenceDiscrete2DFunction(const Discrete2DFunction& function);
int getNumArguments() const;
double evaluate(const double* arguments) const;
double evaluateDerivative(const double* arguments, const int* derivOrder) const;
CustomFunction* clone() const;
private:
const Discrete2DFunction& function;
int xsize, ysize;
std::vector<double> values;
};
/**
* This class adapts a Discrete3DFunction into a Lepton::CustomFunction.
*/
class OPENMM_EXPORT ReferenceDiscrete3DFunction : public Lepton::CustomFunction {
public:
ReferenceDiscrete3DFunction(const Discrete3DFunction& function);
int getNumArguments() const;
double evaluate(const double* arguments) const;
double evaluateDerivative(const double* arguments, const int* derivOrder) const;
CustomFunction* clone() const;
private:
const Discrete3DFunction& function;
int xsize, ysize, zsize;
std::vector<double> values;
};
} // namespace OpenMM
#endif /*OPENMM_REFERENCETABULATEDFUNCTION_H_*/
platforms/reference/src/ReferenceKernels.cpp
View file @ 708c4246
......@@ -55,6 +55,7 @@
#include "ReferenceProperDihedralBond.h"
#include "ReferenceRbDihedralBond.h"
#include "ReferenceStochasticDynamics.h"
#include "ReferenceTabulatedFunction.h"
#include "ReferenceVariableStochasticDynamics.h"
#include "ReferenceVariableVerletDynamics.h"
#include "ReferenceVerletDynamics.h"
......@@ -69,7 +70,6 @@
#include "openmm/internal/CustomNonbondedForceImpl.h"
#include "openmm/internal/CMAPTorsionForceImpl.h"
#include "openmm/internal/NonbondedForceImpl.h"
#include "openmm/internal/SplineFitter.h"
#include "openmm/Integrator.h"
#include "openmm/OpenMMException.h"
#include "SimTKOpenMMUtilities.h"
......@@ -923,38 +923,6 @@ void ReferenceCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& con
dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
}
class ReferenceTabulatedFunction : public Lepton::CustomFunction {
public:
ReferenceTabulatedFunction(double min, double max, const vector<double>& values) :
min(min), max(max), values(values) {
int numValues = values.size();
x.resize(numValues);
for (int i = 0; i < numValues; i++)
x[i] = min+i*(max-min)/(numValues-1);
SplineFitter::createNaturalSpline(x, values, derivs);
}
int getNumArguments() const {
return 1;
}
double evaluate(const double* arguments) const {
double t = arguments[0];
if (t < min || t > max)
return 0.0;
return SplineFitter::evaluateSpline(x, values, derivs, t);
}
double evaluateDerivative(const double* arguments, const int* derivOrder) const {
double t = arguments[0];
if (t < min || t > max)
return 0.0;
return SplineFitter::evaluateSplineDerivative(x, values, derivs, t);
}
CustomFunction* clone() const {
return new ReferenceTabulatedFunction(min, max, values);
}
double min, max;
vector<double> x, values, derivs;
};
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
disposeRealArray(particleParamArray, numParticles);
if (neighborList != NULL)
......@@ -1001,13 +969,8 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
// Create custom functions for the tabulated functions.
map<string, Lepton::CustomFunction*> functions;
for (int i = 0; i < force.getNumFunctions(); i++) {
string name;
vector<double> values;
double min, max;
force.getFunctionParameters(i, name, values, min, max);
functions[name] = new ReferenceTabulatedFunction(min, max, values);
}
for (int i = 0; i < force.getNumFunctions(); i++)
functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
// Parse the various expressions used to calculate the force.
......@@ -1288,13 +1251,8 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
// Create custom functions for the tabulated functions.
map<string, Lepton::CustomFunction*> functions;
for (int i = 0; i < force.getNumFunctions(); i++) {
string name;
vector<double> values;
double min, max;
force.getFunctionParameters(i, name, values, min, max);
functions[name] = new ReferenceTabulatedFunction(min, max, values);
}
for (int i = 0; i < force.getNumFunctions(); i++)
functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
// Parse the expressions for computed values.
......@@ -1507,13 +1465,8 @@ void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const
// Create custom functions for the tabulated functions.
map<string, Lepton::CustomFunction*> functions;
for (int i = 0; i < force.getNumFunctions(); i++) {
string name;
vector<double> values;
double min, max;
force.getFunctionParameters(i, name, values, min, max);
functions[name] = new ReferenceTabulatedFunction(min, max, values);
}
for (int i = 0; i < force.getNumFunctions(); i++)
functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
// Parse the expression and create the object used to calculate the interaction.
......@@ -1609,13 +1562,8 @@ void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system
// Create custom functions for the tabulated functions.
map<string, Lepton::CustomFunction*> functions;
for (int i = 0; i < force.getNumFunctions(); i++) {
string name;
vector<double> values;
double min, max;
force.getFunctionParameters(i, name, values, min, max);
functions[name] = new ReferenceTabulatedFunction(min, max, values);
}
for (int i = 0; i < force.getNumFunctions(); i++)
functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
// Parse the expression and create the object used to calculate the interaction.
......
platforms/reference/src/ReferenceTabulatedFunction.cpp 0 → 100644
View file @ 708c4246
/* -------------------------------------------------------------------------- *
* 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 "ReferenceTabulatedFunction.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/SplineFitter.h"
#include <cmath>
using namespace OpenMM;
using namespace std;
using Lepton::CustomFunction;
extern "C" CustomFunction* createReferenceTabulatedFunction(const TabulatedFunction& function) {
if (dynamic_cast<const Continuous1DFunction*>(&function) != NULL)
return new ReferenceContinuous1DFunction(dynamic_cast<const Continuous1DFunction&>(function));
if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL)
return new ReferenceContinuous2DFunction(dynamic_cast<const Continuous2DFunction&>(function));
if (dynamic_cast<const Continuous3DFunction*>(&function) != NULL)
return new ReferenceContinuous3DFunction(dynamic_cast<const Continuous3DFunction&>(function));
if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL)
return new ReferenceDiscrete1DFunction(dynamic_cast<const Discrete1DFunction&>(function));
if (dynamic_cast<const Discrete2DFunction*>(&function) != NULL)
return new ReferenceDiscrete2DFunction(dynamic_cast<const Discrete2DFunction&>(function));
if (dynamic_cast<const Discrete3DFunction*>(&function) != NULL)
return new ReferenceDiscrete3DFunction(dynamic_cast<const Discrete3DFunction&>(function));
throw OpenMMException("createReferenceTabulatedFunction: Unknown function type");
}
ReferenceContinuous1DFunction::ReferenceContinuous1DFunction(const Continuous1DFunction& function) : function(function) {
function.getFunctionParameters(values, min, max);
int numValues = values.size();
x.resize(numValues);
for (int i = 0; i < numValues; i++)
x[i] = min+i*(max-min)/(numValues-1);
SplineFitter::createNaturalSpline(x, values, derivs);
}
int ReferenceContinuous1DFunction::getNumArguments() const {
return 1;
}
double ReferenceContinuous1DFunction::evaluate(const double* arguments) const {
double t = arguments[0];
if (t < min || t > max)
return 0.0;
return SplineFitter::evaluateSpline(x, values, derivs, t);
}
double ReferenceContinuous1DFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
double t = arguments[0];
if (t < min || t > max)
return 0.0;
return SplineFitter::evaluateSplineDerivative(x, values, derivs, t);
}
CustomFunction* ReferenceContinuous1DFunction::clone() const {
return new ReferenceContinuous1DFunction(function);
}
ReferenceContinuous2DFunction::ReferenceContinuous2DFunction(const Continuous2DFunction& function) : function(function) {
function.getFunctionParameters(xsize, ysize, values, xmin, xmax, ymin, ymax);
x.resize(xsize);
y.resize(ysize);
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);
SplineFitter::create2DNaturalSpline(x, y, values, c);
}
int ReferenceContinuous2DFunction::getNumArguments() const {
return 2;
}
double ReferenceContinuous2DFunction::evaluate(const double* arguments) const {
double u = arguments[0];
if (u < xmin || u > xmax)
return 0.0;
double v = arguments[1];
if (v < ymin || v > ymax)
return 0.0;
return SplineFitter::evaluate2DSpline(x, y, values, c, u, v);
}
double ReferenceContinuous2DFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
double u = arguments[0];
if (u < xmin || u > xmax)
return 0.0;
double v = arguments[1];
if (v < ymin || v > ymax)
return 0.0;
double dx, dy;
SplineFitter::evaluate2DSplineDerivatives(x, y, values, c, u, v, dx, dy);
if (derivOrder[0] == 1 && derivOrder[1] == 0)
return dx;
if (derivOrder[0] == 0 && derivOrder[1] == 1)
return dy;
throw OpenMMException("ReferenceContinuous2DFunction: Unsupported derivative order");
}
CustomFunction* ReferenceContinuous2DFunction::clone() const {
return new ReferenceContinuous2DFunction(function);
}
ReferenceContinuous3DFunction::ReferenceContinuous3DFunction(const Continuous3DFunction& function) : function(function) {
function.getFunctionParameters(xsize, ysize, zsize, values, xmin, xmax, ymin, ymax, zmin, zmax);
x.resize(xsize);
y.resize(ysize);
z.resize(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);
SplineFitter::create3DNaturalSpline(x, y, z, values, c);
}
int ReferenceContinuous3DFunction::getNumArguments() const {
return 3;
}
double ReferenceContinuous3DFunction::evaluate(const double* arguments) const {
double u = arguments[0];
if (u < xmin || u > xmax)
return 0.0;
double v = arguments[1];
if (v < ymin || v > ymax)
return 0.0;
double w = arguments[2];
if (w < zmin || w > zmax)
return 0.0;
return SplineFitter::evaluate3DSpline(x, y, z, values, c, u, v, w);
}
double ReferenceContinuous3DFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
double u = arguments[0];
if (u < xmin || u > xmax)
return 0.0;
double v = arguments[1];
if (v < ymin || v > ymax)
return 0.0;
double w = arguments[2];
if (w < zmin || w > zmax)
return 0.0;
double dx, dy, dz;
SplineFitter::evaluate3DSplineDerivatives(x, y, z, values, c, u, v, w, dx, dy, dz);
if (derivOrder[0] == 1 && derivOrder[1] == 0 && derivOrder[2] == 0)
return dx;
if (derivOrder[0] == 0 && derivOrder[1] == 1 && derivOrder[2] == 0)
return dy;
if (derivOrder[0] == 0 && derivOrder[1] == 0 && derivOrder[2] == 1)
return dz;
throw OpenMMException("ReferenceContinuous3DFunction: Unsupported derivative order");
}
CustomFunction* ReferenceContinuous3DFunction::clone() const {
return new ReferenceContinuous3DFunction(function);
}
ReferenceDiscrete1DFunction::ReferenceDiscrete1DFunction(const Discrete1DFunction& function) : function(function) {
function.getFunctionParameters(values);
}
int ReferenceDiscrete1DFunction::getNumArguments() const {
return 1;
}
double ReferenceDiscrete1DFunction::evaluate(const double* arguments) const {
int i = (int) round(arguments[0]);
if (i < 0 || i >= values.size())
throw OpenMMException("ReferenceDiscrete1DFunction: argument out of range");
return values[i];
}
double ReferenceDiscrete1DFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
return 0.0;
}
CustomFunction* ReferenceDiscrete1DFunction::clone() const {
return new ReferenceDiscrete1DFunction(function);
}
ReferenceDiscrete2DFunction::ReferenceDiscrete2DFunction(const Discrete2DFunction& function) : function(function) {
function.getFunctionParameters(xsize, ysize, values);
}
int ReferenceDiscrete2DFunction::getNumArguments() const {
return 2;
}
double ReferenceDiscrete2DFunction::evaluate(const double* arguments) const {
int i = (int) round(arguments[0]);
int j = (int) round(arguments[1]);
if (i < 0 || i >= xsize || j < 0 || j >= ysize)
throw OpenMMException("ReferenceDiscrete2DFunction: argument out of range");
return values[i+j*xsize];
}
double ReferenceDiscrete2DFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
return 0.0;
}
CustomFunction* ReferenceDiscrete2DFunction::clone() const {
return new ReferenceDiscrete2DFunction(function);
}
ReferenceDiscrete3DFunction::ReferenceDiscrete3DFunction(const Discrete3DFunction& function) : function(function) {
function.getFunctionParameters(xsize, ysize, zsize, values);
}
int ReferenceDiscrete3DFunction::getNumArguments() const {
return 3;
}
double ReferenceDiscrete3DFunction::evaluate(const double* arguments) const {
int i = (int) round(arguments[0]);
int j = (int) round(arguments[1]);
int k = (int) round(arguments[2]);
if (i < 0 || i >= xsize || j < 0 || j >= ysize || k < 0 || k >= zsize)
throw OpenMMException("ReferenceDiscrete3DFunction: argument out of range");
return values[i+(j+k*ysize)*xsize];
}
double ReferenceDiscrete3DFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
return 0.0;
}
CustomFunction* ReferenceDiscrete3DFunction::clone() const {
return new ReferenceDiscrete3DFunction(function);
}
platforms/reference/tests/TestReferenceCustomCompoundBondForce.cpp
View file @ 708c4246
......@@ -167,10 +167,111 @@ void testPositionDependence() {
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->addTabulatedFunction("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->addTabulatedFunction("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() {
try {
testBond();
testPositionDependence();
testContinuous2DFunction();
testContinuous3DFunction();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/reference/tests/TestReferenceCustomGBForce.cpp
View file @ 708c4246
......@@ -281,7 +281,7 @@ void testTabulatedFunction() {
vector<double> table;
for (int i = 0; i < 21; i++)
table.push_back(std::sin(0.25*i));
force->addFunction("fn", table, 1.0, 6.0);
force->addTabulatedFunction("fn", new Continuous1DFunction(table, 1.0, 6.0));
system.addForce(force);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
......
platforms/reference/tests/TestReferenceCustomHbondForce.cpp
View file @ 708c4246
......@@ -217,7 +217,7 @@ void testCustomFunctions() {
vector<double> function(2);
function[0] = 0;
function[1] = 1;
custom->addFunction("foo", function, 0, 10);
custom->addTabulatedFunction("foo", new Continuous1DFunction(function, 0, 10));
system.addForce(custom);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
......
platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp
View file @ 708c4246
......@@ -7,7 +7,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: *
* *
......@@ -227,7 +227,7 @@ void testPeriodic() {
ASSERT_EQUAL_TOL(1.9+1+0.9, state.getPotentialEnergy(), TOL);
}
void testTabulatedFunction() {
void testContinuous1DFunction() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
......@@ -238,21 +238,20 @@ void testTabulatedFunction() {
forceField->addParticle(vector<double>());
vector<double> table;
for (int i = 0; i < 21; i++)
table.push_back(std::sin(0.25*i));
forceField->addFunction("fn", table, 1.0, 6.0);
table.push_back(sin(0.25*i));
forceField->addTabulatedFunction("fn", new Continuous1DFunction(table, 1.0, 6.0));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
double tol = 0.01;
for (int i = 1; i < 30; i++) {
double x = (7.0/30.0)*i;
positions[1] = Vec3(x, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double force = (x < 1.0 || x > 6.0 ? 0.0 : -std::cos(x-1.0));
double energy = (x < 1.0 || x > 6.0 ? 0.0 : std::sin(x-1.0))+1.0;
double force = (x < 1.0 || x > 6.0 ? 0.0 : -cos(x-1.0));
double energy = (x < 1.0 || x > 6.0 ? 0.0 : sin(x-1.0))+1.0;
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.02);
......@@ -262,11 +261,217 @@ void testTabulatedFunction() {
positions[1] = Vec3(x, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double energy = (x < 1.0 || x > 6.0 ? 0.0 : std::sin(x-1.0))+1.0;
double energy = (x < 1.0 || x > 6.0 ? 0.0 : sin(x-1.0))+1.0;
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 1e-4);
}
}
void testContinuous2DFunction() {
const int xsize = 20;
const int ysize = 21;
const double xmin = 0.4;
const double xmax = 1.5;
const double ymin = 0.0;
const double ymax = 2.1;
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r,a)+1");
forceField->addGlobalParameter("a", 0.0);
forceField->addParticle(vector<double>());
forceField->addParticle(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->addTabulatedFunction("fn", new Continuous2DFunction(xsize, ysize, table, xmin, xmax, ymin, ymax));
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) {
positions[1] = Vec3(x, 0, 0);
context.setParameter("a", y);
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) {
energy = sin(0.25*x)*cos(0.33*y)+1.0;
force = -0.25*cos(0.25*x)*cos(0.33*y);
}
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.02);
}
}
}
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);
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->addTabulatedFunction("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() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r)+1");
forceField->addParticle(vector<double>());
forceField->addParticle(vector<double>());
vector<double> table;
for (int i = 0; i < 21; i++)
table.push_back(sin(0.25*i));
forceField->addTabulatedFunction("fn", new Discrete1DFunction(table));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
for (int i = 0; i < (int) table.size(); i++) {
positions[1] = Vec3(i, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[0], 1e-6);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], 1e-6);
ASSERT_EQUAL(table[i]+1.0, state.getPotentialEnergy());
}
}
void testDiscrete2DFunction() {
const int xsize = 10;
const int ysize = 5;
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("fn(r,a)+1");
forceField->addGlobalParameter("a", 0.0);
forceField->addParticle(vector<double>());
forceField->addParticle(vector<double>());
vector<double> table;
for (int i = 0; i < xsize; i++)
for (int j = 0; j < ysize; j++)
table.push_back(sin(0.25*i)+cos(0.33*j));
forceField->addTabulatedFunction("fn", new Discrete2DFunction(xsize, ysize, table));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
for (int i = 0; i < (int) table.size(); i++) {
positions[1] = Vec3(i%xsize, 0, 0);
context.setPositions(positions);
context.setParameter("a", i/xsize);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[0], 1e-6);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], 1e-6);
ASSERT_EQUAL(table[i]+1.0, state.getPotentialEnergy());
}
}
void testDiscrete3DFunction() {
const int xsize = 8;
const int ysize = 5;
const int zsize = 6;
ReferencePlatform platform;
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;
for (int i = 0; i < xsize; i++)
for (int j = 0; j < ysize; j++)
for (int k = 0; k < zsize; k++)
table.push_back(sin(0.25*i)+cos(0.33*j)+0.12345*k);
forceField->addTabulatedFunction("fn", new Discrete3DFunction(xsize, ysize, zsize, table));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
for (int i = 0; i < (int) table.size(); i++) {
positions[1] = Vec3(i%xsize, 0, 0);
context.setPositions(positions);
context.setParameter("a", (i/xsize)%ysize);
context.setParameter("b", i/(xsize*ysize));
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[0], 1e-6);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], 1e-6);
ASSERT_EQUAL(table[i]+1.0, state.getPotentialEnergy());
}
}
void testCoulombLennardJones() {
const int numMolecules = 300;
const int numParticles = numMolecules*2;
......@@ -658,7 +863,12 @@ int main() {
testExclusions();
testCutoff();
testPeriodic();
testTabulatedFunction();
testContinuous1DFunction();
testContinuous2DFunction();
testContinuous3DFunction();
testDiscrete1DFunction();
testDiscrete2DFunction();
testDiscrete3DFunction();
testCoulombLennardJones();
testSwitchingFunction();
testLongRangeCorrection();
......
plugins/amoeba/CMakeLists.txt
View file @ 708c4246
......@@ -101,7 +101,7 @@ SET_SOURCE_FILES_PROPERTIES(${serialization_files} PROPERTIES COMPILE_FLAGS "-DO
IF(OPENMM_BUILD_STATIC_LIB)
ADD_LIBRARY(${STATIC_AMOEBA_TARGET} STATIC ${SOURCE_AMOEBA_FILES} ${SOURCE_AMOEBA_INCLUDE_FILES} ${API_AMOEBA_ABS_INCLUDE_FILES})
SET_TARGET_PROPERTIES(${STATIC_AMOEBA_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_USE_STATIC_LIBRARIES -DOPENMM_BUILDING_STATIC_LIBRARY -DLEPTON_USE_STATIC_LIBRARIES -DLEPTON_BUILDING_STATIC_LIBRARY")
SET_TARGET_PROPERTIES(${STATIC_AMOEBA_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_AMOEBA_BUILDING_STATIC_LIBRARY -DOPENMM_USE_STATIC_LIBRARIES -DOPENMM_BUILDING_STATIC_LIBRARY -DLEPTON_USE_STATIC_LIBRARIES -DLEPTON_BUILDING_STATIC_LIBRARY")
ENDIF(OPENMM_BUILD_STATIC_LIB)
IF(OPENMM_BUILD_C_AND_FORTRAN_WRAPPERS)
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 708c4246
......@@ -893,8 +893,9 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
numMultipoles = force.getNumMultipoles();
CudaArray& posq = cu.getPosq();
float4* posqf = (float4*) cu.getPinnedBuffer();
double4* posqd = (double4*) cu.getPinnedBuffer();
vector<double4> temp(posq.getSize());
float4* posqf = (float4*) &temp[0];
double4* posqd = (double4*) &temp[0];
vector<float2> dampingAndTholeVec;
vector<float> polarizabilityVec;
vector<float> molecularDipolesVec;
......@@ -945,7 +946,7 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
multipoleParticles->upload(multipoleParticlesVec);
molecularDipoles->upload(molecularDipolesVec);
molecularQuadrupoles->upload(molecularQuadrupolesVec);
posq.upload(cu.getPinnedBuffer());
posq.upload(&temp[0]);
// Create workspace arrays.
......
plugins/amoeba/wrappers/generateAmoebaWrappers.py
View file @ 708c4246
......@@ -550,6 +550,8 @@ class CSourceGenerator(WrapperGenerator):
unwrappedType = type[:-1].strip()
if unwrappedType in self.classesByShortName:
unwrappedType = self.classesByShortName[unwrappedType]
if unwrappedType == 'const std::string':
return 'std::string(%s)' % value
return '*'+self.unwrapValue(unwrappedType+'*', value)
if type in self.classesByShortName:
return 'static_cast<%s>(%s)' % (self.classesByShortName[type], value)
......@@ -1134,14 +1136,14 @@ class FortranSourceGenerator(WrapperGenerator):
return type
def isHandleType(self, type):
if type.startswith('OpenMM_'):
return True;
if type == 'Vec3':
return True
if type == 'OpenMM_Vec3':
return False
if type.endswith('*') or type.endswith('&'):
return self.isHandleType(type[:-1].strip())
if type.startswith('const '):
return self.isHandleType(type[6:].strip())
if type.startswith('OpenMM_'):
return True;
return False
def writeOutput(self):
......
plugins/cpupme/CMakeLists.txt
View file @ 708c4246
......@@ -22,7 +22,7 @@ SET(OPENMM_SOURCE_SUBDIRS .)
SET(OPENMMPME_LIBRARY_NAME OpenMMPME)
SET(SHARED_TARGET ${OPENMMPME_LIBRARY_NAME})
SET(STATIC_TARGET ${OPENMMPME_LIBRARY_NAME}_static)
# Ensure that debug libraries have "_d" appended to their names.
# CMake gets this right on Windows automatically with this definition.
......@@ -74,7 +74,7 @@ ENDIF (NOT MSVC)
# Include FFTW related files.
INCLUDE_DIRECTORIES(${FFTW_INCLUDES})
# Build the plugin library.
# Build the shared plugin library.
ADD_LIBRARY(${SHARED_TARGET} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_INCLUDE_FILES})
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
......@@ -90,4 +90,22 @@ SET_TARGET_PROPERTIES(${SHARED_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FL
INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${SHARED_TARGET})
# Build the static plugin library.
IF(OPENMM_BUILD_STATIC_LIB)
ADD_LIBRARY(${STATIC_TARGET} STATIC ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_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} ${FFTW_LIBRARY})
IF (FFTW_THREADS_LIBRARY)
TARGET_LINK_LIBRARIES(${STATIC_TARGET} ${FFTW_THREADS_LIBRARY})
ENDIF (FFTW_THREADS_LIBRARY)
SET_TARGET_PROPERTIES(${STATIC_TARGET} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} -DOPENMM_PME_BUILDING_STATIC_LIBRARY")
INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${STATIC_TARGET})
ENDIF(OPENMM_BUILD_STATIC_LIB)
SUBDIRS(tests)
plugins/cpupme/src/CpuPmeKernelFactory.cpp
View file @ 708c4246
......@@ -32,9 +32,6 @@
using namespace OpenMM;
extern "C" OPENMM_EXPORT_PME void registerPlatforms() {
}
extern "C" OPENMM_EXPORT_PME void registerKernelFactories() {
if (CpuCalcPmeReciprocalForceKernel::isProcessorSupported()) {
CpuPmeKernelFactory* factory = new CpuPmeKernelFactory();
......@@ -43,9 +40,17 @@ extern "C" OPENMM_EXPORT_PME void registerKernelFactories() {
}
}
#ifdef OPENMM_PME_BUILDING_STATIC_LIBRARY
extern "C" void registerCpuPmeKernelFactories() {
registerKernelFactories();
}
#else
extern "C" OPENMM_EXPORT_PME void registerCpuPmeKernelFactories() {
registerKernelFactories();
}
extern "C" OPENMM_EXPORT_PME void registerPlatforms() {
}
#endif
KernelImpl* CpuPmeKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
if (name == CalcPmeReciprocalForceKernel::Name())
......
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