Skip to content

GitLab

Commit 4233befd authored by leeping's avatar leeping
Browse files

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

parents eb9f735a 31d02cdc
Hide whitespace changes
Inline Side-by-side
Showing with 702 additions and 3963 deletions
openmmapi/src/ThreadPool.cpp
View file @ 4233befd
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013 Stanford University and the Authors. *
* Portions copyright (c) 2013-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -60,8 +60,10 @@ static void* threadBody(void* args) {
return 0;
}
ThreadPool::ThreadPool() {
numThreads = getNumProcessors();
ThreadPool::ThreadPool(int numThreads) {
if (numThreads <= 0)
numThreads = getNumProcessors();
this->numThreads = numThreads;
pthread_cond_init(&startCondition, NULL);
pthread_cond_init(&endCondition, NULL);
pthread_mutex_init(&lock, NULL);
......
platforms/cpu/include/CpuPlatform.h
View file @ 4233befd
......@@ -55,27 +55,37 @@ public:
return name;
}
double getSpeed() const;
const std::string& getPropertyValue(const Context& context, const std::string& property) const;
bool supportsDoublePrecision() const;
static bool isProcessorSupported();
void contextCreated(ContextImpl& context, const std::map<std::string, std::string>& properties) const;
void contextDestroyed(ContextImpl& context) const;
/**
* This is the name of the parameter for selecting the number of threads to use.
*/
static const std::string& CpuThreads() {
static const std::string key = "CpuThreads";
return key;
}
/**
* We cannot use the standard mechanism for platform data, because that is already used by the superclass.
* Instead, we maintain a table of ContextImpls to PlatformDatas.
*/
static PlatformData& getPlatformData(ContextImpl& context);
static const PlatformData& getPlatformData(const ContextImpl& context);
private:
static std::map<ContextImpl*, PlatformData*> contextData;
static std::map<const ContextImpl*, PlatformData*> contextData;
};
class CpuPlatform::PlatformData {
public:
PlatformData(int numParticles);
PlatformData(int numParticles, int numThreads);
AlignedArray<float> posq;
std::vector<AlignedArray<float> > threadForce;
ThreadPool threads;
bool isPeriodic;
CpuRandom random;
std::map<std::string, std::string> propertyValues;
};
} // namespace OpenMM
......
platforms/cpu/src/CpuPlatform.cpp
View file @ 4233befd
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013 Stanford University and the Authors. *
* Portions copyright (c) 2013-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -35,6 +35,7 @@
#include "CpuSETTLE.h"
#include "ReferenceConstraints.h"
#include "openmm/internal/hardware.h"
#include <sstream>
using namespace OpenMM;
using namespace std;
......@@ -53,7 +54,7 @@ extern "C" OPENMM_EXPORT_CPU void registerPlatforms() {
}
#endif
map<ContextImpl*, CpuPlatform::PlatformData*> CpuPlatform::contextData;
map<const ContextImpl*, CpuPlatform::PlatformData*> CpuPlatform::contextData;
CpuPlatform::CpuPlatform() {
CpuKernelFactory* factory = new CpuKernelFactory();
......@@ -63,6 +64,23 @@ CpuPlatform::CpuPlatform() {
registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
platformProperties.push_back(CpuThreads());
int threads = getNumProcessors();
char* threadsEnv = getenv("OPENMM_CPU_THREADS");
if (threadsEnv != NULL)
stringstream(threadsEnv) >> threads;
stringstream defaultThreads;
defaultThreads << threads;
setPropertyDefaultValue(CpuThreads(), defaultThreads.str());
}
const string& CpuPlatform::getPropertyValue(const Context& context, const string& property) const {
const ContextImpl& impl = getContextImpl(context);
const PlatformData& data = getPlatformData(impl);
map<string, string>::const_iterator value = data.propertyValues.find(property);
if (value != data.propertyValues.end())
return value->second;
return ReferencePlatform::getPropertyValue(context, property);
}
double CpuPlatform::getSpeed() const {
......@@ -87,7 +105,11 @@ bool CpuPlatform::isProcessorSupported() {
void CpuPlatform::contextCreated(ContextImpl& context, const map<string, string>& properties) const {
ReferencePlatform::contextCreated(context, properties);
PlatformData* data = new PlatformData(context.getSystem().getNumParticles());
const string& threadsPropValue = (properties.find(CpuThreads()) == properties.end() ?
getPropertyDefaultValue(CpuThreads()) : properties.find(CpuThreads())->second);
int numThreads;
stringstream(threadsPropValue) >> numThreads;
PlatformData* data = new PlatformData(context.getSystem().getNumParticles(), numThreads);
contextData[&context] = data;
ReferenceConstraints& constraints = *(ReferenceConstraints*) reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData())->constraints;
if (constraints.settle != NULL) {
......@@ -107,10 +129,17 @@ CpuPlatform::PlatformData& CpuPlatform::getPlatformData(ContextImpl& context) {
return *contextData[&context];
}
CpuPlatform::PlatformData::PlatformData(int numParticles) : posq(4*numParticles) {
int numThreads = threads.getNumThreads();
const CpuPlatform::PlatformData& CpuPlatform::getPlatformData(const ContextImpl& context) {
return *contextData[&context];
}
CpuPlatform::PlatformData::PlatformData(int numParticles, int numThreads) : posq(4*numParticles), threads(numThreads) {
numThreads = threads.getNumThreads();
threadForce.resize(numThreads);
for (int i = 0; i < numThreads; i++)
threadForce[i].resize(4*numParticles);
isPeriodic = false;
stringstream threadsProperty;
threadsProperty << numThreads;
propertyValues[CpuThreads()] = threadsProperty.str();
}
platforms/cpu/src/CpuSETTLE.cpp
View file @ 4233befd
......@@ -41,7 +41,9 @@ public:
inverseMasses(inverseMasses), tolerance(tolerance), threadSettle(threadSettle) {
}
void execute(ThreadPool& threads, int threadIndex) {
threadSettle[threadIndex]->apply(atomCoordinates, atomCoordinatesP, inverseMasses, tolerance);
if (threadIndex < threadSettle.size()) {
threadSettle[threadIndex]->apply(atomCoordinates, atomCoordinatesP, inverseMasses, tolerance);
}
}
vector<OpenMM::RealVec>& atomCoordinates;
vector<OpenMM::RealVec>& atomCoordinatesP;
......
platforms/cuda/tests/CMakeLists.txt
View file @ 4233befd
......@@ -31,35 +31,6 @@ FOREACH(TEST_PROG ${TEST_PROGS})
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS}")
ENDIF (APPLE)
IF( ${TEST_ROOT} STREQUAL "TestCUDAGBSAOBCForce2" )
SET(DEFINE_STRING "-DTEST_PLATFORM=1 ")
IF( INCLUDE_SERIALIZATION )
SET(DEFINE_STRING "${DEFINE_STRING} -DOPENMM_SERIALIZE ")
ENDIF( INCLUDE_SERIALIZATION )
# nonbond
SET(NONBOND_DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=0")
SET(NONBOND_TEST "TestCUDANonbondedForce2")
ADD_EXECUTABLE(${NONBOND_TEST} ${TEST_PROG})
SET_TARGET_PROPERTIES(${NONBOND_TEST} PROPERTIES COMPILE_FLAGS ${NONBOND_DEFINE_STRING} )
ADD_TEST(${NONBOND_TEST} ${EXECUTABLE_OUTPUT_PATH}/${NONBOND_TEST})
# OBC
SET(DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=1")
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES COMPILE_FLAGS ${DEFINE_STRING} )
IF( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
ELSE( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET})
ENDIF( INCLUDE_SERIALIZATION )
ENDIF( ${TEST_ROOT} STREQUAL "TestCUDAGBSAOBCForce2" )
ADD_TEST(${TEST_ROOT}Single ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} single)
IF (OPENMM_BUILD_CUDA_DOUBLE_PRECISION_TESTS)
ADD_TEST(${TEST_ROOT}Mixed ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} mixed)
......
platforms/cuda/tests/TestCudaCustomCompoundBondForce.cpp
View file @ 4233befd
......@@ -296,6 +296,53 @@ void testContinuous3DFunction() {
}
}
void testMultipleBonds() {
// Two compound bonds using Urey-Bradley example from API doc
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(3,
"0.5*(kangle*(angle(p1,p2,p3)-theta0)^2+kbond*(distance(p1,p3)-r0)^2)");
custom->addPerBondParameter("kangle");
custom->addPerBondParameter("kbond");
custom->addPerBondParameter("theta0");
custom->addPerBondParameter("r0");
vector<double> parameters(4);
parameters[0] = 1.0;
parameters[1] = 1.0;
parameters[2] = 2 * M_PI / 3;
parameters[3] = sqrt(3) / 2;
vector<int> particles0(3);
particles0[0] = 0;
particles0[1] = 1;
particles0[2] = 2;
vector<int> particles1(3);
particles1[0] = 1;
particles1[1] = 2;
particles1[2] = 3;
custom->addBond(particles0, parameters);
custom->addBond(particles1, parameters);
customSystem.addForce(custom);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0.5, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0.5, 0, 0);
positions[3] = Vec3(0.6, 0, 0.4);
VerletIntegrator integrator(0.01);
Context context(customSystem, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(0.199, state.getPotentialEnergy(), 1e-3);
vector<Vec3> forces(state.getForces());
ASSERT_EQUAL_VEC(Vec3(-1.160, 0.112, 0.0), forces[0], 1e-3);
ASSERT_EQUAL_VEC(Vec3(0.927, 1.047, -0.638), forces[1], 1e-3);
ASSERT_EQUAL_VEC(Vec3(-0.543, -1.160, 0.721), forces[2], 1e-3);
ASSERT_EQUAL_VEC(Vec3(0.776, 0.0, -0.084), forces[3], 1e-3);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
......@@ -305,6 +352,7 @@ int main(int argc, char* argv[]) {
testParallelComputation();
testContinuous2DFunction();
testContinuous3DFunction();
testMultipleBonds();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/cuda/tests/TestCudaEwald.cpp
View file @ 4233befd
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -252,6 +252,54 @@ void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
}
}
void testPMEParameters() {
// Create a cloud of random point charges.
const int numParticles = 51;
const double boxWidth = 4.7;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
}
force->setNonbondedMethod(NonbondedForce::PME);
// Compute the energy with an error tolerance of 1e-3.
force->setEwaldErrorTolerance(1e-3);
VerletIntegrator integrator1(0.01);
Context context1(system, integrator1, platform);
context1.setPositions(positions);
double energy1 = context1.getState(State::Energy).getPotentialEnergy();
// Try again with an error tolerance of 1e-4.
force->setEwaldErrorTolerance(1e-4);
VerletIntegrator integrator2(0.01);
Context context2(system, integrator2, platform);
context2.setPositions(positions);
double energy2 = context2.getState(State::Energy).getPotentialEnergy();
// Now explicitly set the parameters. These should match the values that were
// used for tolerance 1e-3.
force->setPMEParameters(2.49291157051793, 32, 32, 32);
VerletIntegrator integrator3(0.01);
Context context3(system, integrator3, platform);
context3.setPositions(positions);
double energy3 = context3.getState(State::Energy).getPotentialEnergy();
ASSERT_EQUAL_TOL(energy1, energy3, 1e-6);
ASSERT(fabs((energy1-energy2)/energy1) > 1e-5);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
......@@ -261,6 +309,7 @@ int main(int argc, char* argv[]) {
// testEwald2Ions();
testErrorTolerance(NonbondedForce::Ewald);
testErrorTolerance(NonbondedForce::PME);
testPMEParameters();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/opencl/src/OpenCLFFT3D.cpp
View file @ 4233befd
......@@ -231,7 +231,6 @@ cl::Kernel OpenCLFFT3D::createKernel(int xsize, int ysize, int zsize, int& threa
source<<"if (index < XSIZE*YSIZE)\n";
source<<"out[y*(ZSIZE*XSIZE)+(get_local_id(0)%ZSIZE)*XSIZE+x] = data"<<(stage%2)<<"[get_local_id(0)];\n";
}
source<<"barrier(CLK_GLOBAL_MEM_FENCE);";
map<string, string> replacements;
replacements["XSIZE"] = context.intToString(xsize);
replacements["YSIZE"] = context.intToString(ysize);
......
platforms/opencl/tests/CMakeLists.txt
View file @ 4233befd
......@@ -27,35 +27,6 @@ FOREACH(TEST_PROG ${TEST_PROGS})
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET})
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS}")
IF( ${TEST_ROOT} STREQUAL "TestOpenCLGBSAOBCForce2" )
SET(DEFINE_STRING "-DTEST_PLATFORM=1 ")
IF( INCLUDE_SERIALIZATION )
SET(DEFINE_STRING "${DEFINE_STRING} -DOPENMM_SERIALIZE ")
ENDIF( INCLUDE_SERIALIZATION )
# nonbond
SET(NONBOND_DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=0")
SET(NONBOND_TEST "TestOpenCLNonbondedForce2")
ADD_EXECUTABLE(${NONBOND_TEST} ${TEST_PROG})
SET_TARGET_PROPERTIES(${NONBOND_TEST} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} ${NONBOND_DEFINE_STRING}" )
ADD_TEST(${NONBOND_TEST} ${EXECUTABLE_OUTPUT_PATH}/${NONBOND_TEST})
# OBC
SET(DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=1")
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES LINK_FLAGS "${EXTRA_COMPILE_FLAGS}" COMPILE_FLAGS "${EXTRA_COMPILE_FLAGS} ${DEFINE_STRING}" )
IF( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
ELSE( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET})
ENDIF( INCLUDE_SERIALIZATION )
ENDIF( ${TEST_ROOT} STREQUAL "TestOpenCLGBSAOBCForce2" )
ADD_TEST(${TEST_ROOT}Single ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} single)
IF (OPENMM_BUILD_OPENCL_DOUBLE_PRECISION_TESTS)
ADD_TEST(${TEST_ROOT}Mixed ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} mixed)
......
platforms/opencl/tests/TestOpenCLCustomCompoundBondForce.cpp
View file @ 4233befd
......@@ -296,6 +296,53 @@ void testContinuous3DFunction() {
}
}
void testMultipleBonds() {
// Two compound bonds using Urey-Bradley example from API doc
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(3,
"0.5*(kangle*(angle(p1,p2,p3)-theta0)^2+kbond*(distance(p1,p3)-r0)^2)");
custom->addPerBondParameter("kangle");
custom->addPerBondParameter("kbond");
custom->addPerBondParameter("theta0");
custom->addPerBondParameter("r0");
vector<double> parameters(4);
parameters[0] = 1.0;
parameters[1] = 1.0;
parameters[2] = 2 * M_PI / 3;
parameters[3] = sqrt(3) / 2;
vector<int> particles0(3);
particles0[0] = 0;
particles0[1] = 1;
particles0[2] = 2;
vector<int> particles1(3);
particles1[0] = 1;
particles1[1] = 2;
particles1[2] = 3;
custom->addBond(particles0, parameters);
custom->addBond(particles1, parameters);
customSystem.addForce(custom);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0.5, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0.5, 0, 0);
positions[3] = Vec3(0.6, 0, 0.4);
VerletIntegrator integrator(0.01);
Context context(customSystem, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(0.199, state.getPotentialEnergy(), 1e-3);
vector<Vec3> forces(state.getForces());
ASSERT_EQUAL_VEC(Vec3(-1.160, 0.112, 0.0), forces[0], 1e-3);
ASSERT_EQUAL_VEC(Vec3(0.927, 1.047, -0.638), forces[1], 1e-3);
ASSERT_EQUAL_VEC(Vec3(-0.543, -1.160, 0.721), forces[2], 1e-3);
ASSERT_EQUAL_VEC(Vec3(0.776, 0.0, -0.084), forces[3], 1e-3);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
......@@ -305,6 +352,7 @@ int main(int argc, char* argv[]) {
testParallelComputation();
testContinuous2DFunction();
testContinuous3DFunction();
testMultipleBonds();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/opencl/tests/TestOpenCLEwald.cpp
View file @ 4233befd
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2010 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -252,6 +252,54 @@ void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
}
}
void testPMEParameters() {
// Create a cloud of random point charges.
const int numParticles = 51;
const double boxWidth = 4.7;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
}
force->setNonbondedMethod(NonbondedForce::PME);
// Compute the energy with an error tolerance of 1e-3.
force->setEwaldErrorTolerance(1e-3);
VerletIntegrator integrator1(0.01);
Context context1(system, integrator1, platform);
context1.setPositions(positions);
double energy1 = context1.getState(State::Energy).getPotentialEnergy();
// Try again with an error tolerance of 1e-4.
force->setEwaldErrorTolerance(1e-4);
VerletIntegrator integrator2(0.01);
Context context2(system, integrator2, platform);
context2.setPositions(positions);
double energy2 = context2.getState(State::Energy).getPotentialEnergy();
// Now explicitly set the parameters. These should match the values that were
// used for tolerance 1e-3.
force->setPMEParameters(2.49291157051793, 32, 32, 32);
VerletIntegrator integrator3(0.01);
Context context3(system, integrator3, platform);
context3.setPositions(positions);
double energy3 = context3.getState(State::Energy).getPotentialEnergy();
ASSERT_EQUAL_TOL(energy1, energy3, 1e-6);
ASSERT(fabs((energy1-energy2)/energy1) > 1e-5);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
......@@ -261,6 +309,7 @@ int main(int argc, char* argv[]) {
// testEwald2Ions();
testErrorTolerance(NonbondedForce::Ewald);
testErrorTolerance(NonbondedForce::PME);
testPMEParameters();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/opencl/tests/TestOpenCLGBSAOBCForce2.cpp deleted 100644 → 0
View file @ eb9f735a
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the OpenCL implementations of GBVIForce, GBSAOBCForce and the softcore versions of
* these forces.
*/
#define TEST_NONBONDED 0
#define TEST_OBC 1
#define TEST_GBVI 2
#define TEST_CUDA_PLATFORM 0
#define TEST_OPENCL_PLATFORM 1
#include "openmm/GBVIForce.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/NonbondedForce.h"
#ifdef USE_SOFTCORE
#include "openmm/GBVISoftcoreForce.h"
#include "openmm/GBSAOBCSoftcoreForce.h"
#include "openmm/NonbondedSoftcoreForce.h"
#endif
/**
* Utility methods shared across unit tests
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/System.h"
#include "SimTKOpenMMRealType.h"
#include "OpenMM.h"
#if TEST_PLATFORM == TEST_OPENCL_PLATFORM
#include "ReferencePlatform.h"
#include "OpenCLPlatform.h"
#endif
#if TEST_PLATFORM == TEST_CUDA_PLATFORM
#include "ReferencePlatform.h"
#include "CudaPlatform.h"
#endif
#ifdef USE_SOFTCORE
#include "OpenMMFreeEnergy.h"
#include "openmm/freeEnergyKernels.h"
#endif
#include "sfmt/SFMT.h"
#include "openmm/VerletIntegrator.h"
#ifdef OPENMM_SERIALIZE
#include "openmm/serialization/SerializationProxy.h"
#include "openmm/serialization/SerializationNode.h"
#include "openmm/serialization/XmlSerializer.h"
#endif
#include <iostream>
#include <vector>
#include <algorithm>
#include <cstdio>
#include <typeinfo>
#include <iomanip>
extern "C" void registerFreeEnergyCudaKernelFactories();
using namespace OpenMM;
using namespace std;
static const int NoCutoff_OpenMMTest = 0;
static const int CutoffNonPeriodic_OpenMMTest = 1;
static const int CutoffPeriodic_OpenMMTest = 2;
static const int Ewald_OpenMMTest = 3;
static const int PME_OpenMMTest = 4;
static const int ChargeIndex_OpenMMTest = 0;
static const int SigmaIndex_OpenMMTest = 1;
static const int EpsIndex_OpenMMTest = 2;
static const int GammaIndex_OpenMMTest = 3;
static const int LambdaIndex_OpenMMTest = 4;
static const int Reference_OpenMMTest = 0;
static const int Cuda_OpenMMTest = 1;
static const int OpenCL_OpenMMTest = 2;
class BondInfo_OpenMMTest {
public:
BondInfo_OpenMMTest( int particle1, int particle2, double distance );
int _particle1;
int _particle2;
double _distance;
};
BondInfo_OpenMMTest::BondInfo_OpenMMTest( int particle1, int particle2, double distance ){
_particle1 = particle1;
_particle2 = particle2;
_distance = distance;
}
typedef std::vector<std::string> StringVector;
typedef StringVector::iterator StringVectorI;
typedef StringVector::const_iterator StringVectorCI;
typedef std::vector<StringVector> StringVectorVector;
typedef std::vector<int> IntVector;
typedef IntVector::iterator IntVectorI;
typedef IntVector::const_iterator IntVectorCI;
typedef std::vector<IntVector> IntVectorVector;
typedef std::vector<double> DoubleVector;
typedef DoubleVector::iterator DoubleVectorI;
typedef DoubleVector::const_iterator DoubleVectorCI;
typedef std::vector<DoubleVector> DoubleVectorVector;
// the following are used in parsing parameter file
typedef std::vector<std::string> StringVector;
typedef StringVector::iterator StringVectorI;
typedef StringVector::const_iterator StringVectorCI;
typedef std::vector<StringVector> VectorStringVector;
typedef VectorStringVector::iterator VectorStringVectorI;
typedef VectorStringVector::const_iterator VectorStringVectorCI;
typedef std::vector<std::vector<double> > VectorOfDoubleVectors;
typedef VectorOfDoubleVectors::iterator VectorOfDoubleVectorsI;
typedef VectorOfDoubleVectors::const_iterator VectorOfDoubleVectorsCI;
typedef std::map< int, int> MapIntInt;
typedef MapIntInt::iterator MapIntIntI;
typedef MapIntInt::const_iterator MapIntIntCI;
typedef std::map< double, int> MapDoubleToInt;
typedef MapDoubleToInt::iterator MapDoubleToIntI;
typedef MapDoubleToInt::const_iterator MapDoubleToIntCI;
typedef std::map< std::string, VectorOfDoubleVectors > MapStringVectorOfDoubleVectors;
typedef MapStringVectorOfDoubleVectors::iterator MapStringVectorOfDoubleVectorsI;
typedef MapStringVectorOfDoubleVectors::const_iterator MapStringVectorOfDoubleVectorsCI;
typedef std::map< std::string, StringVector > MapStringStringVector;
typedef MapStringStringVector::iterator MapStringStringVectorI;
typedef MapStringStringVector::const_iterator MapStringStringVectorCI;
typedef std::map< std::string, std::string > MapStringString;
typedef MapStringString::iterator MapStringStringI;
typedef MapStringString::const_iterator MapStringStringCI;
typedef std::map< std::string, std::string > MapStringToInt;
typedef MapStringToInt::iterator MapStringToIntI;
typedef MapStringToInt::const_iterator MapStringToIntCI;
typedef std::vector< std::map< std::string, std::string > > VectorMapStringString;
typedef VectorMapStringString::iterator VectorMapStringStringI;
typedef VectorMapStringString::const_iterator VectorMapStringStringCI;
typedef std::map< std::string, int > MapStringInt;
typedef MapStringInt::iterator MapStringIntI;
typedef MapStringInt::const_iterator MapStringIntCI;
typedef std::map< std::string, std::vector<OpenMM::Vec3> > MapStringVec3;
typedef MapStringVec3::iterator MapStringVec3I;
typedef MapStringVec3::const_iterator MapStringVec3CI;
typedef std::map< std::string, double > MapStringToDouble;
typedef MapStringToDouble::iterator MapStringToDoubleI;
typedef MapStringToDouble::const_iterator MapStringToDoubleCI;
typedef std::vector< MapStringToDouble > VectorOfMapStringToDouble;
typedef std::map< std::string, DoubleVector> MapStringToDoubleVector;
typedef MapStringToDoubleVector::iterator MapStringToDoubleVectorI;
typedef MapStringToDoubleVector::const_iterator MapStringToDoubleVectorCI;
typedef std::map< std::string, DoubleVector > MapStringToDoubleVector;
typedef std::map< std::string, Force*> MapStringForce;
typedef MapStringForce::iterator MapStringForceI;
typedef MapStringForce::const_iterator MapStringForceCI;
typedef std::map< int, IntVector> MapIntIntVector;
typedef MapIntIntVector::const_iterator MapIntIntVectorCI;
typedef std::pair<int, int> IntIntPair;
typedef std::vector<IntIntPair> IntIntPairVector;
typedef IntIntPairVector::iterator IntIntPairVectorI;
typedef IntIntPairVector::const_iterator IntIntPairVectorCI;
typedef std::pair<int, double> IntDoublePair;
typedef std::vector<IntDoublePair> IntDoublePairVector;
typedef IntDoublePairVector::iterator IntDoublePairVectorI;
typedef IntDoublePairVector::const_iterator IntDoublePairVectorCI;
/**
* Predicate for sorting <int,double> pair
*
* @param d1 first IntDoublePair to compare
* @param d2 second IntDoublePair to compare
*
*/
bool TestIntDoublePair( const IntDoublePair& d1, const IntDoublePair& d2 ){
return d1.second < d2.second;
}
class PositionGenerator {
public:
enum GenerationMethod {
/**
* Random positions
*/
Random = 0,
/**
* On grid
*/
SimpleGrid = 1,
};
PositionGenerator( int numMolecules, int numParticlesPerMolecule, double boxSize );
~PositionGenerator( );
/**
* Set logging file reference
*
* @param log log
*
*/
void setLog( FILE* log );
/**
* Set positions
*
* @param method method placement
* @param positions output vector of positions
*
* @return nonzero if error detected; 0 otherwise
*/
int setPositions( GenerationMethod method, std::vector<Vec3>& positions ) const;
/**
* Set positions
*
* @param method method placement
* @param sfmt input random number generator
* @param positions output vector of positions
*
* @return nonzero if error detected; 0 otherwise
*/
int setPositions( GenerationMethod method, OpenMM_SFMT::SFMT& sfmt, std::vector<Vec3>& positions ) const;
/**
* Place particles on a grid
*
* @param origin origin
* @param boxDimensions box dimensions
* @param spacing spacing
* @param sfmt input random number generator
* @param bondDistance input bond distance
* @param array output vector of grid values
*
* @return -1 if particles will not fit on grid; 0 if they do
*/
int setParticlesOnGrid( const Vec3& origin, const Vec3& boxDimensions, const Vec3& spacing,
OpenMM_SFMT::SFMT& sfmt, double bondDistance, std::vector<Vec3>& array ) const;
/**
* Set bond distance
*
* @param bondDistance bond distance
*/
void setBondDistance( double bondDistance );
/**
* Get bond distance
*
* @return bond distance
*/
double getBondDistance( void ) const;
/**
* Get distance
*
* @param index1 index of first particle
* @param index2 index of second particle
* @param positions particle positions
*
* @return distance
*/
double getDistance( int index1, int index2, const std::vector<Vec3>& positions ) const;
/**
* Get distance assumming periodic boundary conditions
*
* @param index1 index of first particle
* @param index2 index of second particle
* @param positions particle positions
*
* @return distance
*/
double getPeriodicDistance( int index1, int index2, const std::vector<Vec3>& positions ) const;
/**
* Get settings
*
* @return info string
*/
std::string getSettings( void ) const;
/**
* Get enclosing box
*
* @param positions input vector of positions
* @param enclosingBox output vector of enclosing box dimensions
*
*/
void getEnclosingBox( const std::vector<Vec3>& positions, Vec3 enclosingBox[2] ) const;
/**
* Get sorted distances from particular position
*
* @param periodicBoundaryConditions if set, apply PBC
* @param positionIndex input position index
* @param positions input vector of positions
* @param sortVector on output sorted IntDoublePairVector
*
*/
void getSortedDistances( int periodicBoundaryConditions, int positionIndex, const std::vector<Vec3>& positions, IntDoublePairVector& sortVector ) const;
/**
* Show min/max distances between all positions
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param showIndex number of min/maxentries to show
*
*/
void showMinMaxDistances( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, int showIndex );
/**
* Show min/max distances between a single position and all remaining positions
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param showIndex number of min/maxentries to show
* @param positionIndexVector list of entries to show min/max distances from
*
*/
void showMinMaxDistances( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, int showIndex,
const IntVector& positionIndexVector );
/**
* Show distances between positions
*
* @param pairs particle indcies for which distance is to be reported
* @param positions input vector of positions
*
*/
void showDistances( const IntIntPairVector& pairs, const std::vector<Vec3>& positions ) const;
/**
* Show particles within a specified distance of a given particle
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param particleIndex particle to check
* @param distanceToCheckFor distance to check for
* @param tolerance distance tolerance
*
*/
void showParticlesWithinDistance( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, unsigned int particleIndex,
double distanceToCheckFor, double tolerance);
private:
int _numMolecules;
int _numParticlesPerMolecule;
int _numParticles;
int _seed;
double _boxSize;
double _bondDistance;
Vec3 _origin;
Vec3 _boxDimensions;
Vec3 _spacing;
FILE* _log;
};
PositionGenerator::PositionGenerator( int numMolecules, int numParticlesPerMolecule, double boxSize ) :
_numMolecules(numMolecules),
_seed(0),
_log(NULL),
_bondDistance(0.1),
_numParticlesPerMolecule(numParticlesPerMolecule),
_numParticles(numMolecules*numParticlesPerMolecule),
_boxSize(boxSize),
_boxDimensions(Vec3(boxSize,boxSize,boxSize)),
_origin(Vec3(0.0,0.0,0.0)) {
double particlesPerDimension = pow( static_cast<double>(_numParticles), (1.0/3.0) );
int particlesPerDimensionI = static_cast<int>(particlesPerDimension+0.999999);
double spacingPerDimension = _boxSize/particlesPerDimension;
_spacing[0] = spacingPerDimension;
_spacing[1] = spacingPerDimension;
_spacing[2] = spacingPerDimension;
}
PositionGenerator::~PositionGenerator( ){};
void PositionGenerator::setBondDistance( double bondDistance ){
_bondDistance = bondDistance;
}
void PositionGenerator::setLog( FILE* log ){
_log = log;
}
double PositionGenerator::getBondDistance( void ) const {
return _bondDistance;
}
double PositionGenerator::getDistance( int index1, int index2, const std::vector<Vec3>& positions ) const {
Vec3 delta = positions[index2] - positions[index1];
return sqrt( delta.dot( delta ) );
}
double PositionGenerator::getPeriodicDistance( int index1, int index2, const std::vector<Vec3>& positions ) const {
Vec3 delta = positions[index2] - positions[index1];
if( _boxSize > 0.0 ){
delta[0] -= floor(delta[0]/_boxSize+0.5f)*_boxSize;
delta[1] -= floor(delta[1]/_boxSize+0.5f)*_boxSize;
delta[2] -= floor(delta[2]/_boxSize+0.5f)*_boxSize;
}
return sqrt( delta.dot( delta ) );
}
/**
* Get positions
*
* @param method method placement
* @param positions output vector of positions
*
* @return nonzero if error detected; 0 otherwise
*/
int PositionGenerator::setPositions( GenerationMethod method, std::vector<Vec3>& positions ) const {
OpenMM_SFMT::SFMT sfmt;
init_gen_rand( _seed, sfmt);
return setPositions( method, sfmt, positions );
}
/**
* Get settings
*
* @return info string
*/
std::string PositionGenerator::getSettings( void ) const {
std::stringstream msg;
msg << "numMolecules " << _numMolecules << std::endl;
msg << "numParticlesPerMolecule " << _numParticlesPerMolecule << std::endl;
msg << "boxSize " << _boxSize << std::endl;
msg << "spacing " << _spacing[0] << std::endl;
msg << "seed " << _seed << std::endl;
return msg.str();
}
/**
* Get positions
*
* @param method method placement
* @param sfmt input random number generator
* @param positions output vector of positions
*
* @return nonzero if error detected; 0 otherwise
*/
int PositionGenerator::setPositions( GenerationMethod method, OpenMM_SFMT::SFMT& sfmt, std::vector<Vec3>& positions ) const {
int errorFlag = 0;
positions.resize( _numParticles );
if( method == Random ){
for( unsigned int ii = 0; ii < static_cast<unsigned int>(_numParticles); ii += _numParticlesPerMolecule ){
positions[ii] = Vec3(_boxSize*genrand_real2(sfmt), _boxSize*genrand_real2(sfmt), _boxSize*genrand_real2(sfmt));
for( unsigned int jj = 1; jj < static_cast<unsigned int>(_numParticlesPerMolecule); jj++) {
positions[ii+jj] = positions[ii] + Vec3(_bondDistance*genrand_real2(sfmt), _bondDistance*genrand_real2(sfmt), _bondDistance*genrand_real2(sfmt));
}
}
} else if( method == SimpleGrid ){
Vec3 origin, boxDimensions, spacing;
std::stringstream msg;
if( _numParticlesPerMolecule > 1 && _bondDistance > 0.0 ){
origin = Vec3(_bondDistance, _bondDistance, _bondDistance );
double particlesPerDimension = pow( static_cast<double>(_numMolecules), (1.0/3.0) );
int particlesPerDimensionI = static_cast<int>(particlesPerDimension+0.999999);
double boxSize = _boxSize;
double spacingPerDimension = (boxSize-_bondDistance)/(particlesPerDimension+1.0);
spacing = Vec3(spacingPerDimension, spacingPerDimension, spacingPerDimension );
boxDimensions = Vec3(boxSize, boxSize, boxSize );
msg << "Bond distance " << _bondDistance << std::endl;
msg << "particlesPerDimension " << particlesPerDimension << std::endl;
msg << "boxSize " << boxSize << std::endl;
msg << "spacingPerDimension " << spacingPerDimension << std::endl;
} else {
origin = _origin;
spacing = _spacing;
boxDimensions = _boxDimensions;
}
msg << getSettings() << std::endl;
if( _log ){
(void) fprintf( _log, "SimpleGrid %s\n", msg.str().c_str() );
}
errorFlag = setParticlesOnGrid( origin, boxDimensions, spacing, sfmt, _bondDistance, positions );
}
return errorFlag;
}
/**
* Place particles on a grid
*
* @param origin origin
* @param boxDimensions box dimensions
* @param spacing spacing
* @param array output vector of grid values
*
* @return -1 if particles will not fit on grid; 0 if they do
*/
int PositionGenerator::setParticlesOnGrid( const Vec3& origin, const Vec3& boxDimensions, const Vec3& spacing, OpenMM_SFMT::SFMT& sfmt,
double bondDistance, std::vector<Vec3>& array ) const {
const double pi = 3.14159265358979323846;
const double pi2 = 2.0*pi;
Vec3 start(origin);
if( array.size() != _numParticles ){
std::stringstream msg;
msg << "PositionGenerator::setParticlesOnGrid position vector size=" << array.size() << " != numParticles=" << _numParticles;
msg << getSettings();
throw OpenMMException( msg.str() );
}
// place molecule centers on grid
for( unsigned int ii = 0; ii < static_cast<unsigned int>(_numParticles); ii += _numParticlesPerMolecule ){
array[ii] = Vec3(start);
bool done = false;
for( unsigned int jj = 0; jj < 3 && !done; jj++ ){
start[jj] += spacing[jj];
if( (start[jj]+4.0*bondDistance) > boxDimensions[jj] ){
start[jj] = origin[jj];
} else {
done = true;
}
}
if( !done ){
std::stringstream msg;
msg << "PositionGenerator::setParticlesOnGrid error in grid settings at molecule=" << ii;
throw OpenMMException( msg.str() );
}
}
// add molecule atoms
for( unsigned int ii = 0; ii < static_cast<unsigned int>(_numMolecules); ii++ ){
int molecularIndex = ii*_numParticlesPerMolecule;
for( unsigned int jj = 1; jj < static_cast<unsigned int>(_numParticlesPerMolecule); jj++ ){
double theta = genrand_real2(sfmt)*pi2;
double phi = genrand_real2(sfmt)*pi;
array[molecularIndex+jj] = array[molecularIndex] + Vec3(_bondDistance*cos(theta)*cos(phi), _bondDistance*cos(theta)*sin(phi), _bondDistance*sin(theta) );
}
}
return 0;
}
/**
* Get enclosing box
*
* @param positions input vector of positions
* @param enclosingBox output Vec3[2] of minimum enclosing box ranges
*
*/
void PositionGenerator::getEnclosingBox( const std::vector<Vec3>& positions, Vec3 enclosingBox[2] ) const {
enclosingBox[0][0] = enclosingBox[1][0] = positions[0][0];
enclosingBox[0][1] = enclosingBox[1][1] = positions[0][1];
enclosingBox[0][2] = enclosingBox[1][2] = positions[0][2];
for( unsigned int ii = 1; ii < positions.size(); ii++ ){
if( enclosingBox[0][0] > positions[ii][0] ){
enclosingBox[0][0] = positions[ii][0];
}
if( enclosingBox[1][0] < positions[ii][0] ){
enclosingBox[1][0] = positions[ii][0];
}
if( enclosingBox[0][1] > positions[ii][1] ){
enclosingBox[0][1] = positions[ii][1];
}
if( enclosingBox[1][1] < positions[ii][1] ){
enclosingBox[1][1] = positions[ii][1];
}
if( enclosingBox[0][2] > positions[ii][2] ){
enclosingBox[0][2] = positions[ii][2];
}
if( enclosingBox[1][2] < positions[ii][2] ){
enclosingBox[1][2] = positions[ii][2];
}
}
return;
}
/**
* Show min/max distances between positions
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param showIndex number of min/maxentries to show
* @param positionIndexVector list of entries to show min/max distances from
*
*/
void PositionGenerator::showMinMaxDistances( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, int showIndex,
const IntVector& positionIndexVector ){
if( !_log )return;
Vec3 box[2];
getEnclosingBox( positions, box );
(void) fprintf( _log, "Enclosing Box (in A): [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e %15.7e]\n",
box[0][0], box[1][0], box[0][1], box[1][1], box[0][2], box[1][2],
(box[1][0] - box[0][0]), (box[1][1] - box[0][1]), (box[1][2] - box[0][2]) );
for( unsigned int ii = 0; ii < positionIndexVector.size(); ii++ ){
if( positionIndexVector[ii] < static_cast<int>(positions.size()) ){
int positionIndex = positionIndexVector[ii];
IntDoublePairVector sortVector;
getSortedDistances( periodicBoundaryConditions, positionIndex, positions, sortVector );
(void) fprintf( _log, "Min/max distance from %6d:\n ", positionIndex );
for( unsigned int jj = 0; jj < sortVector.size() && jj < static_cast<unsigned int>(showIndex); jj++ ){
IntDoublePair pair = sortVector[jj];
(void) fprintf( _log, "[%6d %15.7e] ", pair.first, pair.second);
}
(void) fprintf( _log, "\n " );
for( unsigned int jj = (sortVector.size() - showIndex); jj < sortVector.size() && jj >= 0; jj++ ){
IntDoublePair pair = sortVector[jj];
(void) fprintf( _log, "[%6d %15.7e] ", pair.first, pair.second);
}
(void) fprintf( _log, "\n" );
}
}
return;
}
/**
* Show min/max distances between positions
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param showIndex number of min/maxentries to show
*
*/
void PositionGenerator::showMinMaxDistances( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, int showIndex ){
if( !_log )return;
Vec3 box[2];
getEnclosingBox( positions, box );
(void) fprintf( _log, "Enclosing Box (in A): [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e] [%15.7e %15.7e %15.7e]\n",
box[0][0], box[1][0], box[0][1], box[1][1], box[0][2], box[1][2],
(box[1][0] - box[0][0]), (box[1][1] - box[0][1]), (box[1][2] - box[0][2]) );
IntDoublePairVector hitVector;
double minDistance = 1.0e+30;
double minDistanceCutoff = minDistance*1.1;
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
for( unsigned int jj = ii+1; jj < positions.size(); jj++ ){
double distance = periodicBoundaryConditions ? getPeriodicDistance( jj, ii, positions) :
getDistance( jj, ii, positions);
if( distance < minDistanceCutoff ){
if( distance < minDistance ){
minDistance = distance;
minDistanceCutoff = minDistance*1.1;
}
hitVector.push_back( IntDoublePair(ii*positions.size()+jj,distance ) );
}
}
}
std::sort( hitVector.begin(), hitVector.end(), TestIntDoublePair );
(void) fprintf( _log, "Min distances pbc=%d\n", periodicBoundaryConditions );
for( unsigned int jj = 0; jj < hitVector.size() && jj < static_cast<unsigned int>(showIndex); jj++ ){
IntDoublePair pair = hitVector[jj];
int index = pair.first;
int iIndex = static_cast<int>(index/positions.size());
int jIndex = index - iIndex*positions.size();
(void) fprintf( _log, " [%6d %6d %15.7e]\n", iIndex, jIndex, pair.second);
}
return;
}
/**
* Show particles within a specified distance of a given particle
*
* @param positions input vector of positions
* @param periodicBoundaryConditions if set, use PBC in calculating distances
* @param particleIndex particle to check
* @param distanceToCheckFor distance to check for
* @param tolerance distance tolerance
*
*/
void PositionGenerator::showParticlesWithinDistance( const std::vector<Vec3>& positions,
int periodicBoundaryConditions, unsigned int particleIndex,
double distanceToCheckFor, double tolerance){
if( !_log || particleIndex >= positions.size() )return;
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
double distance = periodicBoundaryConditions ? getPeriodicDistance( particleIndex, ii, positions) :
getDistance( particleIndex, ii, positions);
double delta = fabs( distanceToCheckFor - distance );
if( ii != particleIndex && delta < tolerance ){
(void) fprintf( _log, "Distance=%15.7e between particles %u %u.\n", distance, particleIndex, ii);
}
}
return;
}
/**
* Show distances between positions
*
* @param pairs particle indcies for which distance is to be reported
* @param positions input vector of positions
*
*/
void PositionGenerator::showDistances( const IntIntPairVector& pairs, const std::vector<Vec3>& positions ) const {
for( IntIntPairVectorCI ii = pairs.begin(); ii != pairs.end(); ii++ ){
if( ii->first < static_cast<int>(positions.size()) && ii->second < static_cast<int>(positions.size()) ){
double d = getDistance( ii->first, ii->second, positions );
(void) fprintf( _log, "Distance %6d %6d %15.7e d2=%15.7e\n", ii->first, ii->second, d, d*d );
}
}
return;
}
/**
* Get sorted distances from particular position
*
* @param periodicBoundaryConditions if set, apply PBC
* @param positionIndex input position index
* @param positions input vector of positions
* @param sortVector on output sorted IntDoublePairVector
*
*/
void PositionGenerator::getSortedDistances( int periodicBoundaryConditions, int positionIndex, const std::vector<Vec3>& positions,
IntDoublePairVector& sortVector ) const {
sortVector.resize( 0 );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
if( ii == positionIndex )continue;
double distance = periodicBoundaryConditions ? getPeriodicDistance( positionIndex, ii, positions) : getDistance( positionIndex, ii, positions);
sortVector.push_back( IntDoublePair( ii, distance ) );
}
std::sort( sortVector.begin(), sortVector.end(), TestIntDoublePair );
return;
}
/**---------------------------------------------------------------------------------------
*
* Set string field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setStringFromMap( MapStringString& argumentMap, std::string fieldToCheck, std::string& fieldToSet ){
MapStringStringCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = (*check).second;
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
*
* Set int field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setIntFromMap( MapStringString& argumentMap, std::string fieldToCheck, int& fieldToSet ){
MapStringStringCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = atoi( (*check).second.c_str() );
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
*
* Set int field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setIntFromMapStringToDouble( MapStringToDouble& argumentMap, std::string fieldToCheck, int& fieldToSet ){
// ---------------------------------------------------------------------------------------
MapStringToDoubleCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = static_cast<int>(check->second+0.0000001);
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
* Set float field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setFloatFromMap( MapStringString& argumentMap, std::string fieldToCheck, float& fieldToSet ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "setFloatFromMap";
// ---------------------------------------------------------------------------------------
MapStringStringCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = static_cast<float>(atof( (*check).second.c_str() ));
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
*
* Set double field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setDoubleFromMap( MapStringString& argumentMap, std::string fieldToCheck, double& fieldToSet ){
// ---------------------------------------------------------------------------------------
MapStringStringCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = atof( (*check).second.c_str() );
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
*
* Set double field if in map
*
* @param argumentMap map to check
* @param fieldToCheck key
* @param fieldToSet field to set
*
* @return 1 if argument set, else 0
*
--------------------------------------------------------------------------------------- */
static int setDoubleFromMapStringToDouble( MapStringToDouble& argumentMap, std::string fieldToCheck, double& fieldToSet ){
// ---------------------------------------------------------------------------------------
MapStringToDoubleCI check = argumentMap.find( fieldToCheck );
if( check != argumentMap.end() ){
fieldToSet = check->second;
return 1;
}
return 0;
}
/**---------------------------------------------------------------------------------------
*
* Get relative difference between two forces
*
* @param f1 force1
* @param f2 force2
* @param forceNorm1 output norm of force1
* @param forceNorm2 output norm of force2
* @param relativeDiff output relative difference between force norms
* @param log if set, output forces
*
*
--------------------------------------------------------------------------------------- */
static void getForceRelativeDifference( const Vec3& f1, const Vec3& f2, double& forceNorm1, double& forceNorm2,
double& relativeDiff, FILE* log ) {
double diff = (f1[0] - f2[0])*(f1[0] - f2[0]) +
(f1[1] - f2[1])*(f1[1] - f2[1]) +
(f1[2] - f2[2])*(f1[2] - f2[2]);
forceNorm1 = sqrt( f1[0]*f1[0] + f1[1]*f1[1] + f1[2]*f1[2] );
forceNorm2 = sqrt( f2[0]*f2[0] + f2[1]*f2[1] + f2[2]*f2[2] );
if( forceNorm1 > 0.0 || forceNorm2 > 0.0 ){
relativeDiff = 2.0*sqrt( diff )/(forceNorm1+forceNorm2);
} else {
relativeDiff = 0.0;
}
return;
}
/**---------------------------------------------------------------------------------------
*
* Compare forces from two states
*
* @param state1 state1
* @param state2 state2
* @param relativeTolerance relative tolerance
* @param log if set, output forces
*
* @return number of entries with relative difference > tolerance
*
--------------------------------------------------------------------------------------- */
int compareForcesOfTwoStates( State& state1, State& state2, double relativeTolerance,
DoubleVector& stats, FILE* log ) {
int error = 0;
vector<Vec3> force1 = state1.getForces();
vector<Vec3> force2 = state2.getForces();
double averageRelativeDifference = 0.0;
double count = 0.0;
DoubleVector medians1( force1.size() );
DoubleVector medians2( force1.size() );
IntDoublePairVector relativeDifferences;
for( unsigned int ii = 0; ii < force1.size(); ii++ ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
getForceRelativeDifference( force1[ii], force2[ii], forceNorm1, forceNorm2, relativeDiff, log );
medians1[ii] = forceNorm1;
medians2[ii] = forceNorm2;
relativeDifferences.push_back( IntDoublePair(ii, relativeDiff ) );
averageRelativeDifference += relativeDiff;
count += 1.0;
if( relativeDiff > relativeTolerance ){
error++;
}
if( log ){
(void) fprintf( log, "F %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n", static_cast<unsigned int>(ii),
relativeDiff, force1[ii][0], force1[ii][1], force1[ii][2], force2[ii][0], force2[ii][1], force2[ii][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
// sort relative differences
std::sort( relativeDifferences.begin(), relativeDifferences.end(), TestIntDoublePair );
if( log ){
(void) fprintf( log, "\nEntries w/ largest relative differences.\n" );
for( unsigned int ii = relativeDifferences.size()-1; ii >= relativeDifferences.size()-10 && ii >= 0; ii-- ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
int index = relativeDifferences[ii].first;
getForceRelativeDifference( force1[index], force2[index], forceNorm1, forceNorm2, relativeDiff, log );
(void) fprintf( log, "Fs %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n",
static_cast<unsigned int>(index), relativeDiff,
force1[index][0], force1[index][1], force1[index][2],
force2[index][0], force2[index][1], force2[index][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
if( count > 0.0 ){
averageRelativeDifference /= count;
}
std::sort( medians1.begin(), medians1.end() );
std::sort( medians2.begin(), medians2.end() );
double median1 = medians1[medians1.size()/2];
double median2 = medians2[medians2.size()/2];
stats.resize( 4 );
stats[0] = averageRelativeDifference;
IntDoublePair pair = relativeDifferences[relativeDifferences.size()-1];
stats[1] = pair.second;
stats[2] = static_cast<double>(pair.first);
stats[3] = median1 < median2 ? median1 : median2;
return error;
}
/**
* Create nonbonded force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param reactionFieldDielectric reaction field dielectric
* @param parameterList list of parameters -- used via addParticle()
* @param bonds list of BondInfo_OpenMMTest containing info for exceptions
* @param log logging file (optional -- may be NULL)
*
*/
static NonbondedForce* getNonbondedForce( int nonbondedMethod, double cutoffDistance, double reactionFieldDielectric,
VectorOfDoubleVectors& parameterList, std::vector< BondInfo_OpenMMTest >& bonds, FILE* log ){
NonbondedForce* nonbondedForce = new NonbondedForce();
NonbondedForce::NonbondedMethod method;
switch( nonbondedMethod ){
case NoCutoff_OpenMMTest:
method = NonbondedForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = NonbondedForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = NonbondedForce::CutoffPeriodic;
break;
case Ewald_OpenMMTest:
method = NonbondedForce::Ewald;
break;
case PME_OpenMMTest:
method = NonbondedForce::PME;
break;
default:
method = NonbondedForce::NoCutoff;
}
nonbondedForce->setNonbondedMethod( method );
nonbondedForce->setCutoffDistance( cutoffDistance );
nonbondedForce->setReactionFieldDielectric( reactionFieldDielectric );
// load parameters
for( unsigned int ii = 0; ii < parameterList.size(); ii++ ){
DoubleVector parameters = parameterList[ii];
nonbondedForce->addParticle( parameters[ChargeIndex_OpenMMTest], parameters[SigmaIndex_OpenMMTest], parameters[EpsIndex_OpenMMTest] );
}
// add exceptions
for( unsigned int ii = 0; ii < bonds.size(); ii++ ){
BondInfo_OpenMMTest bond = bonds[ii];
nonbondedForce->addException( bond._particle1, bond._particle2, 0.0f, 1.0, 0.0f );
}
return nonbondedForce;
}
/**
* Create GBVI force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param useQuinticSpline if set use quintic spline for Born radii scaling; else use no scaling
* @param quinticLowerLimitFactor quintic lower limit factor
* @param quinticUpperBornRadiusLimit quintic upper Born radius limit
* @param solventDielectric solvent dielectric
* @param soluteDielectric solute dielectric
* @param parameterList list of parameters -- used via addParticle()
* @param bonds list of BondInfo_OpenMMTest containing info for exceptions
* @param log logging file (optional -- may be NULL)
*
*/
static GBVIForce* getGBVIForce( int nonbondedMethod, double cutoffDistance, int useQuinticSpline,
double quinticLowerLimitFactor, double quinticUpperBornRadiusLimit,
double solventDielectric, double soluteDiecletric,
VectorOfDoubleVectors& parameterList, std::vector< BondInfo_OpenMMTest >& bonds, FILE* log ){
GBVIForce* gbviForce = new GBVIForce();
GBVIForce::NonbondedMethod method;
switch( nonbondedMethod ){
case NoCutoff_OpenMMTest:
method = GBVIForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = GBVIForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = GBVIForce::CutoffPeriodic;
break;
default:
method = GBVIForce::NoCutoff;
}
gbviForce->setNonbondedMethod( method );
gbviForce->setCutoffDistance( cutoffDistance );
gbviForce->setSolventDielectric( solventDielectric );
gbviForce->setSoluteDielectric( soluteDiecletric );
if( useQuinticSpline ){
gbviForce->setBornRadiusScalingMethod( GBVIForce::QuinticSpline );
gbviForce->setQuinticLowerLimitFactor( quinticLowerLimitFactor );
gbviForce->setQuinticUpperBornRadiusLimit( quinticUpperBornRadiusLimit );
} else {
gbviForce->setBornRadiusScalingMethod( GBVIForce::NoScaling );
}
// load parameters
for( unsigned int ii = 0; ii < parameterList.size(); ii++ ){
DoubleVector parameters = parameterList[ii];
gbviForce->addParticle( parameters[ChargeIndex_OpenMMTest], parameters[SigmaIndex_OpenMMTest], parameters[GammaIndex_OpenMMTest] );
}
// add exceptions
for( unsigned int ii = 0; ii < bonds.size(); ii++ ){
BondInfo_OpenMMTest bond = bonds[ii];
gbviForce->addBond( bond._particle1, bond._particle2, bond._distance);
}
return gbviForce;
}
/**
* Create OBC force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param solventDielectric solvent dielectric
* @param soluteDielectric solute dielectric
* @param parameterList list of parameters -- used via addParticle()
* @param log logging file (optional -- may be NULL)
*
*/
static GBSAOBCForce* getGBSAOBCForce( int nonbondedMethod, double cutoffDistance, double solventDielectric, double soluteDiecletric,
VectorOfDoubleVectors& parameterList, FILE* log ){
GBSAOBCForce* obcForce = new GBSAOBCForce();
GBSAOBCForce::NonbondedMethod method;
switch( nonbondedMethod ){
case NoCutoff_OpenMMTest:
method = GBSAOBCForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = GBSAOBCForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = GBSAOBCForce::CutoffPeriodic;
break;
default:
method = GBSAOBCForce::NoCutoff;
}
obcForce->setNonbondedMethod( method );
obcForce->setCutoffDistance( cutoffDistance );
obcForce->setSolventDielectric( solventDielectric );
obcForce->setSoluteDielectric( soluteDiecletric );
// load parameters
for( unsigned int ii = 0; ii < parameterList.size(); ii++ ){
DoubleVector parameters = parameterList[ii];
obcForce->addParticle( parameters[ChargeIndex_OpenMMTest], parameters[SigmaIndex_OpenMMTest], parameters[GammaIndex_OpenMMTest] );
}
return obcForce;
}
/**
* Create nonbonded softcore force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param reactionFieldDielectric reaction field dielectric
* @param log logging file (optional -- may be NULL)
*
*/
#ifdef USE_SOFTCORE
static NonbondedSoftcoreForce* getNonbondedSoftcoreForce( int nonbondedMethod, double cutoffDistance, double reactionFieldDielectric, FILE* log ){
NonbondedSoftcoreForce* nonbondedForce = new NonbondedSoftcoreForce();
NonbondedSoftcoreForce::NonbondedMethod method;
switch( nonbondedMethod ){
case NoCutoff_OpenMMTest:
method = NonbondedSoftcoreForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = NonbondedSoftcoreForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = NonbondedSoftcoreForce::CutoffPeriodic;
break;
default:
method = NonbondedSoftcoreForce::NoCutoff;
}
nonbondedForce->setNonbondedMethod( method );
nonbondedForce->setCutoffDistance( cutoffDistance );
nonbondedForce->setReactionFieldDielectric( reactionFieldDielectric );
return nonbondedForce;
}
/**
* Create GBVI softcore force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param useQuinticSpline if set use quintic spline for Born radii scaling; else use no scaling
* @param quinticLowerLimitFactor quintic lower limit factor
* @param quinticUpperBornRadiusLimit quintic upper Born radius limit
* @param solventDielectric solvent dielectric
* @param soluteDielectric solute dielectric
* @param log logging file (optional -- may be NULL)
*
*/
static GBVISoftcoreForce* getGBVISoftcoreForce( int nonbondedMethod, double cutoffDistance, int useQuinticSpline,
double quinticLowerLimitFactor, double quinticUpperBornRadiusLimit,
double solventDielectric, double soluteDiecletric, FILE* log ){
GBVISoftcoreForce* gbviForce = new GBVISoftcoreForce();
GBVISoftcoreForce::NonbondedMethod method;
switch( gbviMethod ){
case NoCutoff_OpenMMTest:
method = GBVISoftcoreForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = GBVISoftcoreForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = GBVISoftcoreForce::CutoffPeriodic;
break;
default:
method = GBVISoftcoreForce::NoCutoff;
}
gbviForce->setNonbondedMethod( method );
gbviForce->setCutoffDistance( cutoffDistance );
gbviForce->setSolventDielectric( solventDielectric );
gbviForce->setSoluteDielectric( soluteDiecletric );
if( useQuinticSpline ){
gbviForce->setBornRadiusScalingMethod( GBVISoftcoreForce::QuinticSpline );
gbviForce->setQuinticLowerLimitFactor( quinticLowerLimitFactor );
gbviForce->setQuinticUpperBornRadiusLimit( quinticUpperBornRadiusLimit );
} else {
gbviForce->setBornRadiusScalingMethod( GBVISoftcoreForce::NoScaling );
}
return gbviForce;
}
/**
* Create OBC softcore force and set some parameters
*
* @param nonbondedMethod nonbonded method
* @param cutoffDistance cutoff distance
* @param solventDielectric solvent dielectric
* @param soluteDielectric solute dielectric
* @param log logging file (optional -- may be NULL)
*
*/
static GBSAOBCSoftcoreForce* getGBSAOBCSoftcoreForce( int nonbondedMethod, double cutoffDistance,
double solventDielectric, double soluteDiecletric, FILE* log ){
GBSAOBCSoftcoreForce* obcForce = new GBSAOBCSoftcoreForce();
GBSAOBCSoftcoreForce::NonbondedMethod method;
switch( gbviMethod ){
case NoCutoff_OpenMMTest:
method = GBSAOBCSoftcoreForce::NoCutoff;
break;
case CutoffNonPeriodic_OpenMMTest:
method = GBSAOBCSoftcoreForce::CutoffNonPeriodic;
break;
case CutoffPeriodic_OpenMMTest:
method = GBSAOBCSoftcoreForce::CutoffPeriodic;
break;
default:
method = GBSAOBCSoftcoreForce::NoCutoff;
}
obcForce->setNonbondedMethod( method );
obcForce->setCutoffDistance( cutoffDistance );
obcForce->setSolventDielectric( solventDielectric );
obcForce->setSoluteDielectric( soluteDiecletric );
return obcForce;
}
#endif
/**
* Get forces in system
*
* @param system system to serialize
* @param stringForceVector output stringForceVector[forceName] = force index
* @param log logging file (optional -- may be NULL)
*
*/
static void getStringForceMap( System& system, MapStringInt& stringForceVector, FILE* log ){
// print active forces and relevant parameters
for( int ii = 0; ii < system.getNumForces(); ii++ ) {
int hit = 0;
Force& force = system.getForce(ii);
if( !hit ){
try {
CMAPTorsionForce& castForce = dynamic_cast<CMAPTorsionForce&>(force);
stringForceVector["CMAPTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomAngleForce& castForce = dynamic_cast<CustomAngleForce&>(force);
stringForceVector["CustomAngle"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomBondForce& castForce = dynamic_cast<CustomBondForce&>(force);
stringForceVector["CustomBond"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomExternalForce& castForce = dynamic_cast<CustomExternalForce&>(force);
stringForceVector["CustomExternal"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomGBForce& castForce = dynamic_cast<CustomGBForce&>(force);
stringForceVector["CustomGB"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomHbondForce& castForce = dynamic_cast<CustomHbondForce&>(force);
stringForceVector["CustomHbond"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomNonbondedForce& castForce = dynamic_cast<CustomNonbondedForce&>(force);
stringForceVector["CustomNonbonded"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
CustomTorsionForce& castForce = dynamic_cast<CustomTorsionForce&>(force);
stringForceVector["CustomTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
GBSAOBCForce& castForce = dynamic_cast<GBSAOBCForce&>(force);
stringForceVector["GBSAOBC"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
GBVIForce& castForce = dynamic_cast<GBVIForce&>(force);
stringForceVector["GBVI"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
HarmonicAngleForce& castForce = dynamic_cast<HarmonicAngleForce&>(force);
stringForceVector["HarmonicAngle"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
HarmonicBondForce& castForce = dynamic_cast<HarmonicBondForce&>(force);
stringForceVector["HarmonicBond"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
NonbondedForce& castForce = dynamic_cast<NonbondedForce&>(force);
stringForceVector["Nonbonded"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
PeriodicTorsionForce& castForce = dynamic_cast<PeriodicTorsionForce&>(force);
stringForceVector["PeriodicTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
RBTorsionForce& castForce = dynamic_cast<RBTorsionForce&>(force);
stringForceVector["RBTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
MonteCarloBarostat& castForce = dynamic_cast<MonteCarloBarostat&>(force);
stringForceVector["MonteCarloBarostat"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AndersenThermostat& castForce = dynamic_cast<AndersenThermostat&>(force);
stringForceVector["AndersenThermostat"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
#ifdef USE_SOFTCORE
if( !hit ){
try {
GBSAOBCSoftcoreForce& castForce = dynamic_cast<GBSAOBCSoftcoreForce&>(force);
stringForceVector["GBSAOBCSoftcore"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
GBVISoftcoreForce& castForce = dynamic_cast<GBVISoftcoreForce&>(force);
stringForceVector["GBVISoftcore"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
NonbondedSoftcoreForce& castForce = dynamic_cast<NonbondedSoftcoreForce&>(force);
stringForceVector["NonbondedSoftcore"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
#endif
#ifdef INCLUDE_AMOEBA_FORCES
if( !hit ){
try {
AmoebaHarmonicBondForce& castForce = dynamic_cast<AmoebaHarmonicBondForce&>(force);
stringForceVector["AmoebaHarmonicBond"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaHarmonicAngleForce& castForce = dynamic_cast<AmoebaHarmonicAngleForce&>(force);
stringForceVector["AmoebaHarmonicAngle"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaHarmonicInPlaneAngleForce& castForce = dynamic_cast<AmoebaHarmonicInPlaneAngleForce&>(force);
stringForceVector["AmoebaHarmonicInPlaneAngle"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaMultipoleForce& castForce = dynamic_cast<AmoebaMultipoleForce&>(force);
stringForceVector["AmoebaMultipole"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaOutOfPlaneBendForce& castForce = dynamic_cast<AmoebaOutOfPlaneBendForce&>(force);
stringForceVector["AmoebaOutOfPlaneBend"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaPiTorsionForce& castForce = dynamic_cast<AmoebaPiTorsionForce&>(force);
stringForceVector["AmoebaPiTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaStretchBendForce& castForce = dynamic_cast<AmoebaStretchBendForce&>(force);
stringForceVector["AmoebaStretchBend"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaTorsionForce& castForce = dynamic_cast<AmoebaTorsionForce&>(force);
stringForceVector["AmoebaTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaTorsionTorsionForce& castForce = dynamic_cast<AmoebaTorsionTorsionForce&>(force);
stringForceVector["AmoebaTorsionTorsion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaUreyBradleyForce& castForce = dynamic_cast<AmoebaUreyBradleyForce&>(force);
stringForceVector["AmoebaUreyBradley"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaVdwForce& castForce = dynamic_cast<AmoebaVdwForce&>(force);
stringForceVector["AmoebaVdw"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaWcaDispersionForce& castForce = dynamic_cast<AmoebaWcaDispersionForce&>(force);
stringForceVector["AmoebaWcaDispersion"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaGeneralizedKirkwoodForce& castForce = dynamic_cast<AmoebaGeneralizedKirkwoodForce&>(force);
stringForceVector["AmoebaGeneralizedKirkwood"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
if( !hit ){
try {
AmoebaTorsionTorsionForce& castForce = dynamic_cast<AmoebaTorsionTorsionForce&>(force);
stringForceVector["AmoebaTorsionTorsionForce"] = ii;
hit++;
} catch( std::bad_cast ){
}
}
#endif
// COM
if( !hit ){
try {
CMMotionRemover& cMMotionRemover = dynamic_cast<CMMotionRemover&>(force);
hit++;
} catch( std::bad_cast ){
}
}
if( !hit && log ){
(void) fprintf( log, " entry=%2d force not recognized.\n", ii );
}
}
}
/**
* Get forces in system
*
* @param system system to serialize
* @param stringForceVector output stringForceVector[forceName] = force index
* @param log logging file (optional -- may be NULL)
*
*/
static Force* copyForce( const Force& force, FILE* log ){
// print active forces and relevant parameters
Force* forceCopy = NULL;
try {
const CMAPTorsionForce& castForce = dynamic_cast<const CMAPTorsionForce&>(force);
forceCopy = new CMAPTorsionForce( castForce );
} catch( std::bad_cast ){
}
if( forceCopy == NULL ){
try {
const CustomAngleForce& castForce = dynamic_cast<const CustomAngleForce&>(force);
forceCopy = new CustomAngleForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomBondForce& castForce = dynamic_cast<const CustomBondForce&>(force);
forceCopy = new CustomBondForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomExternalForce& castForce = dynamic_cast<const CustomExternalForce&>(force);
forceCopy = new CustomExternalForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomGBForce& castForce = dynamic_cast<const CustomGBForce&>(force);
forceCopy = new CustomGBForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomHbondForce& castForce = dynamic_cast<const CustomHbondForce&>(force);
forceCopy = new CustomHbondForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomNonbondedForce& castForce = dynamic_cast<const CustomNonbondedForce&>(force);
forceCopy = new CustomNonbondedForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const CustomTorsionForce& castForce = dynamic_cast<const CustomTorsionForce&>(force);
forceCopy = new CustomTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const GBSAOBCForce& castForce = dynamic_cast<const GBSAOBCForce&>(force);
forceCopy = new GBSAOBCForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const GBVIForce& castForce = dynamic_cast<const GBVIForce&>(force);
forceCopy = new GBVIForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const HarmonicAngleForce& castForce = dynamic_cast<const HarmonicAngleForce&>(force);
forceCopy = new HarmonicAngleForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const HarmonicBondForce& castForce = dynamic_cast<const HarmonicBondForce&>(force);
forceCopy = new HarmonicBondForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const NonbondedForce& castForce = dynamic_cast<const NonbondedForce&>(force);
forceCopy = new NonbondedForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const PeriodicTorsionForce& castForce = dynamic_cast<const PeriodicTorsionForce&>(force);
forceCopy = new PeriodicTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const RBTorsionForce& castForce = dynamic_cast<const RBTorsionForce&>(force);
forceCopy = new RBTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const MonteCarloBarostat& castForce = dynamic_cast<const MonteCarloBarostat&>(force);
forceCopy = new MonteCarloBarostat( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AndersenThermostat& castForce = dynamic_cast<const AndersenThermostat&>(force);
forceCopy = new AndersenThermostat( castForce );
} catch( std::bad_cast ){
}
}
#ifdef USE_SOFTCORE
if( forceCopy == NULL ){
try {
const GBSAOBCSoftcoreForce& castForce = dynamic_cast<const GBSAOBCSoftcoreForce&>(force);
forceCopy = new GBSAOBCSoftcoreForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const GBVISoftcoreForce& castForce = dynamic_cast<const GBVISoftcoreForce&>(force);
forceCopy = new GBVISoftcoreForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const NonbondedSoftcoreForce& castForce = dynamic_cast<const NonbondedSoftcoreForce&>(force);
forceCopy = new NonbondedSoftcoreForce( castForce );
} catch( std::bad_cast ){
}
}
#endif
#ifdef INCLUDE_AMOEBA_FORCES
if( forceCopy == NULL ){
try {
const AmoebaHarmonicBondForce& castForce = dynamic_cast<const AmoebaHarmonicBondForce&>(force);
forceCopy = new AmoebaHarmonicBondForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaHarmonicAngleForce& castForce = dynamic_cast<const AmoebaHarmonicAngleForce&>(force);
forceCopy = new AmoebaHarmonicAngleForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaHarmonicInPlaneAngleForce& castForce = dynamic_cast<const AmoebaHarmonicInPlaneAngleForce&>(force);
forceCopy = new AmoebaHarmonicInPlaneAngleForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaMultipoleForce& castForce = dynamic_cast<const AmoebaMultipoleForce&>(force);
forceCopy = new AmoebaMultipoleForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaOutOfPlaneBendForce& castForce = dynamic_cast<const AmoebaOutOfPlaneBendForce&>(force);
forceCopy = new AmoebaOutOfPlaneBendForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaPiTorsionForce& castForce = dynamic_cast<const AmoebaPiTorsionForce&>(force);
forceCopy = new AmoebaPiTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaStretchBendForce& castForce = dynamic_cast<const AmoebaStretchBendForce&>(force);
forceCopy = new AmoebaStretchBendForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaTorsionForce& castForce = dynamic_cast<const AmoebaTorsionForce&>(force);
forceCopy = new AmoebaTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaTorsionTorsionForce& castForce = dynamic_cast<const AmoebaTorsionTorsionForce&>(force);
forceCopy = new AmoebaTorsionTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaUreyBradleyForce& castForce = dynamic_cast<const AmoebaUreyBradleyForce&>(force);
forceCopy = new AmoebaUreyBradleyForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaVdwForce& castForce = dynamic_cast<const AmoebaVdwForce&>(force);
forceCopy = new AmoebaVdwForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaWcaDispersionForce& castForce = dynamic_cast<const AmoebaWcaDispersionForce&>(force);
forceCopy = new AmoebaWcaDispersionForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaGeneralizedKirkwoodForce& castForce = dynamic_cast<const AmoebaGeneralizedKirkwoodForce&>(force);
forceCopy = new AmoebaGeneralizedKirkwoodForce( castForce );
} catch( std::bad_cast ){
}
}
if( forceCopy == NULL ){
try {
const AmoebaTorsionTorsionForce& castForce = dynamic_cast<const AmoebaTorsionTorsionForce&>(force);
forceCopy = new AmoebaTorsionTorsionForce( castForce );
} catch( std::bad_cast ){
}
}
#endif
if( log && forceCopy == NULL ){
(void) fprintf( log, " force not recognized.\n" );
}
return forceCopy;
}
/**
* Return copy of system (but not forces)
*
* @param inputSystem system to copy
*
* @return system copy
*
*/
static void copySystem( const System& inputSystem, System& systemCopy, FILE* log ){
// add particle/mass
for( unsigned int ii = 0; ii < static_cast<unsigned int>(inputSystem.getNumParticles()); ii++ ){
systemCopy.addParticle( inputSystem.getParticleMass( static_cast<int>(ii) ) );
}
// box
Vec3 a;
Vec3 b;
Vec3 c;
inputSystem.getDefaultPeriodicBoxVectors( a, b, c );
systemCopy.setDefaultPeriodicBoxVectors( a, b, c );
// copy constraints
for( unsigned int ii = 0; ii < static_cast<unsigned int>(inputSystem.getNumConstraints()); ii++ ){
int particle1, particle2;
double distance;
inputSystem.getConstraintParameters( ii, particle1, particle2, distance);
systemCopy.addConstraint( particle1, particle2, distance);
}
// copy forces
for( unsigned int ii = 0; ii < static_cast<unsigned int>(inputSystem.getNumForces()); ii++ ){
systemCopy.addForce( copyForce( inputSystem.getForce(ii), log) );
}
return;
}
/**
* Randomize parameters
*
* @param parametersLowerBound vector of parameter lower bounds
* @param parametersUpperBound vector of parameter upper bounds
* @param sfmt SFMT random number generator
* @param parameters output vector of randomized parameter values
*
*/
static void randomizeParameters( const std::vector<double>& parametersLowerBound,
const std::vector<double>& parametersUpperBound,
OpenMM_SFMT::SFMT& sfmt, std::vector<double>& parameters ){
if( parametersLowerBound.size() != parametersUpperBound.size() ){
std::stringstream msg;
msg << " randomizeParameters parametersLowerBound size=" << parametersLowerBound.size() << " != parametersUpperBound size=" << parametersUpperBound.size();
throw OpenMMException( msg.str() );
}
if( parametersLowerBound.size() != parameters.size() ){
std::stringstream msg;
msg << " randomizeParameters parametersLowerBound size=" << parametersLowerBound.size() << " != parameter size=" << parameters.size();
throw OpenMMException( msg.str() );
}
for( unsigned int ii = 0; ii < parametersLowerBound.size(); ii++ ){
parameters[ii] = parametersLowerBound[ii] + (parametersUpperBound[ii] - parametersLowerBound[ii])*(genrand_real2(sfmt));
}
return;
}
/**
* Randomize Vec3 vector
*
* @param average mean value
* @param range +/- range
* @param sfmt SFMT random number generator
* @param array output vector of randomized values
*
*/
static void randomizeVec3( double average, double range,
OpenMM_SFMT::SFMT& sfmt, std::vector<Vec3>& array ){
range *= 2.0;
for( unsigned int ii = 0; ii < array.size(); ii++ ){
array[ii] = Vec3( average + range*(genrand_real2(sfmt) - 0.5),
average + range*(genrand_real2(sfmt) - 0.5),
average + range*(genrand_real2(sfmt) - 0.5) );
}
return;
}
/**
* Output contents of MapStringString
*
* @param inputArgumentMap map to output
* @param outputStream output stream
*
*/
static void streamArgumentMap( const MapStringString& inputArgumentMap, std::stringstream& outputStream ){
char buffer[2048];
for( MapStringStringCI ii = inputArgumentMap.begin(); ii != inputArgumentMap.end(); ii++ ){
std::string key = ii->first;
std::string value = ii->second;
(void) sprintf( buffer, " %30s %40s\n", key.c_str(), value.c_str() );
outputStream << buffer;
}
return;
}
/**
* Format argument/value
*
* @param buffer formatted output
* @param key argument name
* @param value argument value
* @param format format string
* @param call if call > 0, skip key name
* @param type type == 0, then use int value; else double
*
*/
static void formatArgument( char* buffer, const std::string& key, double value, const char* format, int call, int type ){
// if call > 0, skip key name
unsigned int index = 0;
while( index < key.size() ){
buffer[index] = call ? ' ' : key[index];
index++;
}
// add blank
buffer[index++] = ' ';
buffer[index] = static_cast<char>(NULL);
if( type == 0 ){
int valueInt = static_cast<int>(value+0.00001);
(void) sprintf( buffer + index, format, valueInt );
} else {
(void) sprintf( buffer + index, format, value );
}
return;
}
/**
* Output contents of MapStringString w/ all argument on one line
*
* @param inputArgumentMap map to output
* @param exclude map of keys to exclude from output
* @param outputStream output stream
*
*/
static void streamArgumentMapOneLine( const MapStringToDouble& inputArgumentMap, const MapStringToInt& exclude,
const StringVector& printFirst, int callId, std::stringstream& outputStream ){
char buffer[2048];
MapStringToInt excludeAll(exclude);
for( unsigned int ii = 0; ii < printFirst.size(); ii++ ){
MapStringToDoubleCI iter = inputArgumentMap.find( printFirst[ii] );
if( iter != inputArgumentMap.end() ){
std::string key = iter->first;
if( exclude.find( key ) == exclude.end() ){
double value = iter->second;
if( key == "numMolecules" ){
formatArgument( buffer, key, value, "%6d ", callId, 0 );
} else if( key == "nonbondedMethod" ){
formatArgument( buffer, key, value, "%1d ", callId, 0 );
} else if( key == "lambda1" || key == "lambda2" ){
formatArgument( buffer, key, value, "%4.2f ", callId, 1 );
} else if( key == "boxSize" ){
formatArgument( buffer, key, value, "%6.2f ", callId, 1 );
} else if( key == "cutoffDistance" ){
formatArgument( buffer, key, value, "%7.3f ", callId, 1 );
} else if( key == "relativeTolerance" ){
formatArgument( buffer, key, value, "%8.1e ", callId, 1 );
} else if( key == "positionPlacementMethod" || key == "applyAssert" || key == "serialize" ){
formatArgument( buffer, key, value, "%1d ", callId, 1 );
} else {
formatArgument( buffer, key, value, "%15.7e ", callId, 1 );
}
outputStream << buffer;
excludeAll[key] = 1;
}
}
}
for( MapStringToDoubleCI ii = inputArgumentMap.begin(); ii != inputArgumentMap.end(); ii++ ){
std::string key = ii->first;
if( excludeAll.find( key ) == excludeAll.end() ){
double value = ii->second;
int valueInt = static_cast<int>(value+0.00001);
double valueDouble = static_cast<double>(valueInt);
if( key == "numMolecules" ){
(void) sprintf( buffer, "%s=%6d ", key.c_str(), valueInt );
} else if( key == "nonbondedMethod" || key == "positionPlacementMethod" || key == "applyAssert" || key == "serialize" ){
(void) sprintf( buffer, "%s=%1d ", key.c_str(), valueInt );
} else if( key == "lambda1" || key == "lambda2" ){
(void) sprintf( buffer, "%s=%4.2f ", key.c_str(), value );
} else if( key == "boxSize" || key == "cutoffDistance" ){
(void) sprintf( buffer, "%s=%6.2f ", key.c_str(), value );
} else if( key == "relativeTolerance" ){
(void) sprintf( buffer, "%s=%8.1e ", key.c_str(), value );
} else if( valueDouble == value ){
(void) sprintf( buffer, "%s=%6d ", key.c_str(), valueInt );
} else {
(void) sprintf( buffer, "%s=%15.7e ", key.c_str(), value );
}
outputStream << buffer;
}
}
outputStream << std::endl;
return;
}
/**
* Get signature of a MapStringToDouble object
*
* @param inputArgumentMap map
* @return signature
*
*/
static double getMapStringToDoubleSignature( const MapStringToDouble& inputArgumentMap ){
double signature = 0.0;
double offset = 0.1;
for( MapStringToDoubleCI ii = inputArgumentMap.begin(); ii != inputArgumentMap.end(); ii++ ){
signature += (offset + ii->second);
offset += 0.1;
}
return signature;
}
/**
* Compare two MapStringToDouble to see if they have the same (key,value) pairs
*
* @param inputArgumentMap1 map 1
* @param inputArgumentMap2 map 2
*
* @return true if maps have same (key,value) pairs; otherwise false
*
*/
static bool compareMapStringToDoubles( const MapStringToDouble& inputArgumentMap1, const MapStringToDouble& inputArgumentMap2 ){
if( inputArgumentMap1.size() != inputArgumentMap1.size() ){
return false;
}
for( MapStringToDoubleCI ii = inputArgumentMap1.begin(); ii != inputArgumentMap1.end(); ii++ ){
MapStringToDoubleCI jj = inputArgumentMap2.find( (*ii).first );
if( jj == inputArgumentMap2.end() || jj->second != ii->second ){
return false;
}
}
return true;
}
/**
* Generate collection of inputArguments maps given
* list of DoubleVectors for each argument
*
* @param inputArguments map[argumentKey] = vector of double parameter values
* @param argumentMaps output vector of generated maps
*
*/
static void generateInputArgumentMapsFromStringVectors( const MapStringToDoubleVector& inputArguments,
VectorOfMapStringToDouble& argumentMaps ){
for( MapStringToDoubleVectorCI ii = inputArguments.begin(); ii != inputArguments.end(); ii++ ){
std::string argumentName = (*ii).first;
DoubleVector arguments = (*ii).second;
unsigned int initialArgumentMapSize = argumentMaps.size();
// generate signature map for each argument map
MapDoubleToInt signatures;
for( unsigned int kk = 0; kk < initialArgumentMapSize; kk++ ){
double signature = getMapStringToDoubleSignature( argumentMaps[kk] );
signatures[signature] = 1;
}
// for each current argumment map, add a new argument map w/ (key,value)
// check that no existing map has the same arguments before adding to the
// vector of argument maps
for( unsigned int kk = 0; kk < initialArgumentMapSize; kk++ ){
for( unsigned int jj = 0; jj < arguments.size(); jj++ ){
MapStringToDouble inputArgumentMap = MapStringToDouble(argumentMaps[kk]);
inputArgumentMap[argumentName] = arguments[jj];
double signature = getMapStringToDoubleSignature( inputArgumentMap );
if( signatures.find( signature ) == signatures.end() ){
argumentMaps.push_back( inputArgumentMap );
} else {
bool match = 0;
for( unsigned int mm = 0; mm < initialArgumentMapSize && !match; mm++ ){
match = compareMapStringToDoubles( inputArgumentMap, argumentMaps[mm] );
}
if( !match ){
argumentMaps.push_back( inputArgumentMap );
}
}
}
}
}
return;
}
/**
* Predicate for sorting map[string] = double
*
* @param d1 first MapStringToDouble to compare
* @param d2 second MapStringToDouble to compare
*
*/
bool TestMapSortPredicate( const MapStringToDouble& d1, const MapStringToDouble& d2 ){
StringVector sortOrder;
sortOrder.push_back( "numMolecules" );
sortOrder.push_back( "nonbondedMethod" );
sortOrder.push_back( "lambda2" );
sortOrder.push_back( "boxSize" );
for( unsigned int ii = 0; ii < sortOrder.size(); ii++ ){
if( d1.find( sortOrder[ii] ) != d1.end() &&
d2.find( sortOrder[ii] ) != d2.end() ){
MapStringToDoubleCI d1i = d1.find( sortOrder[ii] );
MapStringToDoubleCI d2i = d2.find( sortOrder[ii] );
if( d1i->second != d2i->second ){
return d1i->second < d2i->second;
}
}
}
return false;
}
#ifdef USE_SOFTCORE
static CustomNonbondedForce* buildCustomNonbondedSoftcoreForce( const NonbondedSoftcoreForce& nonbondedSoftcoreForce ){
CustomNonbondedForce* customNonbonded;
if( nonbondedSoftcoreForce.getNonbondedMethod() == NoCutoff ){
customNonbonded = new CustomNonbondedForce("lambda*4*eps*(dem^2-dem)+138.935456*q/r;"
"q=q1*q2;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda);"
"sigma=0.5*(sigma1+sigma2);"
"eps=sqrt(eps1*eps2);"
"lambda=min(lambda1,lambda2)");
customNonbonded->setNonbondedMethod( CustomNonbondedForce::NoCutoff );
} else {
customNonbonded = new CustomNonbondedForce("lambda*4*eps*(dem^2-dem)+138.935456*q*(1.0/r+(krf*r*r)-crf);"
"q=q1*q2;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda);"
"sigma=0.5*(sigma1+sigma2);"
"eps=sqrt(eps1*eps2);"
"lambda=min(lambda1,lambda2)");
customNonbonded->setCutoffDistance( nonbondedSoftcoreForce.getCutoffDistance() );
if( nonbondedSoftcoreForce.getNonbondedMethod() == CutoffNonPeriodic ){
customNonbonded->setNonbondedMethod( CustomNonbondedForce::CutoffNonPeriodic );
} else {
customNonbonded->setNonbondedMethod( CustomNonbondedForce::CutoffPeriodic );
}
double cutoffDistance = nonbondedSoftcoreForce.getCutoffDistance();
double reactionFieldDielectric = nonbondedSoftcoreForce.getReactionFieldDielectric();
double eps2 = (reactionFieldDielectric - 1.0)/(2.0*reactionFieldDielectric+1.0);
double kValue = eps2/(cutoffDistance*cutoffDistance*cutoffDistance);
customNonbonded->addGlobalParameter("krf", kValue );
double cValue = (1.0/cutoffDistance)*(3.0*reactionFieldDielectric)/(2.0*reactionFieldDielectric + 1.0);
customNonbonded->addGlobalParameter("crf", cValue );
}
customNonbonded->addPerParticleParameter("q");
customNonbonded->addPerParticleParameter("sigma");
customNonbonded->addPerParticleParameter("eps");
customNonbonded->addPerParticleParameter("lambda");
vector<double> nonbondedParams(4);
for( unsigned int ii = 0; ii < nonbondedSoftcoreForce.getNumParticles(); ii++ ){
double charge;
double sigma;
double epsilon;
double softcoreLJLambda;
nonbondedSoftcoreForce.getParticleParameters(ii, charge, sigma, epsilon, softcoreLJLambda);
nonbondedParams[0] = charge;
nonbondedParams[1] = sigma;
nonbondedParams[2] = epsilon;
nonbondedParams[3] = softcoreLJLambda;
customNonbonded->addParticle( nonbondedParams );
}
return customNonbonded;
}
CustomBondForce* buildCustomBondForceForNonbondedExceptions( const NonbondedSoftcoreForce& nonbondedSoftcoreForce ){
CustomBondForce* customBond;
if( nonbondedSoftcoreForce.getNonbondedMethod() == NoCutoff ){
customBond = new CustomBondForce("lambda*4*eps*(dem^2-dem)+138.935456*q/r;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda)");
} else {
customBond = new CustomBondForce("withinCutoff*(lambda*4*eps*(dem^2-dem)+138.935456*q*(1.0/r+(krf*r*r)-crf));"
"withinCutoff=step(cutoff-r);"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda)");
double cutoffDistance = nonbondedSoftcoreForce.getCutoffDistance();
double reactionFieldDielectric = nonbondedSoftcoreForce.getReactionFieldDielectric();
double eps2 = (reactionFieldDielectric - 1.0)/(2.0*reactionFieldDielectric+1.0);
double kValue = eps2/(cutoffDistance*cutoffDistance*cutoffDistance);
customBond->addGlobalParameter("krf", kValue );
double cValue = (1.0/cutoffDistance)*(3.0*reactionFieldDielectric)/(2.0*reactionFieldDielectric + 1.0);
customBond->addGlobalParameter("crf", cValue );
customBond->addGlobalParameter("cutoff", cutoffDistance );
}
customBond->addPerBondParameter("q");
customBond->addPerBondParameter("sigma");
customBond->addPerBondParameter("eps");
customBond->addPerBondParameter("lambda");
for( unsigned int ii = 0; ii < nonbondedSoftcoreForce.getNumExceptions(); ii++ ){
int particle1, particle2;
double chargeProd;
double sigma;
double epsilon;
double softcoreLJLambda;
nonbondedSoftcoreForce.getExceptionParameters( ii, particle1, particle2, chargeProd, sigma, epsilon, softcoreLJLambda );
vector<double> bondParams(4);
bondParams[0] = chargeProd;
bondParams[1] = sigma;
bondParams[2] = epsilon;
bondParams[3] = softcoreLJLambda;
customBond->addBond( particle1, particle2, bondParams );
}
return customBond;
}
#endif
/**
* Load plugins
*
* @param pluginDirectory plugin directory; if OPENMM_PLUGIN_DIR use ENV variable
* @param log logging file (optional -- may be NULL)
*
*/
static int loadPlugins( const std::string& pluginDirectory, std::vector<std::string>& loaded, FILE* log ){
const char* openmmPluginDirectory;
int envVariableIsSet = 0;
if( pluginDirectory.compare( "OPENMM_PLUGIN_DIR") == 0 ){
openmmPluginDirectory = getenv( "OPENMM_PLUGIN_DIR" );
} else {
openmmPluginDirectory = pluginDirectory.c_str();
}
try {
if( openmmPluginDirectory ){
envVariableIsSet = 1;
if( log ){
(void) fprintf( log, "openmmPluginDirectory=%s\n", openmmPluginDirectory );
(void) fflush( log );
}
loaded = Platform::loadPluginsFromDirectory( openmmPluginDirectory );
if( log ){
(void) fprintf( log, "\nLoaded following %u lib(s) from %s:\n", static_cast<unsigned int>(loaded.size()), openmmPluginDirectory ); (void) fflush( log );
for( unsigned int ii = 0; ii < loaded.size(); ii++ ){
(void) fprintf( log, " %s\n", loaded[ii].c_str() );
}
(void) fprintf( log, "\n" ); (void) fflush( log );
}
} else {
if( log && pluginDirectory == "OPENMM_PLUGIN_DIR" ){
(void) fprintf( log, "Env variable OPENMM_PLUGIN_DIR is not set.\n" );
(void) fflush( log );
}
}
} catch(const exception& e) {
(void) fprintf( log, "Exception: %s\n", e.what() );
(void) fflush( log );
}
return envVariableIsSet;
}
/**
* Set device id
*
* @param platform platform
* @param deviceId device id
* @param log logging file (optional -- may be NULL)
*
*/
static void setDeviceId( Platform& platform, int deviceId, FILE* log ){
std::stringstream deviceIdStr;
deviceIdStr << deviceId;
int wasSet = 0;
if( platform.getName().compare( "Cuda" ) == 0 ){
platform.setPropertyDefaultValue( "CudaDevice", deviceIdStr.str() );
wasSet = 1;
} else if( platform.getName().compare( "OpenCL" ) == 0 ){
platform.setPropertyDefaultValue( "OpenCLDeviceIndex", deviceIdStr.str());
wasSet = 1;
}
if( log && wasSet ){
(void) fprintf( log, "Set deviceId to %d\n", deviceId );
(void) fflush( log );
}
return;
}
/**
* Set device id
*
* @param platform platform
* @param deviceId device id
* @param log logging file (optional -- may be NULL)
*
*/
static void setDeviceIdUsingEnvVariable( Platform& platform, FILE* log ){
const char* deviceId = getenv( "GPU_DEVICE_ID" );
if( deviceId == NULL ){
return;
}
int wasSet = 0;
if( platform.getName().compare( "Cuda" ) == 0 ){
platform.setPropertyDefaultValue( "CudaDevice", deviceId );
wasSet = 1;
} else if( platform.getName().compare( "OpenCL" ) == 0 ){
platform.setPropertyDefaultValue( "OpenCLDeviceIndex", deviceId);
wasSet = 1;
}
if( log && wasSet ){
(void) fprintf( log, "Set deviceId to %s based on env variable GPU_DEVICE_ID setting.\n", deviceId );
(void) fflush( log );
}
return;
}
/**
* Get platform name
*
* @param platformId platformId( 0=Reference, 1=Cuda, 2=OpenCL)
* @param platformName output platform name
*
*/
static void getPlatformName( int platformId, std::string& platformName ){
switch( platformId ){
case Reference_OpenMMTest:
platformName = "Reference";
break;
case Cuda_OpenMMTest:
platformName = "Cuda";
break;
case OpenCL_OpenMMTest:
platformName = "OpenCL";
break;
default:
platformName = "NA";
break;
}
return;
}
/**
* Get lib name
*
* @param libPrefix lib prefix (lib or "")
* @param libSuffix lib suffix (.so, .dylib, .dll)
* @param baseName base name
*
* @return libname
*
*/
static std::string getLibName( const std::string& libPrefix, const std::string& libSuffix, const std::string& baseName ){
std::string fullName = libPrefix;
fullName.append( baseName );
fullName.append( libSuffix );
return fullName;
}
/**
* Get nonbonded method name
*
* @param nonbondedMethod nonbonded method flag
* @return nonbonded method name
*
*/
static std::string getNonbondedMethodName( int nonbondedMethod ){
switch( nonbondedMethod ){
case NoCutoff_OpenMMTest:
return "NoCutoff";
case CutoffNonPeriodic_OpenMMTest:
return "CutoffNonPeriodic";
case CutoffPeriodic_OpenMMTest:
return "CutoffPeriodic";
case Ewald_OpenMMTest:
return "Ewald";
case PME_OpenMMTest:
return "PME";
default:
return "NA";
}
}
/**
* Check if required libs are available
*
* @param requiredLibs list of required libs
* @param loadedLibs list of available libs
* @param log optional logging reference
*
* @return 1 if all required libs are loaded; else 0
*
*/
static int checkRequiredLibsAreAvailable( const StringVector& requiredLibs, const StringVector& loadedLibs, FILE* log ){
unsigned int matchCount = 0;
for( unsigned int kk = 0; kk < requiredLibs.size(); kk++ ){
unsigned int match = 0;
for( unsigned int ii = 0; ii < loadedLibs.size() && match == 0; ii++ ){
if( loadedLibs[ii].compare( requiredLibs[kk] ) == 0 ){
match = 1;
}
}
if( log && !match ){
(void) fprintf( log, "Missing lib %s\n", requiredLibs[kk].c_str() );
}
matchCount += match;
}
int allPresent;
if( matchCount < requiredLibs.size() ){
allPresent = 0;
if( log ){
(void) fprintf( log, "Aborting tests due to missing libs.\n" );
}
} else {
allPresent = 1;
}
return allPresent;
}
/**
* Perform comparison of energies/forces for two systems
*
* @param system1 first system
* @param system2 second system
* @param platform1 first platform name (Reference, Cuda, OpenCL)
* @param platform2 second platform name (Reference, Cuda, OpenCL)
* @param positions positions
* @param inputArgumentMap arguments/flags (relativeTolerance, applyAssert, ...)
* @param idString id string
* @param log logging file (optional -- may be NULL)
*
*/
void runSystemComparisonTest( System& system1, System& system2,
const std::vector<Vec3>& positions, MapStringToDouble& inputArgumentMap,
const std::string& idString, const std::string& precision, FILE* log ){
int applyAssert = 0;
int platformId1 = 0;
int platformId2 = 0;
int deviceId1 = 0;
int deviceId2 = 0;
double relativeTolerance = 1.0e-04;
setDoubleFromMapStringToDouble( inputArgumentMap, "relativeTolerance", relativeTolerance );
setIntFromMapStringToDouble( inputArgumentMap, "applyAssert", applyAssert ) ;
setIntFromMapStringToDouble( inputArgumentMap, "platformId1", platformId1 ) ;
setIntFromMapStringToDouble( inputArgumentMap, "platformId2", platformId2 ) ;
setIntFromMapStringToDouble( inputArgumentMap, "deviceId1", deviceId1 ) ;
setIntFromMapStringToDouble( inputArgumentMap, "deviceId2", deviceId2 ) ;
std::string platformName1;
std::string platformName2;
getPlatformName( platformId1, platformName1 );
getPlatformName( platformId2, platformName2 );
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
if( log ){
(void) fprintf( log, "System1: particles=%d forces=%d System2: particles=%d forces=%d\n",
system1.getNumParticles(), system1.getNumForces(),
system2.getNumParticles(), system2.getNumForces() );
(void) fprintf( log, "Positions=%u\n", static_cast<unsigned int>(positions.size()) );
(void) fprintf( log, "Platform1=%s Platform2=%s\n", platformName1.c_str(), platformName2.c_str() );
(void) fprintf( log, "deviceId1=%d deviceId2=%d\n", deviceId1, deviceId2 );
(void) fprintf( log, "relativeTolerance=%8.2e applyAssert=%d\n", relativeTolerance, applyAssert );
MapStringInt stringForceVector1;
MapStringInt stringForceVector2;
getStringForceMap( system1, stringForceVector1, log );
(void) fprintf( log, "Forces in system 1: [" );
for( MapStringIntCI ii = stringForceVector1.begin(); ii != stringForceVector1.end(); ii++ ){
(void) fprintf( log, " %s ", ii->first.c_str() );
}
getStringForceMap( system2, stringForceVector2, log );
(void) fprintf( log, "]\nForces in system 2: [" );
for( MapStringIntCI ii = stringForceVector2.begin(); ii != stringForceVector2.end(); ii++ ){
(void) fprintf( log, " %s ", ii->first.c_str() );
}
(void) fprintf( log, "]\n" );
}
if( system1.getNumParticles() != system2.getNumParticles() ){
std::stringstream msg;
msg << "Number of particles for systems to be compared are unequal: " << system1.getNumParticles() << " != " << system2.getNumParticles();
throw OpenMMException( msg.str() );
}
if( system1.getNumParticles() != static_cast<int>(positions.size()) ){
std::stringstream msg;
msg << "Number of particles for system does not equal size of position array: " << system1.getNumParticles() << " != " << positions.size();
throw OpenMMException( msg.str() );
}
#if TEST_PLATFORM == TEST_OPENCL_PLATFORM
ReferencePlatform platform1;
OpenCLPlatform platform2;
platform2.setPropertyDefaultValue("OpenCLPrecision", precision);
#elif TEST_PLATFORM == TEST_CUDA_PLATFORM
ReferencePlatform platform1;
CudaPlatform platform2;
platform2.setPropertyDefaultValue("CudaPrecision", precision);
#else
Platform& platform1 = Platform::getPlatformByName( platformName1 );
if( deviceId1 ){
setDeviceId( platform1, deviceId1, log );
}
setDeviceIdUsingEnvVariable( platform1, log );
Platform& platform2 = Platform::getPlatformByName( platformName2 );
if( deviceId2 ){
setDeviceId( platform2, deviceId2, log );
}
setDeviceIdUsingEnvVariable( platform2, log );
#endif
Context context1( system1, integrator1, platform1 );
context1.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
Context context2( system2, integrator2, platform2 );
context2.setPositions(positions);
State state2 = context2.getState(State::Forces | State::Energy);
double energyDiff = 0.0;
if( fabs( state1.getPotentialEnergy() ) > 0.0 || fabs( state2.getPotentialEnergy()) > 0.0 ){
energyDiff = fabs( state1.getPotentialEnergy() - state2.getPotentialEnergy() )/( fabs( state1.getPotentialEnergy() ) + fabs( state2.getPotentialEnergy() ) );
}
if( log ){
DoubleVector stats;
compareForcesOfTwoStates( state1, state2, relativeTolerance, stats, log );
(void) fprintf( log, "%s %6d eDff=%15.7e fMx=%15.7e fAvg=%15.7e fMed=%15.7e eCd=%15.7e eRf=%15.7e mxFIdx=%d\n",
idString.c_str(), system1.getNumParticles(), energyDiff,
stats[1], stats[0], stats[3], state1.getPotentialEnergy(), state2.getPotentialEnergy(), static_cast<int>(stats[2]+0.0001));
(void) fflush( log );
}
if( applyAssert ){
ASSERT( energyDiff < relativeTolerance );
for( int ii = 0; ii < system1.getNumParticles(); ii++ ){
Vec3 f1 = state1.getForces()[ii];
Vec3 f2 = state2.getForces()[ii];
double f1N = sqrt( (f1[0]*f1[0]) + (f1[1]*f1[1]) + (f1[2]*f1[2]) );
double f2N = sqrt( (f2[0]*f2[0]) + (f2[1]*f2[1]) + (f2[2]*f2[2]) );
double diff = (f1[0]-f2[0])*(f1[0]-f2[0]) +
(f1[1]-f2[1])*(f1[1]-f2[1]) +
(f1[2]-f2[2])*(f1[2]-f2[2]);
if( f1N > 0.0 || f1N > 0.0 ){
diff = 2.0*sqrt( diff )/(f1N + f2N);
}
ASSERT( diff < relativeTolerance );
}
}
}
/**
* Serialize system
*
* @param system system to serialize
* @param serializeFileName file name for xml output
* @param log logging file (optional -- may be NULL)
*
*/
void serializeSystem( System& system, const std::string& serializeFileName, FILE* log ){
#ifdef OPENMM_SERIALIZE
//registerAmoebaSerializationProxies();
std::stringstream buffer;
XmlSerializer::serialize<System>(&system, "System", buffer);
FILE* filePtr = fopen( serializeFileName.c_str(), "w" );
if( filePtr == NULL ){
if( log ){
(void) fprintf( log, "Unable to open xml file %s\n", serializeFileName.c_str() );
return;
}
}
(void) fprintf( filePtr, "%s", buffer.str().c_str() );
(void) fclose( filePtr );
if( log ){
(void) fprintf( log, "Wrote system to xml file %s\n", serializeFileName.c_str() );
}
#endif
return;
}
/**
* Output vector of Vec3 to file
*
* @param positions system to serialize
* @param fileName file name for output
* @param log logging file (optional -- may be NULL)
*
*/
void serializeVectorOfVec3( const std::vector<Vec3>& positions, std::string fileName, FILE* log ){
#ifdef OPENMM_SERIALIZE
FILE* filePtr = fopen( fileName.c_str(), "w" );
if( filePtr == NULL ){
if( log ){
(void) fprintf( log, "Unable to open Vec3 file %s\n", fileName.c_str() );
return;
}
}
(void) fprintf( filePtr, "Positions %u\n", static_cast<unsigned int>(positions.size()) );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( filePtr, "%9u %17.10e %17.10e %17.10e\n", ii, positions[ii][0], positions[ii][1], positions[ii][2] );
}
(void) fclose( filePtr );
if( log ){
(void) fprintf( log, "Wrote to file %s\n", fileName.c_str() );
}
#endif
return;
}
/**
* Serialize system and positions
*
* @param system system to serialize
* @param positions positions to output
* @param baseFileName base file name for xml/txt output
* @param log logging file (optional -- may be NULL)
*
*/
void serializeSystemAndPositions( System& system, const std::vector<Vec3>& positions, const std::string& baseFileName, FILE* log ){
std::stringstream xmlfileName;
xmlfileName << baseFileName << ".xml";
serializeSystem( system, xmlfileName.str(), log );
std::stringstream posfileName;
posfileName << baseFileName << ".txt";
serializeVectorOfVec3( positions, posfileName.str(), log );
return;
}
void runTests( MapStringToDouble& inputArgumentMap, const std::string& precision, FILE* log ){
double lambda1 = 1.0;
double lambda2 = 1.0;
int nonbondedMethod = 0;
int numMolecules = 1;
int numParticlesPerMolecule = 2;
int useQuinticSpline = 1;
int applyAssert = 1;
int positionPlacementMethod = 0;
int serialize = 0;
double boxSize = 10.0;
double relativeTolerance = 1.0e-04;
double quinticLowerLimitFactor = 0.8;
double quinticUpperBornRadiusLimit = 2.0;
std::stringstream baseFileName;
setDoubleFromMapStringToDouble( inputArgumentMap, "lambda1", lambda1 );
setDoubleFromMapStringToDouble( inputArgumentMap, "lambda2", lambda2 );
setDoubleFromMapStringToDouble( inputArgumentMap, "boxSize", boxSize );
double cutoffDistance = boxSize*0.4;
setDoubleFromMapStringToDouble( inputArgumentMap, "cutoffDistance", cutoffDistance);
setDoubleFromMapStringToDouble( inputArgumentMap, "relativeTolerance", relativeTolerance );
baseFileName << "Nb";
#if IMPLICIT_SOLVENT == TEST_GBVI
setDoubleFromMapStringToDouble( inputArgumentMap, "quinticLowerLimitFactor", quinticLowerLimitFactor );
setDoubleFromMapStringToDouble( inputArgumentMap, "quinticUpperBornRadiusLimit", quinticUpperBornRadiusLimit );
baseFileName << "Gbvi";
#endif
#if IMPLICIT_SOLVENT == TEST_OBC
baseFileName << "Obc";
#endif
setIntFromMapStringToDouble( inputArgumentMap, "positionPlacementMethod", positionPlacementMethod ) ;
setIntFromMapStringToDouble( inputArgumentMap, "nonbondedMethod", nonbondedMethod );
setIntFromMapStringToDouble( inputArgumentMap, "numMolecules", numMolecules );
setIntFromMapStringToDouble( inputArgumentMap, "numParticlesPerMolecule", numParticlesPerMolecule );
setIntFromMapStringToDouble( inputArgumentMap, "serialize", serialize );
setIntFromMapStringToDouble( inputArgumentMap, "applyAssert", applyAssert );
double bondDistance = 0.1;
double minDistance = 0.1;
double cellSize = 2.0*bondDistance + minDistance;
double boxLength = cellSize*pow( static_cast<double>(numMolecules), 0.333333 );
if( positionPlacementMethod == 1 && boxLength > boxSize ){
boxSize = boxLength;
if( log ){
// (void) fprintf( log, "Updated box size: bL=%6.3f cell=%6.2e bond=%5.2f separation=%5.2f\n", boxLength, cellSize, bondDistance, minDistance );
}
}
if( nonbondedMethod >= 2 && cutoffDistance > boxSize*0.5 ){
cutoffDistance = boxSize*0.49;
}
int numParticles = numMolecules*numParticlesPerMolecule;
int includeGbvi = 1;
double reactionFieldDielectric = 80.0;
if( log ){
double particleDensity = static_cast<double>(numParticles)/(boxSize*boxSize*boxSize);
double particleCube = pow( particleDensity, (-1.0/3.0) );
(void) fprintf( log, "\n--------------------------------------------------------------------------------------\n" );
(void) fprintf( log, "Input arguments\n" );
(void) fflush( log );
//(void) fprintf( log, " includeGbvi %d\n", includeGbvi );
(void) fprintf( log, " nonbondedMethod %d\n", nonbondedMethod );
(void) fprintf( log, " numParticles %d\n", numParticles );
(void) fprintf( log, " numMolecules %d\n", numMolecules );
(void) fprintf( log, " numParticlesPerMolecule %d\n", numParticlesPerMolecule );
(void) fprintf( log, " positionPlacementMethod %d\n", positionPlacementMethod);
(void) fprintf( log, " boxSize %8.3f\n", boxSize );
(void) fprintf( log, " cutoffDistance %15.7e\n", cutoffDistance );
(void) fprintf( log, " reactionFieldDielectric %8.3f\n", reactionFieldDielectric );
#if IMPLICIT_SOLVENT == TEST_GBVI
(void) fprintf( log, " useQuinticSpline %d\n", useQuinticSpline );
(void) fprintf( log, " quinticLowerLimitFactor %8.3f\n", quinticLowerLimitFactor );
(void) fprintf( log, " quinticUpperBornRadiusLimit %8.3f\n", quinticUpperBornRadiusLimit );
#endif
#ifdef USE_SOFTCORE
(void) fprintf( log, " lambda1 %8.3f\n", lambda1 );
(void) fprintf( log, " lambda2 %8.3f\n", lambda2 );
#endif
(void) fprintf( log, " relativeTolerance %8.1e\n", relativeTolerance );
(void) fprintf( log, " particleDensity %8.2e\n", particleDensity );
(void) fprintf( log, " particleCube %8.2e\n", particleCube );
}
// Create two systems: one with GbviSoftcoreForce NonbondedSoftcoreForce forces, and one using a CustomNonbondedForce, CustomGBVI force to implement the same interaction.
System standardSystem;
for (int i = 0; i < numParticles; i++) {
standardSystem.addParticle(1.0);
}
standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
double solventDielectric = 78.3; // 1.0 or 1.0e+10
double soluteDiecletric = 1.0;
std::vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
PositionGenerator positionGenerator( numMolecules, numParticlesPerMolecule, boxSize );
if( log ){
positionGenerator.setLog( log );
}
if( positionPlacementMethod == 1 ){
positionGenerator.setPositions( PositionGenerator::SimpleGrid, sfmt, positions );
} else {
positionGenerator.setBondDistance( 0.3 );
positionGenerator.setPositions( PositionGenerator::Random, sfmt, positions );
}
// show info on particle positions
if( log ){
int periodicBoundaryConditions = (nonbondedMethod > CutoffNonPeriodic_OpenMMTest) ? 1 : 0;
int showIndex = 5;
double distanceTolerance = 1.0e-04;
IntVector positionIndexVector;
positionIndexVector.push_back( 0 );
positionIndexVector.push_back( 5713 );
positionIndexVector.push_back( 6291 );
positionIndexVector.push_back( 3191 );
positionIndexVector.push_back( 3769 );
positionIndexVector.push_back( static_cast<int>(positions.size())-1 );
positionGenerator.showMinMaxDistances( positions, periodicBoundaryConditions, showIndex, positionIndexVector);
positionGenerator.showMinMaxDistances( positions, periodicBoundaryConditions, showIndex );
positionGenerator.showParticlesWithinDistance( positions, periodicBoundaryConditions, 5713, cutoffDistance, distanceTolerance );
positionGenerator.showParticlesWithinDistance( positions, periodicBoundaryConditions, 6291, cutoffDistance, distanceTolerance );
IntIntPairVector pairs;
pairs.push_back( IntIntPair( 5713, 6291 ) );
pairs.push_back( IntIntPair( 5713, 3191 ) );
pairs.push_back( IntIntPair( 5713, 3769 ) );
pairs.push_back( IntIntPair( 6291, 3191 ) );
positionGenerator.showDistances( pairs, positions );
}
const int numberOfParameters = 5;
std::vector<double> parameterLowerBound( numberOfParameters, 0.0 );
double fixedCharge = 0.1;
parameterLowerBound[ChargeIndex_OpenMMTest] = fixedCharge; // charge
parameterLowerBound[SigmaIndex_OpenMMTest] = 0.1; // sigma
parameterLowerBound[EpsIndex_OpenMMTest] = 0.5; // eps
parameterLowerBound[GammaIndex_OpenMMTest] = 0.1; // gamma
parameterLowerBound[LambdaIndex_OpenMMTest] = lambda1; // lambda
std::vector<double> parameterUpperBound( parameterLowerBound );
parameterUpperBound[ChargeIndex_OpenMMTest] = fixedCharge; // charge
parameterUpperBound[SigmaIndex_OpenMMTest] = 0.3; // sigma
parameterUpperBound[EpsIndex_OpenMMTest] = 40.0; // eps
parameterUpperBound[GammaIndex_OpenMMTest] = 40.0; // gamma
#if IMPLICIT_SOLVENT == TEST_OBC
parameterLowerBound[GammaIndex_OpenMMTest] = 0.1; // overlap factor
parameterUpperBound[GammaIndex_OpenMMTest] = 1.5;
#endif
std::vector<double> parameters( numberOfParameters );
VectorOfDoubleVectors parameterList;
std::vector< BondInfo_OpenMMTest > bonds;
double charge = fixedCharge;
for( int ii = 0; ii < numMolecules; ii++) {
charge *= -1.0;
double lambda = ii < (numMolecules/2) ? lambda1 : lambda2;
randomizeParameters( parameterLowerBound, parameterUpperBound, sfmt, parameters );
parameters[ChargeIndex_OpenMMTest] = charge;
parameters[LambdaIndex_OpenMMTest] = lambda;
parameterList.push_back( parameters );
int baseParticleIndex = ii*numParticlesPerMolecule;
for( int jj = 1; jj < numParticlesPerMolecule; jj++) {
// alternate charges
charge *= -1.0;
randomizeParameters( parameterLowerBound, parameterUpperBound, sfmt, parameters );
parameters[ChargeIndex_OpenMMTest] = charge;
parameters[LambdaIndex_OpenMMTest] = lambda;
parameterList.push_back( parameters );
double bondDistance = positionGenerator.getDistance( baseParticleIndex, baseParticleIndex+jj, positions );
bonds.push_back( BondInfo_OpenMMTest( baseParticleIndex, baseParticleIndex+jj, bondDistance ) );
}
// alternate charge if numParticlesPerMolecule is odd
if( (numParticlesPerMolecule % 2) ){
charge *= -1.0;
}
}
#ifdef USE_SOFTCORE
baseFileName << "Softcore";
NonbondedSoftcoreForce* nonbondedSoftcoreForce = getNonbondedSoftcoreForce( nonbondedMethod, cutoffDistance, reactionFieldDielectric,
parameterList, bonds, log);
#if IMPLICIT_SOLVENT == TEST_GBVI
GBVISoftcoreForce* gbviSoftcoreForce = getGBVISoftcoreForce( nonbondedMethod, nonbondedSoftcoreForce->getCutoffDistance(),
useQuinticSpline, quinticLowerLimitFactor, quinticUpperBornRadiusLimit,
solventDielectric, soluteDiecletric,
parameterList, bonds, log );
#endif
#if IMPLICIT_SOLVENT == TEST_OBC
GBSAOBCSoftcoreForce* gbviSoftcoreForce = getGBSAOBCSoftcoreForce( nonbondedMethod, nonbondedSoftcoreForce->getCutoffDistance(),
solventDielectric, soluteDiecletric,
parameterList, log );
#endif
#else
NonbondedForce* nonbondedSoftcoreForce = getNonbondedForce( nonbondedMethod, cutoffDistance, reactionFieldDielectric,
parameterList, bonds, log);
#if IMPLICIT_SOLVENT == TEST_GBVI
GBVIForce* gbviSoftcoreForce = getGBVIForce( nonbondedMethod, nonbondedSoftcoreForce->getCutoffDistance(),
useQuinticSpline, quinticLowerLimitFactor, quinticUpperBornRadiusLimit,
solventDielectric, soluteDiecletric, parameterList, bonds, log );
#endif
#if IMPLICIT_SOLVENT == TEST_OBC
GBSAOBCForce* gbviSoftcoreForce = getGBSAOBCForce( nonbondedMethod, nonbondedSoftcoreForce->getCutoffDistance(),
solventDielectric, soluteDiecletric, parameterList, log );
#endif
#endif
standardSystem.addForce(nonbondedSoftcoreForce);
#if IMPLICIT_SOLVENT > 0
if( includeGbvi ){
standardSystem.addForce(gbviSoftcoreForce);
}
#endif
// copy system and forces
System systemCopy;
copySystem( standardSystem, systemCopy, log );
// serialize
baseFileName << "_N" << positions.size();
baseFileName << "_Nb" << nonbondedMethod;
serializeSystemAndPositions( standardSystem, positions, baseFileName.str(), log);
// perform comparison
std::stringstream idString;
idString << "Nb " << nonbondedMethod << " l2 " << std::fixed << setprecision(2) << lambda2;
runSystemComparisonTest( standardSystem, systemCopy, positions, inputArgumentMap, idString.str(), precision, log );
}
int main(int argc, char* argv[]) {
std::string precision;
if (argc > 1)
precision = std::string(argv[1]);
else
precision = "single";
try {
#ifdef USE_SOFTCORE
registerFreeEnergyCudaKernelFactories( );
#endif
VectorOfMapStringToDouble vectorOfMapStringToDouble;
MapStringToDouble inputArgumentMap;
MapStringToDoubleVector generativeArgumentMaps;
//FILE* log = stderr;
FILE* log = NULL;
std::vector<std::string> loadedLibs;
int envVariableIsSet = loadPlugins( "OPENMM_PLUGIN_DIR", loadedLibs, log );
inputArgumentMap["platformId1"] = 0;
inputArgumentMap["platformId2"] = 1;
inputArgumentMap["deviceId1"] = 0;
inputArgumentMap["deviceId2"] = 0;
inputArgumentMap["lambda2"] = 1.0;
inputArgumentMap["nonbondedMethod"] = 0;
inputArgumentMap["numMolecules"] = 10;
inputArgumentMap["boxSize"] = 5.0;
inputArgumentMap["positionPlacementMethod"] = 1;
inputArgumentMap["cutoffDistance"] = 0.301*inputArgumentMap["boxSize"];
//inputArgumentMap["cutoffDistance"] = 1.0;
inputArgumentMap["relativeTolerance"] = 5.0e-04;
inputArgumentMap["applyAssert"] = 1;
inputArgumentMap["serialize"] = 1;
inputArgumentMap["numParticlesPerMolecule"] = 2;
#ifdef USE_SOFTCORE
DoubleVector lamda2;
lamda2.push_back( 1.0 );
lamda2.push_back( 0.5 );
lamda2.push_back( 0.0 );
if( lamda2.size() > 0 ){
generativeArgumentMaps["lambda2"] = lamda2;
inputArgumentMap["lambda2"] = lamda2[0];
}
#endif
DoubleVector numberOfMolecules;
//numberOfMolecules.push_back( 10 );
#if IMPLICIT_SOLVENT != TEST_NONBONDED
numberOfMolecules.push_back( 100 );
#endif
numberOfMolecules.push_back( 1000 );
numberOfMolecules.push_back( 2000 );
numberOfMolecules.push_back( 4000 );
//numberOfMolecules.push_back( 8000 );
if( numberOfMolecules.size() > 0 ){
generativeArgumentMaps["numMolecules"] = numberOfMolecules;
inputArgumentMap["numMolecules"] = numberOfMolecules[0];
}
DoubleVector nonbondedMethod;
nonbondedMethod.push_back( NoCutoff_OpenMMTest );
nonbondedMethod.push_back( CutoffNonPeriodic_OpenMMTest );
nonbondedMethod.push_back( CutoffPeriodic_OpenMMTest );
#if IMPLICIT_SOLVENT == TEST_NONBONDED
nonbondedMethod.push_back( Ewald_OpenMMTest );
nonbondedMethod.push_back( PME_OpenMMTest );
#endif
if( nonbondedMethod.size() > 0 ){
generativeArgumentMaps["nonbondedMethod"] = nonbondedMethod;
inputArgumentMap["nonbondedMethod"] = nonbondedMethod[0];
}
DoubleVector platformId2s;
#if TEST_PLATFORM == TEST_OPENCL_PLATFORM
platformId2s.push_back( OpenCL_OpenMMTest );
#elif TEST_PLATFORM == TEST_CUDA_PLATFORM
platformId2s.push_back( Cuda_OpenMMTest );
#else
platformId2s.push_back( Cuda_OpenMMTest );
#endif
// check that required libs are available for platform to be tested
// if unavailable, skip tests
std::string libPrefix = "lib";
std::string libSuffix = ".so";
#ifdef _MSC_VER
libPrefix = "";
libSuffix = ".dll";
#endif
#ifdef __APPLE__
libSuffix = ".dylib";
#endif
StringVector requiredLibs;
for( unsigned int kk = 0; kk < platformId2s.size(); kk++ ){
if( platformId2s[kk] == OpenCL_OpenMMTest ){
requiredLibs.push_back( getLibName( libPrefix, libSuffix, "OpenMMOpenCL" ) );
}
if( platformId2s[kk] == Cuda_OpenMMTest ){
requiredLibs.push_back( getLibName( libPrefix, libSuffix, "OpenMMCuda") );
#ifdef USE_SOFTCORE
requiredLibs.push_back( getLibName( libPrefix, libSuffix, "OpenMMFreeEnergy" ) );
requiredLibs.push_back( getLibName( libPrefix, libSuffix, "OpenMMFreeEnergyCuda" ) );
#endif
}
}
// if TEST_PLATFORM is not set, then check that required libs are available
#if TEST_PLATFORM != TEST_OPENCL_PLATFORM && TEST_PLATFORM != TEST_CUDA_PLATFORM
envVariableIsSet = checkRequiredLibsAreAvailable( requiredLibs, loadedLibs, log );
if( envVariableIsSet == 0 && log ){
(void) fprintf( log, "Aborting tests due to missing libs.\n" );
}
#else
// unit test path: force tests to run
envVariableIsSet = 1;
#endif
if( platformId2s.size() > 0 ){
generativeArgumentMaps["platformId2"] = platformId2s;
inputArgumentMap["platformId2"] = platformId2s[0];
}
vectorOfMapStringToDouble.push_back( inputArgumentMap );
generateInputArgumentMapsFromStringVectors( generativeArgumentMaps, vectorOfMapStringToDouble );
// modify relative tolerance for large systems
// case: Distance 433 669 1.5000001e+00 d2= 2.2500002e+00 w/ cutoff=1.500
for( unsigned int kk = 0; kk < vectorOfMapStringToDouble.size(); kk++ ){
int numMolecules = 0;
int numParticlesPerMolecule = 2;
setIntFromMapStringToDouble( vectorOfMapStringToDouble[kk], "numMolecules", numMolecules );
setIntFromMapStringToDouble( vectorOfMapStringToDouble[kk], "numParticlesPerMolecule", numParticlesPerMolecule );
if( numMolecules*numParticlesPerMolecule > 1000 ){
vectorOfMapStringToDouble[kk]["relativeTolerance"] = 6.0e-03;
}
}
if( log ){
MapStringToInt exclude;
exclude["lambda1"] = 1;
exclude["platformId1"] = 1;
exclude["platformId2"] = 1;
exclude["deviceId1"] = 1;
exclude["deviceId2"] = 1;
exclude["numParticlesPerMolecule"] = 1;
std::stringstream outputStream;
std::sort( vectorOfMapStringToDouble.begin(), vectorOfMapStringToDouble.end(), TestMapSortPredicate);
StringVector printOrder;
printOrder.push_back( "numMolecules" );
printOrder.push_back( "nonbondedMethod" );
printOrder.push_back( "lambda2" );
printOrder.push_back( "boxSize" );
for( unsigned int kk = 0; kk < vectorOfMapStringToDouble.size(); kk++ ){
streamArgumentMapOneLine( vectorOfMapStringToDouble[kk], exclude, printOrder, kk, outputStream );
}
(void) fprintf( log, "Initial argument maps: %u\n%s", static_cast<unsigned int>(vectorOfMapStringToDouble.size()), outputStream.str().c_str() );
}
// run tests
if( envVariableIsSet ){
int wasException = 0;
for( unsigned int kk = 0; kk < vectorOfMapStringToDouble.size() && wasException < 3; kk++ ){
try {
runTests( vectorOfMapStringToDouble[kk], precision, log );
} catch(const exception& e) {
std::stringstream msg;
#if IMPLICIT_SOLVENT == TEST_NONBONDED
msg << "Nonbonded";
#elif IMPLICIT_SOLVENT == TEST_OBC
msg << "GBSAOBC";
#elif IMPLICIT_SOLVENT == TEST_GBVI
msg << "GBVI";
#endif
int numMolecules = 0;
int numParticlesPerMolecule = 0;
int nonbondedMethod = 0;
setIntFromMapStringToDouble( vectorOfMapStringToDouble[kk], "numMolecules", numMolecules );
setIntFromMapStringToDouble( vectorOfMapStringToDouble[kk], "numParticlesPerMolecule", numParticlesPerMolecule );
setIntFromMapStringToDouble( vectorOfMapStringToDouble[kk], "nonbondedMethod", nonbondedMethod);
msg << " test: system size=" << numMolecules*numParticlesPerMolecule << " nonbonded method=" << getNonbondedMethodName( nonbondedMethod );
msg << " exception: " << e.what() << endl;
// msg << "Note cases have been encountered for nonbonded methods with cutoffs where the error was due to particles being within 1.0e-05 of the cutoff." << endl;
cout << msg.str();
wasException += 1;
}
}
if( wasException ){
return 1;
}
}
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/reference/src/ReferenceTabulatedFunction.cpp
View file @ 4233befd
/* -------------------------------------------------------------------------- *
* 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);
}
/* -------------------------------------------------------------------------- *
* 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"
#ifdef _MSC_VER
/**
* We need to define this ourselves, since Visual Studio is missing round() from cmath.
*/
static int round(double x) {
return (int) (x+0.5);
}
#else
#include <cmath>
#endif
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/src/SimTKReference/ReferenceCustomCompoundBondIxn.cpp
View file @ 4233befd
......@@ -124,7 +124,7 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<RealVec>&
// Compute all of the variables the energy can depend on.
const vector<int>& atoms = bondAtoms[0];
const vector<int>& atoms = bondAtoms[bond];
for (int i = 0; i < (int) particleTerms.size(); i++) {
const ParticleTermInfo& term = particleTerms[i];
variables[term.name] = atomCoordinates[term.atom][term.component];
......
platforms/reference/tests/TestReferenceCustomCompoundBondForce.cpp
View file @ 4233befd
......@@ -266,12 +266,61 @@ void testContinuous3DFunction() {
}
}
void testMultipleBonds() {
// Two compound bonds using Urey-Bradley example from API doc
ReferencePlatform platform;
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(3,
"0.5*(kangle*(angle(p1,p2,p3)-theta0)^2+kbond*(distance(p1,p3)-r0)^2)");
custom->addPerBondParameter("kangle");
custom->addPerBondParameter("kbond");
custom->addPerBondParameter("theta0");
custom->addPerBondParameter("r0");
vector<double> parameters(4);
parameters[0] = 1.0;
parameters[1] = 1.0;
parameters[2] = 2 * M_PI / 3;
parameters[3] = sqrt(3) / 2;
vector<int> particles0(3);
particles0[0] = 0;
particles0[1] = 1;
particles0[2] = 2;
vector<int> particles1(3);
particles1[0] = 1;
particles1[1] = 2;
particles1[2] = 3;
custom->addBond(particles0, parameters);
custom->addBond(particles1, parameters);
customSystem.addForce(custom);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0.5, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0.5, 0, 0);
positions[3] = Vec3(0.6, 0, 0.4);
VerletIntegrator integrator(0.01);
Context context(customSystem, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(0.199, state.getPotentialEnergy(), 1e-3);
vector<Vec3> forces(state.getForces());
ASSERT_EQUAL_VEC(Vec3(-1.160, 0.112, 0.0), forces[0], 1e-3);
ASSERT_EQUAL_VEC(Vec3(0.927, 1.047, -0.638), forces[1], 1e-3);
ASSERT_EQUAL_VEC(Vec3(-0.543, -1.160, 0.721), forces[2], 1e-3);
ASSERT_EQUAL_VEC(Vec3(0.776, 0.0, -0.084), forces[3], 1e-3);
}
int main() {
try {
testBond();
testPositionDependence();
testContinuous2DFunction();
testContinuous3DFunction();
testMultipleBonds();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp
View file @ 4233befd
......@@ -398,7 +398,7 @@ void testDiscrete1DFunction() {
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());
ASSERT_EQUAL_TOL(table[i]+1.0, state.getPotentialEnergy(), 1e-6);
}
}
......@@ -431,7 +431,7 @@ void testDiscrete2DFunction() {
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());
ASSERT_EQUAL_TOL(table[i]+1.0, state.getPotentialEnergy(), 1e-6);
}
}
......@@ -468,7 +468,7 @@ void testDiscrete3DFunction() {
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());
ASSERT_EQUAL_TOL(table[i]+1.0, state.getPotentialEnergy(), 1e-6);
}
}
......
platforms/reference/tests/TestReferenceEwald.cpp
View file @ 4233befd
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2010 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -354,6 +354,55 @@ void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
}
}
void testPMEParameters() {
// Create a cloud of random point charges.
const int numParticles = 51;
const double boxWidth = 4.7;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
}
force->setNonbondedMethod(NonbondedForce::PME);
ReferencePlatform platform;
// Compute the energy with an error tolerance of 1e-3.
force->setEwaldErrorTolerance(1e-3);
VerletIntegrator integrator1(0.01);
Context context1(system, integrator1, platform);
context1.setPositions(positions);
double energy1 = context1.getState(State::Energy).getPotentialEnergy();
// Try again with an error tolerance of 1e-4.
force->setEwaldErrorTolerance(1e-4);
VerletIntegrator integrator2(0.01);
Context context2(system, integrator2, platform);
context2.setPositions(positions);
double energy2 = context2.getState(State::Energy).getPotentialEnergy();
// Now explicitly set the parameters. These should match the values that were
// used for tolerance 1e-3.
force->setPMEParameters(2.49291157051793, 32, 32, 32);
VerletIntegrator integrator3(0.01);
Context context3(system, integrator3, platform);
context3.setPositions(positions);
double energy3 = context3.getState(State::Energy).getPotentialEnergy();
ASSERT_EQUAL_TOL(energy1, energy3, 1e-6);
ASSERT(fabs((energy1-energy2)/energy1) > 1e-5);
}
int main() {
try {
testEwaldExact();
......@@ -362,6 +411,7 @@ int main() {
// testWaterSystem();
testErrorTolerance(NonbondedForce::Ewald);
testErrorTolerance(NonbondedForce::PME);
testPMEParameters();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
wrappers/python/simtk/openmm/app/desmonddmsfile.py
View file @ 4233befd
'''
"""
desmonddmsfile.py: Load Desmond dms files
Portions copyright (c) 2013 Stanford University and the Authors
......@@ -23,7 +23,7 @@ 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.
'''
"""
import os
import math
......@@ -38,16 +38,16 @@ from simtk.unit import (nanometer, angstrom, dalton, radian,
class DesmondDMSFile(object):
'''DesmondDMSFile parses a Desmond DMS (desmond molecular system) and
"""DesmondDMSFile parses a Desmond DMS (desmond molecular system) and
constructs a topology and (optionally) an OpenMM System from it
'''
"""
def __init__(self, file):
'''Load a DMS file
"""Load a DMS file
Parameters:
- file (string) the name of the file to load
'''
- file (string) the name of the file to load
"""
# sqlite3 is included in the standard lib, but at python
# compile time, you can disable support (I think), so it's
......@@ -73,8 +73,8 @@ class DesmondDMSFile(object):
# build the provenance string
provenance = []
q = '''SELECT id, user, timestamp, version, workdir, cmdline, executable
FROM provenance'''
q = """SELECT id, user, timestamp, version, workdir, cmdline, executable
FROM provenance"""
#for id, user, timestamp, version, workdir, cmdline, executable in self._conn.execute(q):
for row in self._conn.execute('SELECT * FROM provenance'):
rowdict = dict(zip(self._tables['provenance'], row))
......@@ -89,23 +89,23 @@ class DesmondDMSFile(object):
self._angleConstraints = [{} for x in range(len(self._topologyAtoms))]
def getPositions(self):
'''Get the positions of each atom in the system
'''
"""Get the positions of each atom in the system
"""
return self.positions
def getTopology(self):
'''Get the topology of the system
'''
"""Get the topology of the system
"""
return self.topology
def getProvenance(self):
'''Get the provenance string of this system
'''
"""Get the provenance string of this system
"""
return self.provenance
def _createTopology(self):
'''Build the topology of the system
'''
"""Build the topology of the system
"""
top = Topology()
positions = []
......@@ -119,8 +119,8 @@ class DesmondDMSFile(object):
lastChain = None
lastResId = None
c = top.addChain()
q = '''SELECT id, name, anum, resname, resid, chain, x, y, z
FROM particle'''
q = """SELECT id, name, anum, resname, resid, chain, x, y, z
FROM particle"""
for (atomId, atomName, atomNumber, resName, resId, chain, x, y, z) in self._conn.execute(q):
newChain = False
if chain != lastChain:
......@@ -156,17 +156,17 @@ class DesmondDMSFile(object):
def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*nanometer,
ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None):
'''Construct an OpenMM System representing the topology described by this dms file
"""Construct an OpenMM System representing the topology described by this dms file
Parameters:
- nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
- nonbondedCutoff (distance=1*nanometer) The cutoff distance to use for nonbonded interactions
- ewaldErrorTolerance (float=0.0005) The error tolerance to use if nonbondedMethod is Ewald or PME.
- removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System
- hydrogenMass (mass=None) The mass to use for hydrogen atoms bound to heavy atoms. Any mass added to a hydrogen is
subtracted from the heavy atom to keep their total mass the same.
'''
- nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
- nonbondedCutoff (distance=1*nanometer) The cutoff distance to use for nonbonded interactions
- ewaldErrorTolerance (float=0.0005) The error tolerance to use if nonbondedMethod is Ewald or PME.
- removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System
- hydrogenMass (mass=None) The mass to use for hydrogen atoms bound to heavy atoms. Any mass added to a hydrogen is
subtracted from the heavy atom to keep their total mass the same.
"""
self._checkForUnsupportedTerms()
sys = mm.System()
......@@ -219,14 +219,14 @@ class DesmondDMSFile(object):
return sys
def _addBondsToSystem(self, sys):
'''Create the harmonic bonds
'''
"""Create the harmonic bonds
"""
bonds = mm.HarmonicBondForce()
sys.addForce(bonds)
q = '''SELECT p0, p1, r0, fc, constrained
q = """SELECT p0, p1, r0, fc, constrained
FROM stretch_harm_term INNER JOIN stretch_harm_param
ON stretch_harm_term.param=stretch_harm_param.id'''
ON stretch_harm_term.param=stretch_harm_param.id"""
for p0, p1, r0, fc, constrained in self._conn.execute(q):
if constrained:
sys.addConstraint(p0, p1, r0*angstrom)
......@@ -242,15 +242,15 @@ class DesmondDMSFile(object):
return bonds
def _addAnglesToSystem(self, sys):
'''Create the harmonic angles
'''
"""Create the harmonic angles
"""
angles = mm.HarmonicAngleForce()
sys.addForce(angles)
degToRad = math.pi/180
q = '''SELECT p0, p1, p2, theta0, fc, constrained
q = """SELECT p0, p1, p2, theta0, fc, constrained
FROM angle_harm_term INNER JOIN angle_harm_param
ON angle_harm_term.param=angle_harm_param.id'''
ON angle_harm_term.param=angle_harm_param.id"""
for p0, p1, p2, theta0, fc, constrained in self._conn.execute(q):
if constrained:
l1 = self._atomBonds[p1][p0]
......@@ -267,15 +267,15 @@ class DesmondDMSFile(object):
return angles
def _addConstraintsToSystem(self, sys):
'''Add constraints to system. Normally these should already be
"""Add constraints to system. Normally these should already be
added by the bonds table, but we want to make sure that there's
no extra information in the constraints table that we're not
including in the system'''
including in the system"""
for term_table in [n for n in self._tables.keys() if n.startswith('constraint_a') and n.endswith('term')]:
param_table = term_table.replace('term', 'param')
q = '''SELECT p0, p1, r1
q = """SELECT p0, p1, r1
FROM %(term)s INNER JOIN %(param)s
ON %(term)s.param=%(param)s.id''' % \
ON %(term)s.param=%(param)s.id""" % \
{'term': term_table, 'param': param_table}
for p0, p1, r1 in self._conn.execute(q):
if not p1 in self._atomBonds[p0]:
......@@ -285,9 +285,9 @@ class DesmondDMSFile(object):
if 'constraint_hoh_term' in self._tables:
degToRad = math.pi/180
q = '''SELECT p0, p1, p2, r1, r2, theta
q = """SELECT p0, p1, p2, r1, r2, theta
FROM constraint_hoh_term INNER JOIN constraint_hoh_param
ON constraint_hoh_term.param=constraint_hoh_param.id'''
ON constraint_hoh_term.param=constraint_hoh_param.id"""
for p0, p1, p2, r1, r2, theta in self._conn.execute(q):
# Here, p0 is the heavy atom and p1 and p2 are the H1 and H2
# wihth O-H1 and O-H2 distances r1 and r2
......@@ -296,14 +296,14 @@ class DesmondDMSFile(object):
sys.addConstraint(p1, p2, length)
def _addPeriodicTorsionsToSystem(self, sys):
'''Create the torsion terms
'''
"""Create the torsion terms
"""
periodic = mm.PeriodicTorsionForce()
sys.addForce(periodic)
q = '''SELECT p0, p1, p2, p3, phi0, fc0, fc1, fc2, fc3, fc4, fc5, fc6
q = """SELECT p0, p1, p2, p3, phi0, fc0, fc1, fc2, fc3, fc4, fc5, fc6
FROM dihedral_trig_term INNER JOIN dihedral_trig_param
ON dihedral_trig_term.param=dihedral_trig_param.id'''
ON dihedral_trig_term.param=dihedral_trig_param.id"""
for p0, p1, p2, p3, phi0, fc0, fc1, fc2, fc3, fc4, fc5, fc6 in self._conn.execute(q):
for order, fc in enumerate([fc0, fc1, fc2, fc3, fc4, fc5, fc6]):
if fc == 0:
......@@ -312,8 +312,8 @@ class DesmondDMSFile(object):
def _addImproperHarmonicTorsionsToSystem(self, sys):
'''Create the improper harmonic torsion terms
'''
"""Create the improper harmonic torsion terms
"""
if not self._hasTable('improper_harm_term'):
return
......@@ -322,15 +322,15 @@ class DesmondDMSFile(object):
harmonicTorsion.addPerTorsionParameter('k')
sys.addForce(harmonicTorsion)
q = '''SELECT p0, p1, p2, p3, phi0, fc
q = """SELECT p0, p1, p2, p3, phi0, fc
FROM improper_harm_term INNER JOIN improper_harm_param
ON improper_harm_term.param=improper_harm_param.id'''
ON improper_harm_term.param=improper_harm_param.id"""
for p0, p1, p2, p3, phi0, fc in self._conn.execute(q):
harmonicTorsion.addTorsion(p0, p1, p2, p3, [phi0*degree, fc*kilocalorie_per_mole])
def _addCMAPToSystem(self, sys):
'''Create the CMAP terms
'''
"""Create the CMAP terms
"""
if not self._hasTable('torsiontorsion_cmap_term'):
return
......@@ -354,30 +354,30 @@ class DesmondDMSFile(object):
index = cmap.addMap(size, map*kilocalorie_per_mole)
cmap_indices[name] = index
q = '''SELECT p0, p1, p2, p3, p4, p5, p6, p7, cmapid
q = """SELECT p0, p1, p2, p3, p4, p5, p6, p7, cmapid
FROM torsiontorsion_cmap_term INNER JOIN torsiontorsion_cmap_param
ON torsiontorsion_cmap_term.param=torsiontorsion_cmap_param.id'''
ON torsiontorsion_cmap_term.param=torsiontorsion_cmap_param.id"""
for p0, p1, p2, p3, p4, p5, p6, p7, cmapid in self._conn.execute(q):
cmap.addTorsion(cmap_indices[cmapid], p0, p1, p2, p3, p4, p5, p6, p7)
def _addNonbondedForceToSystem(self, sys):
'''Create the nonbonded force
'''
"""Create the nonbonded force
"""
nb = mm.NonbondedForce()
sys.addForce(nb)
q = '''SELECT charge, sigma, epsilon
q = """SELECT charge, sigma, epsilon
FROM particle INNER JOIN nonbonded_param
ON particle.nbtype=nonbonded_param.id'''
ON particle.nbtype=nonbonded_param.id"""
for charge, sigma, epsilon in self._conn.execute(q):
nb.addParticle(charge, sigma*angstrom, epsilon*kilocalorie_per_mole)
for p0, p1 in self._conn.execute('SELECT p0, p1 FROM exclusion'):
nb.addException(p0, p1, 0.0, 1.0, 0.0)
q = '''SELECT p0, p1, aij, bij, qij
q = """SELECT p0, p1, aij, bij, qij
FROM pair_12_6_es_term INNER JOIN pair_12_6_es_param
ON pair_12_6_es_term.param=pair_12_6_es_param.id;'''
ON pair_12_6_es_term.param=pair_12_6_es_param.id;"""
for p0, p1, a_ij, b_ij, q_ij in self._conn.execute(q):
a_ij = (a_ij*kilocalorie_per_mole*(angstrom**12)).in_units_of(kilojoule_per_mole*(nanometer**12))
b_ij = (b_ij*kilocalorie_per_mole*(angstrom**6)).in_units_of(kilojoule_per_mole*(nanometer**6))
......@@ -391,41 +391,41 @@ class DesmondDMSFile(object):
new_sigma = (a_ij / b_ij)**(1.0/6.0)
nb.addException(p0, p1, q_ij, new_sigma, new_epsilon, True)
n_total = self._conn.execute('''SELECT COUNT(*) FROM pair_12_6_es_term''').fetchone()
n_in_exclusions = self._conn.execute('''SELECT COUNT(*)
n_total = self._conn.execute("""SELECT COUNT(*) FROM pair_12_6_es_term""").fetchone()
n_in_exclusions = self._conn.execute("""SELECT COUNT(*)
FROM exclusion INNER JOIN pair_12_6_es_term
ON exclusion.p0==pair_12_6_es_term.p0 AND exclusion.p1==pair_12_6_es_term.p1''').fetchone()
ON exclusion.p0==pair_12_6_es_term.p0 AND exclusion.p1==pair_12_6_es_term.p1""").fetchone()
if not n_total == n_in_exclusions:
raise NotImplementedError('All pair_12_6_es_terms must have a corresponding exclusion')
return nb
def _addVirtualSitesToSystem(self, sys):
'''Create any virtual sites in the systempy
'''
"""Create any virtual sites in the systempy
"""
if not any(t.startswith('virtual_') for t in self._tables.keys()):
return
if 'virtual_lc2_term' in self._tables:
q = '''SELECT p0, p1, p2, c1
q = """SELECT p0, p1, p2, c1
FROM virtual_lc2_term INNER JOIN virtual_lc2_param
ON virtual_lc2_term.param=virtual_lc2_param.id'''
ON virtual_lc2_term.param=virtual_lc2_param.id"""
for p0, p1, p2, c1 in self._conn.execute(q):
vsite = mm.TwoParticleAverageSite(p1, p2, (1-c1), c1)
sys.setVirtualSite(p0, vsite)
if 'virtual_lc3_term' in self._tables:
q = '''SELECT p0, p1, p2, p3, c1, c2
q = """SELECT p0, p1, p2, p3, c1, c2
FROM virtual_lc3_term INNER JOIN virtual_lc3_param
ON virtual_lc3_term.param=virtual_lc3_param.id'''
ON virtual_lc3_term.param=virtual_lc3_param.id"""
for p0, p1, p2, p3, c1, c2 in self._conn.execute(q):
vsite = mm.ThreeParticleAverageSite(p1, p2, p3, (1-c1-c2), c1, c2)
sys.setVirtualSite(p0, vsite)
if 'virtual_out3_term' in self._tables:
q = '''SELECT p0, p1, p2, p3, c1, c2, c3
q = """SELECT p0, p1, p2, p3, c1, c2, c3
FROM virtual_out3_term INNER JOIN virtual_out3_param
ON virtual_out3_term.param=virtual_out3_param.id'''
ON virtual_out3_term.param=virtual_out3_param.id"""
for p0, p1, p2, p3, c1, c2, c3 in self._conn.execute(q):
vsite = mm.OutOfPlaneSite(p1, p2, p3, c1, c2, c3)
sys.setVirtualSite(p0, vsite)
......@@ -436,14 +436,14 @@ class DesmondDMSFile(object):
def _hasTable(self, table_name):
'''Does our DMS file contain this table?
'''
"""Does our DMS file contain this table?
"""
return table_name in self._tables
def _readSchemas(self):
'''Read the schemas of each of the tables in the dms file, populating
"""Read the schemas of each of the tables in the dms file, populating
the `_tables` instance attribute
'''
"""
tables = {}
for table in self._conn.execute("SELECT name FROM sqlite_master WHERE type='table'"):
names = []
......@@ -453,9 +453,9 @@ class DesmondDMSFile(object):
self._tables = tables
def _checkForUnsupportedTerms(self):
'''Check the file for forcefield terms that are not currenty supported,
"""Check the file for forcefield terms that are not currenty supported,
raising a NotImplementedError
'''
"""
if 'posre_harm_term' in self._tables:
raise NotImplementedError('Position restraints are not implemented.')
flat_bottom_potential_terms = ['stretch_fbhw_term', 'angle_fbhw_term',
......@@ -488,8 +488,8 @@ class DesmondDMSFile(object):
raise NotImplementedError('Drude particles are not currently supported')
def close(self):
'''Close the SQL connection
'''
"""Close the SQL connection
"""
if self._open:
self._conn.close()
......
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