Commit 86fc3001 authored by Christopher Bruns's avatar Christopher Bruns
Browse files

Started making some improvements to installer behavior.

Including installing Test<Platform>HarmonicBondForce tests.
parent 14125aa9
...@@ -425,7 +425,10 @@ ELSE (EXECUTABLE_OUTPUT_PATH) ...@@ -425,7 +425,10 @@ ELSE (EXECUTABLE_OUTPUT_PATH)
SET (TEST_PATH .) SET (TEST_PATH .)
ENDIF (EXECUTABLE_OUTPUT_PATH) ENDIF (EXECUTABLE_OUTPUT_PATH)
INCLUDE(ApiDoxygen.cmake) # Build and install documentation
INCLUDE(cmake_modules/ApiDoxygen.cmake)
install(FILES docs/UserGuide/OpenMMUserGuide.pdf
DESTINATION docs)
ADD_SUBDIRECTORY(tests) ADD_SUBDIRECTORY(tests)
ADD_SUBDIRECTORY(examples) ADD_SUBDIRECTORY(examples)
......
...@@ -57,4 +57,6 @@ SET(INCLUDE_CUDA_STATIC 1) ...@@ -57,4 +57,6 @@ SET(INCLUDE_CUDA_STATIC 1)
# ADD_TEST(${TEST_STATIC} "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFileStatic" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/bptiMdRfPbcParameters.txt" " +checkEnergyForceConsistent -checkForces" ) # ADD_TEST(${TEST_STATIC} "${EXECUTABLE_OUTPUT_PATH}/TestCudaUsingParameterFileStatic" "-parameterFileName" "${CMAKE_CURRENT_SOURCE_DIR}/bptiMdRfPbcParameters.txt" " +checkEnergyForceConsistent -checkForces" )
#ENDIF(INCLUDE_CUDA_STATIC) #ENDIF(INCLUDE_CUDA_STATIC)
install(TARGETS TestCudaHarmonicBondForce
RUNTIME DESTINATION bin)
...@@ -75,12 +75,14 @@ void testBonds() { ...@@ -75,12 +75,14 @@ void testBonds() {
int main() { int main() {
try { try {
cout << "Running test..." << endl;
testBonds(); testBonds();
} }
catch(const exception& e) { catch(const exception& e) {
cout << "exception: " << e.what() << endl; cout << "exception: " << e.what() << endl;
cout << "FAIL - ERROR. Test failed." << endl;
return 1; return 1;
} }
cout << "Done" << endl; cout << "PASS - Test succeeded." << endl;
return 0; return 0;
} }
...@@ -9,22 +9,13 @@ FILE(GLOB TEST_PROGS "*Test*.cpp") ...@@ -9,22 +9,13 @@ FILE(GLOB TEST_PROGS "*Test*.cpp")
FOREACH(TEST_PROG ${TEST_PROGS}) FOREACH(TEST_PROG ${TEST_PROGS})
GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE) GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE)
# All tests use shared libraries
# Link with shared library
ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG}) ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET}) TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET})
ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT}) ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
# Link with static library
# SET(TEST_STATIC ${TEST_ROOT}Static)
# ADD_EXECUTABLE(${TEST_STATIC} ${TEST_PROG})
# SET_TARGET_PROPERTIES(${TEST_STATIC}
# PROPERTIES
# COMPILE_FLAGS "-DOPENMM_USE_STATIC_LIBRARIES"
# )
# TARGET_LINK_LIBRARIES(${TEST_STATIC} ${STATIC_TARGET})
# ADD_TEST(${TEST_STATIC} ${EXECUTABLE_OUTPUT_PATH}/${TEST_STATIC})
ENDFOREACH(TEST_PROG ${TEST_PROGS}) ENDFOREACH(TEST_PROG ${TEST_PROGS})
install(TARGETS TestReferenceHarmonicBondForce
RUNTIME DESTINATION bin)
...@@ -75,12 +75,14 @@ void testBonds() { ...@@ -75,12 +75,14 @@ void testBonds() {
int main() { int main() {
try { try {
cout << "Running test..." << endl;
testBonds(); testBonds();
} }
catch(const exception& e) { catch(const exception& e) {
cout << "exception: " << e.what() << endl; cout << "exception: " << e.what() << endl;
cout << "FAIL - ERROR. Test failed." << endl;
return 1; return 1;
} }
cout << "Done" << endl; cout << "PASS - Test succeeded." << endl;
return 0; return 0;
} }
...@@ -8,25 +8,9 @@ ENABLE_TESTING() ...@@ -8,25 +8,9 @@ ENABLE_TESTING()
FILE(GLOB TEST_PROGS "*Test*.cpp") FILE(GLOB TEST_PROGS "*Test*.cpp")
FOREACH(TEST_PROG ${TEST_PROGS}) FOREACH(TEST_PROG ${TEST_PROGS})
GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE) GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE)
# All tests use shared library only
# Link with shared library
ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG}) ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET}) TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET})
ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT}) ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
# Link with static library
# SET(TEST_STATIC ${TEST_ROOT}Static)
# ADD_EXECUTABLE(${TEST_STATIC} ${TEST_PROG})
# SET_TARGET_PROPERTIES(${TEST_STATIC}
# PROPERTIES
# COMPILE_FLAGS "-DOPENMM_USE_STATIC_LIBRARIES"
# )
# TARGET_LINK_LIBRARIES(${TEST_STATIC} ${STATIC_TARGET})
# ADD_TEST(${TEST_STATIC} ${EXECUTABLE_OUTPUT_PATH}/${TEST_STATIC})
ENDFOREACH(TEST_PROG ${TEST_PROGS}) ENDFOREACH(TEST_PROG ${TEST_PROGS})
#ADD_SUBDIRECTORY(platforms)
#
# Test each supported platform with TestPlatform*.cpp programs
#
ENABLE_TESTING()
# For a particular platform, for example "ReferencePlatform", there must be:
# 1) A C++ object with that name derived from OpenMM::Platform
# 2) A header file with name like "ReferencePlatform.h"
SET(SUPPORTED_PLATFORMS "ReferencePlatform" CACHE STRING "Hardware platforms supported on this computer")
# Automatically create tests using files named "TestPlatform*.cpp.in"
FILE(GLOB PLATFORM_GENERIC_PROGS "TestPlatform*.cpp.in")
FOREACH(PLATFORM_GENERIC_PROG ${PLATFORM_GENERIC_PROGS})
# Create platform-specific test for each supported platform
FOREACH(OPENMM_PLATFORM_NAME ${SUPPORTED_PLATFORMS})
STRING(REGEX REPLACE "TestPlatform([^\\/]*)\\.cpp.in" "\\1" TEST_NAME ${PLATFORM_GENERIC_PROG})
GET_FILENAME_COMPONENT(TEST_NAME ${TEST_NAME} NAME_WE)
SET(PLATFORM_SPECIFIC_PROG "${CMAKE_CURRENT_BINARY_DIR}/Test${OPENMM_PLATFORM_NAME}${TEST_NAME}.cpp")
CONFIGURE_FILE(${PLATFORM_GENERIC_PROG} ${PLATFORM_SPECIFIC_PROG} @ONLY)
GET_FILENAME_COMPONENT(TEST_EXECUTABLE ${PLATFORM_SPECIFIC_PROG} NAME_WE)
# Link with shared library
ADD_EXECUTABLE(${TEST_EXECUTABLE} ${PLATFORM_SPECIFIC_PROG})
TARGET_LINK_LIBRARIES(${TEST_EXECUTABLE} ${SHARED_TARGET})
ADD_TEST(${TEST_EXECUTABLE} ${EXECUTABLE_OUTPUT_PATH}/${TEST_EXECUTABLE})
ENDFOREACH(OPENMM_PLATFORM_NAME ${SUPPORTED_PLATFORMS})
ENDFOREACH(PLATFORM_GENERIC_PROG ${PLATFORM_GENERIC_PROGS})
/* -------------------------------------------------------------------------- *
* 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 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 reference implementation of AndersenThermostat.
*/
#include "@OPENMM_PLATFORM_NAME@.h"
#include "../../../tests/AssertionUtilities.h"
#include "AndersenThermostat.h"
#include "OpenMMContext.h"
#include "StandardMMForceField.h"
#include "System.h"
#include "VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "../src/sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
void testTemperature(Platform& platform)
{
const int numAtoms = 8;
const double temp = 100.0;
const double collisionFreq = 10.0;
System system(numAtoms, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 2.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
AndersenThermostat* thermstat = new AndersenThermostat(temp, collisionFreq);
system.addForce(thermstat);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
for (int i = 0; i < numAtoms; ++i)
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
}
ke /= 1000;
double expected = 0.5*numAtoms*3*BOLTZ*temp;
ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
}
int main()
{
try {
@OPENMM_PLATFORM_NAME@ platform;
testTemperature(platform);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
cout << "*** test FAILED ***" << endl;
return 1;
}
cout << "Done" << endl;
cout << "test passed." << endl;
return 0;
}
/* -------------------------------------------------------------------------- *
* 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 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 reference implementation of BrownianIntegrator.
*/
#include "@OPENMM_PLATFORM_NAME@.h"
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "StandardMMForceField.h"
#include "System.h"
#include "BrownianIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "../src/sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testSingleBond(Platform& platform) {
System system(2, 0);
system.setAtomMass(0, 2.0);
system.setAtomMass(1, 2.0);
double dt = 0.01;
BrownianIntegrator integrator(0, 0.1, dt);
StandardMMForceField* forceField = new StandardMMForceField(2, 1, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1.5, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
// This is simply an overdamped harmonic oscillator, so compare it to the analytical solution.
double rate = 2*1.0/0.1;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Velocities);
double time = state.getTime();
double expectedDist = 1.5+0.5*std::exp(-rate*time);
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedDist, 0, 0), state.getPositions()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedDist, 0, 0), state.getPositions()[1], 0.02);
if (i > 0) {
double expectedSpeed = -0.5*rate*std::exp(-rate*(time-0.5*dt));
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedSpeed, 0, 0), state.getVelocities()[0], 0.11);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedSpeed, 0, 0), state.getVelocities()[1], 0.11);
}
integrator.step(1);
}
}
void testTemperature(Platform& platform) {
const int numAtoms = 8;
const int numBonds = numAtoms-1;
const double temp = 100.0;
System system(numAtoms, 0);
BrownianIntegrator integrator(temp, 2.0, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, numBonds, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 2.0);
// forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
for (int i = 0; i < numBonds; ++i)
forceField->setBondParameters(i, i, i+1, 1.0, i);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
for (int i = 0; i < numAtoms; ++i)
positions[i] = Vec3(i, 0, 0);
context.setPositions(positions);
// Let it equilibrate.
integrator.step(1000);
// Now run it for a while and see if the temperature is correct.
double pe = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Energy);
pe += state.getPotentialEnergy();
integrator.step(1);
}
pe /= 1000;
double expected = 0.5*numBonds*BOLTZ*temp;
ASSERT_EQUAL_TOL(expected, pe, 3*expected/std::sqrt(1000.0));
}
void testConstraints(Platform& platform) {
const int numAtoms = 8;
const double temp = 100.0;
System system(numAtoms, numAtoms-1);
BrownianIntegrator integrator(temp, 2.0, 0.001);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 10.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numAtoms-1; ++i)
system.setConstraintParameters(i, i, i+1, 1.0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
vector<Vec3> velocities(numAtoms);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions);
for (int j = 0; j < numAtoms-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
ASSERT_EQUAL_TOL(1.0, dist, 2e-4);
}
integrator.step(1);
}
}
int main() {
try {
@OPENMM_PLATFORM_NAME@ platform;
testSingleBond(platform);
testTemperature(platform);
testConstraints(platform);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
/* -------------------------------------------------------------------------- *
* 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 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 reference implementation of AndersenThermostat.
*/
#include "@OPENMM_PLATFORM_NAME@.h"
#include "../../../tests/AssertionUtilities.h"
#include "CMMotionRemover.h"
#include "OpenMMContext.h"
#include "StandardMMForceField.h"
#include "System.h"
#include "VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "../src/sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
Vec3 calcCM(const vector<Vec3>& values, System& system) {
Vec3 cm;
for (int j = 0; j < system.getNumAtoms(); ++j) {
cm[0] += values[j][0]*system.getAtomMass(j);
cm[1] += values[j][1]*system.getAtomMass(j);
cm[2] += values[j][2]*system.getAtomMass(j);
}
return cm;
}
void testMotionRemoval(Platform& platform) {
const int numAtoms = 8;
const double temp = 100.0;
const double collisionFreq = 10.0;
System system(numAtoms, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 1, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, i+1);
forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
forceField->setBondParameters(0, 2, 3, 2.0, 0.5);
system.addForce(forceField);
CMMotionRemover* remover = new CMMotionRemover();
system.addForce(remover);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
vector<Vec3> velocities(numAtoms);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i) {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Now run it for a while and see if the center of mass remains fixed.
Vec3 cmPos = calcCM(context.getState(State::Positions).getPositions(), system);
for (int i = 0; i < 1000; ++i) {
integrator.step(1);
State state = context.getState(State::Positions | State::Velocities);
Vec3 pos = calcCM(state.getPositions(), system);
ASSERT_EQUAL_VEC(cmPos, pos, 1e-2);
Vec3 vel = calcCM(state.getVelocities(), system);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), vel, 1e-2);
}
}
int main() {
try {
@OPENMM_PLATFORM_NAME@ platform;
testMotionRemoval(platform);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
/* -------------------------------------------------------------------------- *
* 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 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 platform implementation of GBSAOBCForceField.
*/
#include "@OPENMM_PLATFORM_NAME@.h"
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "GBSAOBCForceField.h"
#include "System.h"
#include "LangevinIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "../src/sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testSingleAtom(Platform& platform) {
System system(1, 0);
system.setAtomMass(0, 2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForceField* forceField = new GBSAOBCForceField(1);
forceField->setAtomParameters(0, 0.5, 0.15, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = 0.15-0.009; // dielectric offset
double eps0 = EPSILON0;
double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric())/bornRadius;
double extendedRadius = bornRadius+0.14; // probe radius
double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from? Just copied it from CpuImplicitSolvent.cpp
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
void testForce(Platform& platform) {
const int numAtoms = 10;
System system(numAtoms, 0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForceField* forceField = new GBSAOBCForceField(numAtoms);
for (int i = 0; i < numAtoms; ++i)
forceField->setAtomParameters(i, i%2 == 0 ? -1 : 1, 0.15, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
// Set random positions for all the atoms.
vector<Vec3> positions(numAtoms);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i)
positions[i] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
// Take a small step in the direction of the energy gradient.
double norm = 0.0;
for (int i = 0; i < numAtoms; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-3;
double step = delta/norm;
for (int i = 0; i < numAtoms; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
// See whether the potential energy changed by the expected amount.
State state2 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta, 0.01)
}
int main() {
try {
@OPENMM_PLATFORM_NAME@ platform;
testSingleAtom(platform);
testForce(platform);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
/* -------------------------------------------------------------------------- *
* 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 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 platform implementation of the kernel to calculate kinetic energy.
*/
#include "@OPENMM_PLATFORM_NAME@.h"
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "System.h"
#include "VerletIntegrator.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testCalcKE(Platform& platform) {
System system(4, 0);
for (int i = 0; i < 4; ++i)
system.setAtomMass(i, i+1);
VerletIntegrator integrator(0.01);
OpenMMContext context(system, integrator, platform);
vector<Vec3> velocities(4);
velocities[0] = Vec3(1, 0, 0);
velocities[1] = Vec3(0, 1, 0);
velocities[2] = Vec3(0, 0, 2);
velocities[3] = Vec3(std::sqrt(2.0), 0, std::sqrt(2.0));
context.setVelocities(velocities);
State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(0.5*(1+2+4*3+4*4), state.getKineticEnergy(), TOL);
}
int main() {
try {
@OPENMM_PLATFORM_NAME@ platform;
testCalcKE(platform);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
/* -------------------------------------------------------------------------- *
* 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 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 platform implementation of LangevinIntegrator.
*/
#include "@OPENMM_PLATFORM_NAME@.h"
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "StandardMMForceField.h"
#include "System.h"
#include "LangevinIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "../src/sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testSingleBond(Platform& platform) {
System system(2, 0);
system.setAtomMass(0, 2.0);
system.setAtomMass(1, 2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 1, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1.5, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
// This is simply a damped harmonic oscillator, so compare it to the analytical solution.
double freq = std::sqrt(1-0.05*0.05);
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Velocities);
double time = state.getTime();
double expectedDist = 1.5+0.5*std::exp(-0.05*time)*std::cos(freq*time);
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedDist, 0, 0), state.getPositions()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedDist, 0, 0), state.getPositions()[1], 0.02);
double expectedSpeed = -0.5*std::exp(-0.05*time)*(0.05*std::cos(freq*time)+freq*std::sin(freq*time));
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedSpeed, 0, 0), state.getVelocities()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedSpeed, 0, 0), state.getVelocities()[1], 0.02);
integrator.step(1);
}
// Not set the friction to a tiny value and see if it conserves energy.
integrator.setFriction(5e-5);
context.setPositions(positions);
State state = context.getState(State::Energy);
double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 1000; ++i) {
state = context.getState(State::Energy);
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testTemperature(Platform& platform) {
const int numAtoms = 8;
const double temp = 100.0;
System system(numAtoms, 0);
LangevinIntegrator integrator(temp, 2.0, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 2.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
for (int i = 0; i < numAtoms; ++i)
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
}
ke /= 1000;
double expected = 0.5*numAtoms*3*BOLTZ*temp;
ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
}
void testConstraints(Platform& platform) {
const int numAtoms = 8;
const double temp = 100.0;
System system(numAtoms, numAtoms-1);
LangevinIntegrator integrator(temp, 2.0, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 10.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numAtoms-1; ++i)
system.setConstraintParameters(i, i, i+1, 1.0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
vector<Vec3> velocities(numAtoms);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions);
for (int j = 0; j < numAtoms-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
ASSERT_EQUAL_TOL(1.0, dist, 2e-4);
}
integrator.step(1);
}
}
int main() {
try {
@OPENMM_PLATFORM_NAME@ platform;
testSingleBond(platform);
testTemperature(platform);
testConstraints(platform);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
/* -------------------------------------------------------------------------- *
* 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 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 all the different force terms in the platform implementation of StandardMMForceField.
*/
#include "@OPENMM_PLATFORM_NAME@.h"
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "StandardMMForceField.h"
#include "System.h"
#include "VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testBonds(Platform& platform) {
System system(3, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 2, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1.5, 0.8);
forceField->setBondParameters(1, 1, 2, 1.2, 0.7);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 2, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.5, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.7*0.2, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
ASSERT_EQUAL_TOL(0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
}
void testAngles(Platform& platform) {
System system(4, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 2, 0, 0);
forceField->setAngleParameters(0, 0, 1, 2, PI_M/3, 1.1);
forceField->setAngleParameters(1, 1, 2, 3, PI_M/2, 1.2);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(2, 1, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double torque1 = 1.1*PI_M/6;
double torque2 = 1.2*PI_M/4;
ASSERT_EQUAL_VEC(Vec3(torque1, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-0.5*torque2, 0.5*torque2, 0), forces[3], TOL); // reduced by sqrt(2) due to the bond length, another sqrt(2) due to the angle
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(0.5*1.1*(PI_M/6)*(PI_M/6) + 0.5*1.2*(PI_M/4)*(PI_M/4), state.getPotentialEnergy(), TOL);
}
void testPeriodicTorsions(Platform& platform) {
System system(4, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 1, 0);
forceField->setPeriodicTorsionParameters(0, 0, 1, 2, 3, 2, PI_M/3, 1.1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(1, 0, 2);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double torque = -2*1.1*std::sin(2*PI_M/3);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), TOL);
}
void testRBTorsions(Platform& platform) {
System system(4, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 0, 1);
forceField->setRBTorsionParameters(0, 0, 1, 2, 3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(1, 1, 1);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double psi = 0.25*PI_M - PI_M;
double torque = 0.0;
for (int i = 1; i < 6; ++i) {
double c = 0.1*(i+1);
torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
}
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, -0.5*torque), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
double energy = 0.0;
for (int i = 0; i < 6; ++i) {
double c = 0.1*(i+1);
energy += c*std::pow(std::cos(psi), i);
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
void testCoulomb(Platform& platform) {
System system(2, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0);
forceField->setAtomParameters(0, 0.5, 1, 0);
forceField->setAtomParameters(1, -1.5, 1, 0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double force = 138.935485*(-0.75)/4.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(138.935485*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
}
void testLJ(Platform& platform) {
System system(2, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0);
forceField->setAtomParameters(0, 0, 1.2, 1);
forceField->setAtomParameters(1, 0, 1.4, 2);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double x = 1.3/2.0;
double eps = SQRT_TWO;
double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/2.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0)), state.getPotentialEnergy(), TOL);
}
void testExclusionsAnd14(Platform& platform) {
System system(5, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1, 0);
forceField->setBondParameters(1, 1, 2, 1, 0);
forceField->setBondParameters(2, 2, 3, 1, 0);
forceField->setBondParameters(3, 3, 4, 1, 0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(5);
const double r = 1.0;
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(r, 0, 0);
positions[2] = Vec3(r, 0, 0);
positions[3] = Vec3(r, 0, 0);
positions[4] = Vec3(r, 0, 0);
for (int i = 1; i < 5; ++i) {
// Test LJ forces
forceField->setAtomParameters(0, 0, 1.5, 1);
for (int j = 1; j < 5; ++j)
forceField->setAtomParameters(j, 0, 1.5, 0);
forceField->setAtomParameters(i, 0, 1.5, 1);
context.reinitialize();
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double x = 1.5/r;
double eps = 1.0;
double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
double energy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
if (i == 3) {
force *= 0.5;
energy *= 0.5;
}
if (i < 3) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
// Test Coulomb forces
forceField->setAtomParameters(0, 2, 1.5, 0);
forceField->setAtomParameters(i, 2, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces();
force = 138.935485*4/(r*r);
energy = 138.935485*4/r;
if (i == 3) {
force /= 1.2;
energy /= 1.2;
}
if (i < 3) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
}
int main() {
try {
@OPENMM_PLATFORM_NAME@ platform;
testBonds(platform);
testAngles(platform);
testPeriodicTorsions(platform);
testRBTorsions(platform);
testCoulomb(platform);
testLJ(platform);
testExclusionsAnd14(platform);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
/* -------------------------------------------------------------------------- *
* 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 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 platform implementation of VerletIntegrator.
*/
#include "@OPENMM_PLATFORM_NAME@.h"
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "StandardMMForceField.h"
#include "System.h"
#include "VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "../src/sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testSingleBond(Platform& platform) {
System system(2, 0);
system.setAtomMass(0, 2.0);
system.setAtomMass(1, 2.0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 1, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1.5, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
// This is simply a harmonic oscillator, so compare it to the analytical solution.
const double freq = 1.0;;
State state = context.getState(State::Energy);
const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 1000; ++i) {
state = context.getState(State::Positions | State::Velocities | State::Energy);
double time = state.getTime();
double expectedDist = 1.5+0.5*std::cos(freq*time);
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedDist, 0, 0), state.getPositions()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedDist, 0, 0), state.getPositions()[1], 0.02);
double expectedSpeed = -0.5*freq*std::sin(freq*time);
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedSpeed, 0, 0), state.getVelocities()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedSpeed, 0, 0), state.getVelocities()[1], 0.02);
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testConstraints(Platform& platform) {
const int numAtoms = 8;
const double temp = 100.0;
System system(numAtoms, numAtoms-1);
VerletIntegrator integrator(0.002);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 10.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numAtoms-1; ++i)
system.setConstraintParameters(i, i, i+1, 1.0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
vector<Vec3> velocities(numAtoms);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy);
for (int j = 0; j < numAtoms-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
ASSERT_EQUAL_TOL(1.0, dist, 2e-4);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
int main() {
try {
@OPENMM_PLATFORM_NAME@ platform;
testSingleBond(platform);
testConstraints(platform);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment