Skip to content

GitLab

Commit ad75a390 authored by Peter Eastman's avatar Peter Eastman
Browse files

Split StandardMMForceField into separate classes for each force term.

parent 8078c417
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1026 additions and 439 deletions
platforms/cuda/tests/TestCudaStandardMMForceField.cpp → platforms/cuda/tests/TestCudaNonbondedForce.cpp
View file @ ad75a390
......@@ -30,13 +30,14 @@
* -------------------------------------------------------------------------- */
/**
* This tests all the different force terms in the reference implementation of StandardMMForceField.
* This tests all the different force terms in the reference implementation of NonbondedForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "CudaPlatform.h"
#include "StandardMMForceField.h"
#include "HarmonicBondForce.h"
#include "NonbondedForce.h"
#include "System.h"
#include "LangevinIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
......@@ -48,114 +49,11 @@ using namespace std;
const double TOL = 1e-5;
void testBonds() {
CudaPlatform platform;
System system(3, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 2, 0, 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() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 2, 0, 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() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 1, 0, 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() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 0, 1, 0);
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() {
CudaPlatform platform;
System system(2, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0, 0);
NonbondedForce* forceField = new NonbondedForce(2, 0);
forceField->setAtomParameters(0, 0.5, 1, 0);
forceField->setAtomParameters(1, -1.5, 1, 0);
system.addForce(forceField);
......@@ -176,7 +74,7 @@ void testLJ() {
CudaPlatform platform;
System system(2, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0, 0);
NonbondedForce* forceField = new NonbondedForce(2, 0);
forceField->setAtomParameters(0, 0, 1.2, 1);
forceField->setAtomParameters(1, 0, 1.4, 2);
system.addForce(forceField);
......@@ -199,12 +97,14 @@ void testExclusionsAnd14() {
CudaPlatform platform;
System system(5, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0, 2);
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);
HarmonicBondForce* bonds = new HarmonicBondForce(4);
bonds->setBondParameters(0, 0, 1, 1, 0);
bonds->setBondParameters(1, 1, 2, 1, 0);
bonds->setBondParameters(2, 2, 3, 1, 0);
bonds->setBondParameters(3, 3, 4, 1, 0);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(5, 2);
system.addForce(nonbonded);
for (int i = 1; i < 5; ++i) {
// Test LJ forces
......@@ -212,13 +112,13 @@ void testExclusionsAnd14() {
vector<Vec3> positions(5);
const double r = 1.0;
for (int j = 0; j < 5; ++j) {
forceField->setAtomParameters(j, 0, 1.5, 0);
nonbonded->setAtomParameters(j, 0, 1.5, 0);
positions[j] = Vec3(0, j, 0);
}
forceField->setAtomParameters(0, 0, 1.5, 1);
forceField->setAtomParameters(i, 0, 1.5, 1);
forceField->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
nonbonded->setAtomParameters(0, 0, 1.5, 1);
nonbonded->setAtomParameters(i, 0, 1.5, 1);
nonbonded->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
positions[i] = Vec3(r, 0, 0);
OpenMMContext context(system, integrator, platform);
context.setPositions(positions);
......@@ -242,10 +142,10 @@ void testExclusionsAnd14() {
// Test Coulomb forces
forceField->setAtomParameters(0, 2, 1.5, 0);
forceField->setAtomParameters(i, 2, 1.5, 0);
forceField->setNonbonded14Parameters(0, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
nonbonded->setAtomParameters(0, 2, 1.5, 0);
nonbonded->setAtomParameters(i, 2, 1.5, 0);
nonbonded->setNonbonded14Parameters(0, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
OpenMMContext context2(system, integrator, platform);
context2.setPositions(positions);
state = context2.getState(State::Forces | State::Energy);
......@@ -270,11 +170,11 @@ void testCutoff() {
CudaPlatform platform;
System system(3, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 0, 0, 0, 0, 0);
NonbondedForce* forceField = new NonbondedForce(3, 0);
forceField->setAtomParameters(0, 1.0, 1, 0);
forceField->setAtomParameters(1, 1.0, 1, 0);
forceField->setAtomParameters(2, 1.0, 1, 0);
forceField->setNonbondedMethod(StandardMMForceField::CutoffNonPeriodic);
forceField->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
const double cutoff = 2.9;
forceField->setCutoffDistance(cutoff);
system.addForce(forceField);
......@@ -303,15 +203,17 @@ void testCutoff14() {
CudaPlatform platform;
System system(5, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0, 2);
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);
forceField->setNonbondedMethod(StandardMMForceField::CutoffNonPeriodic);
HarmonicBondForce* bonds = new HarmonicBondForce(4);
bonds->setBondParameters(0, 0, 1, 1, 0);
bonds->setBondParameters(1, 1, 2, 1, 0);
bonds->setBondParameters(2, 2, 3, 1, 0);
bonds->setBondParameters(3, 3, 4, 1, 0);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(5, 2);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
const double cutoff = 3.5;
forceField->setCutoffDistance(cutoff);
system.addForce(forceField);
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(5);
positions[0] = Vec3(0, 0, 0);
......@@ -323,12 +225,12 @@ void testCutoff14() {
// Test LJ forces
forceField->setAtomParameters(0, 0, 1.5, 1);
nonbonded->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);
forceField->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
nonbonded->setAtomParameters(j, 0, 1.5, 0);
nonbonded->setAtomParameters(i, 0, 1.5, 1);
nonbonded->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
context.reinitialize();
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
......@@ -353,10 +255,10 @@ void testCutoff14() {
// Test Coulomb forces
const double q = 0.7;
forceField->setAtomParameters(0, q, 1.5, 0);
forceField->setAtomParameters(i, q, 1.5, 0);
forceField->setNonbonded14Parameters(0, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
nonbonded->setAtomParameters(0, q, 1.5, 0);
nonbonded->setAtomParameters(i, q, 1.5, 0);
nonbonded->setNonbonded14Parameters(0, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
......@@ -384,16 +286,18 @@ void testPeriodic() {
CudaPlatform platform;
System system(3, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 1, 0, 0, 0, 0);
forceField->setAtomParameters(0, 1.0, 1, 0);
forceField->setAtomParameters(1, 1.0, 1, 0);
forceField->setAtomParameters(2, 1.0, 1, 0);
forceField->setBondParameters(0, 0, 1, 1.0, 0.0);
forceField->setNonbondedMethod(StandardMMForceField::CutoffPeriodic);
HarmonicBondForce* bonds = new HarmonicBondForce(1);
bonds->setBondParameters(0, 0, 1, 1, 0);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(3, 0);
nonbonded->setAtomParameters(0, 1.0, 1, 0);
nonbonded->setAtomParameters(1, 1.0, 1, 0);
nonbonded->setAtomParameters(2, 1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 2.0;
forceField->setCutoffDistance(cutoff);
forceField->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
system.addForce(forceField);
nonbonded->setCutoffDistance(cutoff);
nonbonded->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
system.addForce(nonbonded);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
......@@ -414,10 +318,6 @@ void testPeriodic() {
int main() {
try {
testBonds();
testAngles();
testPeriodicTorsions();
testRBTorsions();
testCoulomb();
testLJ();
testExclusionsAnd14();
......
platforms/cuda/tests/TestCudaPeriodicTorsionForce.cpp 0 → 100644
View file @ ad75a390
/* -------------------------------------------------------------------------- *
* 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 Cuda implementation of PeriodicTorsionForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "CudaPlatform.h"
#include "PeriodicTorsionForce.h"
#include "System.h"
#include "LangevinIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testPeriodicTorsions() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
PeriodicTorsionForce* forceField = new PeriodicTorsionForce(1);
forceField->setTorsionParameters(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);
}
int main() {
try {
testPeriodicTorsions();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda/tests/TestCudaRBTorsionForce.cpp 0 → 100644
View file @ ad75a390
/* -------------------------------------------------------------------------- *
* 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 Cuda implementation of RBTorsionForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "CudaPlatform.h"
#include "RBTorsionForce.h"
#include "System.h"
#include "LangevinIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testRBTorsions() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
RBTorsionForce* forceField = new RBTorsionForce(1);
forceField->setTorsionParameters(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);
}
int main() {
try {
testRBTorsions();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/reference/src/ReferenceKernelFactory.cpp
View file @ ad75a390
......@@ -36,8 +36,18 @@
using namespace OpenMM;
KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Platform& platform, OpenMMContextImpl& context) const {
if (name == CalcStandardMMForceFieldKernel::Name())
return new ReferenceCalcStandardMMForceFieldKernel(name, platform);
if (name == CalcNonbondedForceKernel::Name())
return new ReferenceCalcNonbondedForceKernel(name, platform);
else if (name == CalcHarmonicBondForceKernel::Name())
return new ReferenceCalcHarmonicBondForceKernel(name, platform);
else if (name == CalcHarmonicAngleForceKernel::Name())
return new ReferenceCalcHarmonicAngleForceKernel(name, platform);
else if (name == CalcHarmonicAngleForceKernel::Name())
return new ReferenceCalcHarmonicAngleForceKernel(name, platform);
else if (name == CalcPeriodicTorsionForceKernel::Name())
return new ReferenceCalcPeriodicTorsionForceKernel(name, platform);
else if (name == CalcRBTorsionForceKernel::Name())
return new ReferenceCalcRBTorsionForceKernel(name, platform);
else if (name == CalcGBSAOBCForceFieldKernel::Name())
return new ReferenceCalcGBSAOBCForceFieldKernel(name, platform);
else if (name == IntegrateVerletStepKernel::Name())
......@@ -52,7 +62,7 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
return new ReferenceCalcKineticEnergyKernel(name, platform);
else if (name == RemoveCMMotionKernel::Name())
return new ReferenceRemoveCMMotionKernel(name, platform);
else {
throw OpenMMException( (std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str() );
}
else {
throw OpenMMException( (std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str() );
}
}
platforms/reference/src/ReferenceKernels.cpp
View file @ ad75a390
......@@ -29,9 +29,6 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "System.h"
#include "ReferenceKernels.h"
#include "ReferenceFloatStreamImpl.h"
#include "gbsa/CpuObc.h"
......@@ -107,41 +104,15 @@ void disposeRealArray(RealOpenMM** array, int size) {
}
}
ReferenceCalcStandardMMForceFieldKernel::~ReferenceCalcStandardMMForceFieldKernel() {
ReferenceCalcHarmonicBondForceKernel::~ReferenceCalcHarmonicBondForceKernel() {
disposeIntArray(bondIndexArray, numBonds);
disposeRealArray(bondParamArray, numBonds);
disposeIntArray(angleIndexArray, numAngles);
disposeRealArray(angleParamArray, numAngles);
disposeIntArray(periodicTorsionIndexArray, numPeriodicTorsions);
disposeRealArray(periodicTorsionParamArray, numPeriodicTorsions);
disposeIntArray(rbTorsionIndexArray, numRBTorsions);
disposeRealArray(rbTorsionParamArray, numRBTorsions);
disposeRealArray(atomParamArray, numAtoms);
disposeIntArray(exclusionArray, numAtoms);
disposeIntArray(bonded14IndexArray, num14);
disposeRealArray(bonded14ParamArray, num14);
if (neighborList != NULL)
delete neighborList;
}
void ReferenceCalcStandardMMForceFieldKernel::initialize(const System& system, const StandardMMForceField& force, const std::vector<std::set<int> >& exclusions) {
numAtoms = force.getNumAtoms();
void ReferenceCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
numBonds = force.getNumBonds();
numAngles = force.getNumAngles();
numPeriodicTorsions = force.getNumPeriodicTorsions();
numRBTorsions = force.getNumRBTorsions();
num14 = force.getNumNonbonded14();
bondIndexArray = allocateIntArray(numBonds, 2);
bondParamArray = allocateRealArray(numBonds, 2);
angleIndexArray = allocateIntArray(numAngles, 3);
angleParamArray = allocateRealArray(numAngles, 2);
periodicTorsionIndexArray = allocateIntArray(numPeriodicTorsions, 4);
periodicTorsionParamArray = allocateRealArray(numPeriodicTorsions, 3);
rbTorsionIndexArray = allocateIntArray(numRBTorsions, 4);
rbTorsionParamArray = allocateRealArray(numRBTorsions, 6);
bonded14IndexArray = allocateIntArray(num14, 2);
bonded14ParamArray = allocateRealArray(num14, 3);
atomParamArray = allocateRealArray(numAtoms, 3);
for (int i = 0; i < force.getNumBonds(); ++i) {
int atom1, atom2;
double length, k;
......@@ -151,6 +122,40 @@ void ReferenceCalcStandardMMForceFieldKernel::initialize(const System& system, c
bondParamArray[i][0] = length;
bondParamArray[i][1] = k;
}
}
void ReferenceCalcHarmonicBondForceKernel::executeForces(OpenMMContextImpl& context) {
RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
ReferenceBondForce refBondForce;
ReferenceHarmonicBondIxn harmonicBond;
refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
}
double ReferenceCalcHarmonicBondForceKernel::executeEnergy(OpenMMContextImpl& context) {
RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumAtoms(), 3);
RealOpenMM* energyArray = new RealOpenMM[numBonds];
RealOpenMM energy = 0;
ReferenceBondForce refBondForce;
ReferenceHarmonicBondIxn harmonicBond;
for (int i = 0; i < numBonds; ++i)
energyArray[i] = 0;
refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, energyArray, 0, &energy, harmonicBond);
disposeRealArray(forceData, context.getSystem().getNumAtoms());
delete[] energyArray;
return energy;
}
ReferenceCalcHarmonicAngleForceKernel::~ReferenceCalcHarmonicAngleForceKernel() {
disposeIntArray(angleIndexArray, numAngles);
disposeRealArray(angleParamArray, numAngles);
}
void ReferenceCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
numAngles = force.getNumAngles();
angleIndexArray = allocateIntArray(numAngles, 3);
angleParamArray = allocateRealArray(numAngles, 2);
for (int i = 0; i < force.getNumAngles(); ++i) {
int atom1, atom2, atom3;
double angle, k;
......@@ -161,33 +166,141 @@ void ReferenceCalcStandardMMForceFieldKernel::initialize(const System& system, c
angleParamArray[i][0] = angle;
angleParamArray[i][1] = k;
}
for (int i = 0; i < force.getNumPeriodicTorsions(); ++i) {
}
void ReferenceCalcHarmonicAngleForceKernel::executeForces(OpenMMContextImpl& context) {
RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
ReferenceBondForce refBondForce;
ReferenceAngleBondIxn angleBond;
refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, angleBond);
}
double ReferenceCalcHarmonicAngleForceKernel::executeEnergy(OpenMMContextImpl& context) {
RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumAtoms(), 3);
RealOpenMM* energyArray = new RealOpenMM[numAngles];
RealOpenMM energy = 0;
ReferenceBondForce refBondForce;
ReferenceAngleBondIxn angleBond;
for (int i = 0; i < numAngles; ++i)
energyArray[i] = 0;
refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, angleBond);
disposeRealArray(forceData, context.getSystem().getNumAtoms());
delete[] energyArray;
return energy;
}
ReferenceCalcPeriodicTorsionForceKernel::~ReferenceCalcPeriodicTorsionForceKernel() {
disposeIntArray(torsionIndexArray, numTorsions);
disposeRealArray(torsionParamArray, numTorsions);
}
void ReferenceCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
numTorsions = force.getNumTorsions();
torsionIndexArray = allocateIntArray(numTorsions, 4);
torsionParamArray = allocateRealArray(numTorsions, 3);
for (int i = 0; i < force.getNumTorsions(); ++i) {
int atom1, atom2, atom3, atom4, periodicity;
double phase, k;
force.getPeriodicTorsionParameters(i, atom1, atom2, atom3, atom4, periodicity, phase, k);
periodicTorsionIndexArray[i][0] = atom1;
periodicTorsionIndexArray[i][1] = atom2;
periodicTorsionIndexArray[i][2] = atom3;
periodicTorsionIndexArray[i][3] = atom4;
periodicTorsionParamArray[i][0] = k;
periodicTorsionParamArray[i][1] = phase;
periodicTorsionParamArray[i][2] = periodicity;
force.getTorsionParameters(i, atom1, atom2, atom3, atom4, periodicity, phase, k);
torsionIndexArray[i][0] = atom1;
torsionIndexArray[i][1] = atom2;
torsionIndexArray[i][2] = atom3;
torsionIndexArray[i][3] = atom4;
torsionParamArray[i][0] = k;
torsionParamArray[i][1] = phase;
torsionParamArray[i][2] = periodicity;
}
for (int i = 0; i < force.getNumRBTorsions(); ++i) {
}
void ReferenceCalcPeriodicTorsionForceKernel::executeForces(OpenMMContextImpl& context) {
RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
ReferenceBondForce refBondForce;
ReferenceProperDihedralBond periodicTorsionBond;
refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, periodicTorsionBond);
}
double ReferenceCalcPeriodicTorsionForceKernel::executeEnergy(OpenMMContextImpl& context) {
RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumAtoms(), 3);
RealOpenMM* energyArray = new RealOpenMM[numTorsions];
RealOpenMM energy = 0;
ReferenceBondForce refBondForce;
ReferenceProperDihedralBond periodicTorsionBond;
for (int i = 0; i < numTorsions; ++i)
energyArray[i] = 0;
refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, energyArray, 0, &energy, periodicTorsionBond);
disposeRealArray(forceData, context.getSystem().getNumAtoms());
delete[] energyArray;
return energy;
}
ReferenceCalcRBTorsionForceKernel::~ReferenceCalcRBTorsionForceKernel() {
disposeIntArray(torsionIndexArray, numTorsions);
disposeRealArray(torsionParamArray, numTorsions);
}
void ReferenceCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
numTorsions = force.getNumTorsions();
torsionIndexArray = allocateIntArray(numTorsions, 4);
torsionParamArray = allocateRealArray(numTorsions, 6);
for (int i = 0; i < force.getNumTorsions(); ++i) {
int atom1, atom2, atom3, atom4;
double c0, c1, c2, c3, c4, c5;
force.getRBTorsionParameters(i, atom1, atom2, atom3, atom4, c0, c1, c2, c3, c4, c5);
rbTorsionIndexArray[i][0] = atom1;
rbTorsionIndexArray[i][1] = atom2;
rbTorsionIndexArray[i][2] = atom3;
rbTorsionIndexArray[i][3] = atom4;
rbTorsionParamArray[i][0] = c0;
rbTorsionParamArray[i][1] = c1;
rbTorsionParamArray[i][2] = c2;
rbTorsionParamArray[i][3] = c3;
rbTorsionParamArray[i][4] = c4;
rbTorsionParamArray[i][5] = c5;
force.getTorsionParameters(i, atom1, atom2, atom3, atom4, c0, c1, c2, c3, c4, c5);
torsionIndexArray[i][0] = atom1;
torsionIndexArray[i][1] = atom2;
torsionIndexArray[i][2] = atom3;
torsionIndexArray[i][3] = atom4;
torsionParamArray[i][0] = c0;
torsionParamArray[i][1] = c1;
torsionParamArray[i][2] = c2;
torsionParamArray[i][3] = c3;
torsionParamArray[i][4] = c4;
torsionParamArray[i][5] = c5;
}
}
void ReferenceCalcRBTorsionForceKernel::executeForces(OpenMMContextImpl& context) {
RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
ReferenceBondForce refBondForce;
ReferenceRbDihedralBond rbTorsionBond;
refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, rbTorsionBond);
}
double ReferenceCalcRBTorsionForceKernel::executeEnergy(OpenMMContextImpl& context) {
RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumAtoms(), 3);
RealOpenMM* energyArray = new RealOpenMM[numTorsions];
RealOpenMM energy = 0;
ReferenceBondForce refBondForce;
ReferenceRbDihedralBond rbTorsionBond;
for (int i = 0; i < numTorsions; ++i)
energyArray[i] = 0;
refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, energyArray, 0, &energy, rbTorsionBond);
disposeRealArray(forceData, context.getSystem().getNumAtoms());
delete[] energyArray;
return energy;
}
ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
disposeRealArray(atomParamArray, numAtoms);
disposeIntArray(exclusionArray, numAtoms);
disposeIntArray(bonded14IndexArray, num14);
disposeRealArray(bonded14ParamArray, num14);
if (neighborList != NULL)
delete neighborList;
}
void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force, const std::vector<std::set<int> >& exclusions) {
numAtoms = force.getNumAtoms();
num14 = force.getNumNonbonded14();
bonded14IndexArray = allocateIntArray(num14, 2);
bonded14ParamArray = allocateRealArray(num14, 3);
atomParamArray = allocateRealArray(numAtoms, 3);
RealOpenMM sqrtEps = static_cast<RealOpenMM>( std::sqrt(138.935485) );
for (int i = 0; i < numAtoms; ++i) {
double charge, radius, depth;
......@@ -215,7 +328,7 @@ void ReferenceCalcStandardMMForceFieldKernel::initialize(const System& system, c
bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge*sqrtEps*sqrtEps);
}
nonbondedMethod = CalcStandardMMForceFieldKernel::NonbondedMethod(force.getNonbondedMethod());
nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
Vec3 boxVectors[3];
force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
......@@ -229,18 +342,9 @@ void ReferenceCalcStandardMMForceFieldKernel::initialize(const System& system, c
}
void ReferenceCalcStandardMMForceFieldKernel::executeForces(OpenMMContextImpl& context) {
void ReferenceCalcNonbondedForceKernel::executeForces(OpenMMContextImpl& context) {
RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
ReferenceBondForce refBondForce;
ReferenceHarmonicBondIxn harmonicBond;
refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
ReferenceAngleBondIxn angleBond;
refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, angleBond);
ReferenceProperDihedralBond periodicTorsionBond;
refBondForce.calculateForce(numPeriodicTorsions, periodicTorsionIndexArray, posData, periodicTorsionParamArray, forceData, 0, 0, 0, periodicTorsionBond);
ReferenceRbDihedralBond rbTorsionBond;
refBondForce.calculateForce(numRBTorsions, rbTorsionIndexArray, posData, rbTorsionParamArray, forceData, 0, 0, 0, rbTorsionBond);
ReferenceLJCoulombIxn clj;
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (nonbondedMethod != NoCutoff) {
......@@ -250,35 +354,17 @@ void ReferenceCalcStandardMMForceFieldKernel::executeForces(OpenMMContextImpl& c
if (periodic)
clj.setPeriodic(periodicBoxSize);
clj.calculatePairIxn(numAtoms, posData, atomParamArray, exclusionArray, 0, forceData, 0, 0);
ReferenceBondForce refBondForce;
ReferenceLJCoulomb14 nonbonded14;
if (nonbondedMethod != NoCutoff)
nonbonded14.setUseCutoff(nonbondedCutoff, 78.3f);
refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
}
double ReferenceCalcStandardMMForceFieldKernel::executeEnergy(OpenMMContextImpl& context) {
double ReferenceCalcNonbondedForceKernel::executeEnergy(OpenMMContextImpl& context) {
RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
RealOpenMM** forceData = allocateRealArray(numAtoms, 3);
int arraySize = max(max(max(max(numAtoms, numBonds), numAngles), numPeriodicTorsions), numRBTorsions);
RealOpenMM* energyArray = new RealOpenMM[arraySize];
RealOpenMM energy = 0;
ReferenceBondForce refBondForce;
ReferenceHarmonicBondIxn harmonicBond;
for (int i = 0; i < arraySize; ++i)
energyArray[i] = 0;
refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, energyArray, 0, &energy, harmonicBond);
ReferenceAngleBondIxn angleBond;
for (int i = 0; i < arraySize; ++i)
energyArray[i] = 0;
refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, angleBond);
ReferenceProperDihedralBond periodicTorsionBond;
for (int i = 0; i < arraySize; ++i)
energyArray[i] = 0;
refBondForce.calculateForce(numPeriodicTorsions, periodicTorsionIndexArray, posData, periodicTorsionParamArray, forceData, energyArray, 0, &energy, periodicTorsionBond);
ReferenceRbDihedralBond rbTorsionBond;
for (int i = 0; i < arraySize; ++i)
energyArray[i] = 0;
refBondForce.calculateForce(numRBTorsions, rbTorsionIndexArray, posData, rbTorsionParamArray, forceData, energyArray, 0, &energy, rbTorsionBond);
ReferenceLJCoulombIxn clj;
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (nonbondedMethod != NoCutoff) {
......@@ -288,10 +374,12 @@ double ReferenceCalcStandardMMForceFieldKernel::executeEnergy(OpenMMContextImpl&
if (periodic)
clj.setPeriodic(periodicBoxSize);
clj.calculatePairIxn(numAtoms, posData, atomParamArray, exclusionArray, 0, forceData, 0, &energy);
ReferenceBondForce refBondForce;
ReferenceLJCoulomb14 nonbonded14;
if (nonbondedMethod != NoCutoff)
nonbonded14.setUseCutoff(nonbondedCutoff, 78.3f);
for (int i = 0; i < arraySize; ++i)
RealOpenMM* energyArray = new RealOpenMM[numAtoms];
for (int i = 0; i < numAtoms; ++i)
energyArray[i] = 0;
refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);
disposeRealArray(forceData, numAtoms);
......
platforms/reference/src/ReferenceKernels.h
View file @ ad75a390
......@@ -46,23 +46,159 @@ class ReferenceVerletDynamics;
namespace OpenMM {
/**
* This kernel is invoked by StandardMMForceField to calculate the forces acting on the system.
* This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
*/
class ReferenceCalcStandardMMForceFieldKernel : public CalcStandardMMForceFieldKernel {
class ReferenceCalcHarmonicBondForceKernel : public CalcHarmonicBondForceKernel {
public:
ReferenceCalcStandardMMForceFieldKernel(std::string name, const Platform& platform) : CalcStandardMMForceFieldKernel(name, platform) {
ReferenceCalcHarmonicBondForceKernel(std::string name, const Platform& platform) : CalcHarmonicBondForceKernel(name, platform) {
}
~ReferenceCalcStandardMMForceFieldKernel();
~ReferenceCalcHarmonicBondForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the StandardMMForceField this kernel will be used for
* @param force the HarmonicBondForce this kernel will be used for
*/
void initialize(const System& system, const HarmonicBondForce& force);
/**
* Execute the kernel to calculate the forces.
*
* @param context the context in which to execute this kernel
*/
void executeForces(OpenMMContextImpl& context);
/**
* Execute the kernel to calculate the energy.
*
* @param context the context in which to execute this kernel
* @return the potential energy due to the HarmonicBondForce
*/
double executeEnergy(OpenMMContextImpl& context);
private:
int numBonds;
int **bondIndexArray;
RealOpenMM **bondParamArray;
};
/**
* This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
*/
class ReferenceCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
public:
ReferenceCalcHarmonicAngleForceKernel(std::string name, const Platform& platform) : CalcHarmonicAngleForceKernel(name, platform) {
}
~ReferenceCalcHarmonicAngleForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the HarmonicAngleForce this kernel will be used for
*/
void initialize(const System& system, const HarmonicAngleForce& force);
/**
* Execute the kernel to calculate the forces.
*
* @param context the context in which to execute this kernel
*/
void executeForces(OpenMMContextImpl& context);
/**
* Execute the kernel to calculate the energy.
*
* @param context the context in which to execute this kernel
* @return the potential energy due to the HarmonicAngleForce
*/
double executeEnergy(OpenMMContextImpl& context);
private:
int numAngles;
int **angleIndexArray;
RealOpenMM **angleParamArray;
};
/**
* This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
class ReferenceCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
public:
ReferenceCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform) : CalcPeriodicTorsionForceKernel(name, platform) {
}
~ReferenceCalcPeriodicTorsionForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the PeriodicTorsionForce this kernel will be used for
*/
void initialize(const System& system, const PeriodicTorsionForce& force);
/**
* Execute the kernel to calculate the forces.
*
* @param context the context in which to execute this kernel
*/
void executeForces(OpenMMContextImpl& context);
/**
* Execute the kernel to calculate the energy.
*
* @param context the context in which to execute this kernel
* @return the potential energy due to the PeriodicTorsionForce
*/
double executeEnergy(OpenMMContextImpl& context);
private:
int numTorsions;
int **torsionIndexArray;
RealOpenMM **torsionParamArray;
};
/**
* This kernel is invoked by RBTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
class ReferenceCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
public:
ReferenceCalcRBTorsionForceKernel(std::string name, const Platform& platform) : CalcRBTorsionForceKernel(name, platform) {
}
~ReferenceCalcRBTorsionForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the RBTorsionForce this kernel will be used for
*/
void initialize(const System& system, const RBTorsionForce& force);
/**
* Execute the kernel to calculate the forces.
*
* @param context the context in which to execute this kernel
*/
void executeForces(OpenMMContextImpl& context);
/**
* Execute the kernel to calculate the energy.
*
* @param context the context in which to execute this kernel
* @return the potential energy due to the RBTorsionForce
*/
double executeEnergy(OpenMMContextImpl& context);
private:
int numTorsions;
int **torsionIndexArray;
RealOpenMM **torsionParamArray;
};
/**
* This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
*/
class ReferenceCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
public:
ReferenceCalcNonbondedForceKernel(std::string name, const Platform& platform) : CalcNonbondedForceKernel(name, platform) {
}
~ReferenceCalcNonbondedForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the NonbondedForce this kernel will be used for
* @param exclusions the i'th element lists the indices of all atoms with which the i'th atom should not interact through
* nonbonded forces. Bonded 1-4 pairs are also included in this list, since they should be omitted from
* the standard nonbonded calculation.
*/
void initialize(const System& system, const StandardMMForceField& force, const std::vector<std::set<int> >& exclusions);
void initialize(const System& system, const NonbondedForce& force, const std::vector<std::set<int> >& exclusions);
/**
* Execute the kernel to calculate the forces.
*
......@@ -73,13 +209,13 @@ public:
* Execute the kernel to calculate the energy.
*
* @param context the context in which to execute this kernel
* @return the potential energy due to the StandardMMForceField
* @return the potential energy due to the NonbondedForce
*/
double executeEnergy(OpenMMContextImpl& context);
private:
int numAtoms, numBonds, numAngles, numPeriodicTorsions, numRBTorsions, num14;
int **bondIndexArray, **angleIndexArray, **periodicTorsionIndexArray, **rbTorsionIndexArray, **exclusionArray, **bonded14IndexArray;
RealOpenMM **bondParamArray, **angleParamArray, **periodicTorsionParamArray, **rbTorsionParamArray, **atomParamArray, **bonded14ParamArray;
int numAtoms, num14;
int **exclusionArray, **bonded14IndexArray;
RealOpenMM **atomParamArray, **bonded14ParamArray;
RealOpenMM nonbondedCutoff, periodicBoxSize[3];
std::vector<std::set<int> > exclusions;
NonbondedMethod nonbondedMethod;
......
platforms/reference/src/ReferencePlatform.cpp
View file @ ad75a390
......@@ -38,7 +38,11 @@ using namespace OpenMM;
ReferencePlatform::ReferencePlatform() {
ReferenceKernelFactory* factory = new ReferenceKernelFactory();
registerKernelFactory(CalcStandardMMForceFieldKernel::Name(), factory);
registerKernelFactory(CalcHarmonicBondForceKernel::Name(), factory);
registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);
registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
registerKernelFactory(CalcGBSAOBCForceFieldKernel::Name(), factory);
registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
......
platforms/reference/tests/TestReferenceAndersenThermostat.cpp
View file @ ad75a390
......@@ -37,7 +37,7 @@
#include "AndersenThermostat.h"
#include "OpenMMContext.h"
#include "ReferencePlatform.h"
#include "StandardMMForceField.h"
#include "NonbondedForce.h"
#include "System.h"
#include "VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
......@@ -55,7 +55,7 @@ void testTemperature() {
ReferencePlatform platform;
System system(numAtoms, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0, 0);
NonbondedForce* forceField = new NonbondedForce(numAtoms, 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);
......
platforms/reference/tests/TestReferenceBrownianIntegrator.cpp
View file @ ad75a390
......@@ -36,7 +36,8 @@
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "ReferencePlatform.h"
#include "StandardMMForceField.h"
#include "HarmonicBondForce.h"
#include "NonbondedForce.h"
#include "System.h"
#include "BrownianIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
......@@ -56,7 +57,7 @@ void testSingleBond() {
system.setAtomMass(1, 2.0);
double dt = 0.01;
BrownianIntegrator integrator(0, 0.1, dt);
StandardMMForceField* forceField = new StandardMMForceField(2, 1, 0, 0, 0, 0);
HarmonicBondForce* forceField = new HarmonicBondForce(1);
forceField->setBondParameters(0, 0, 1, 1.5, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
......@@ -90,7 +91,7 @@ void testTemperature() {
ReferencePlatform platform;
System system(numAtoms, 0);
BrownianIntegrator integrator(temp, 2.0, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, numBonds, 0, 0, 0, 0);
HarmonicBondForce* forceField = new HarmonicBondForce(numBonds);
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);
......@@ -127,7 +128,7 @@ void testConstraints() {
ReferencePlatform platform;
System system(numAtoms, numAtoms-1);
BrownianIntegrator integrator(temp, 2.0, 0.001);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0, 0);
NonbondedForce* forceField = new NonbondedForce(numAtoms, 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);
......
platforms/reference/tests/TestReferenceCMMotionRemover.cpp
View file @ ad75a390
......@@ -37,7 +37,8 @@
#include "CMMotionRemover.h"
#include "OpenMMContext.h"
#include "ReferencePlatform.h"
#include "StandardMMForceField.h"
#include "HarmonicBondForce.h"
#include "NonbondedForce.h"
#include "System.h"
#include "VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
......@@ -65,13 +66,15 @@ void testMotionRemoval() {
ReferencePlatform platform;
System system(numAtoms, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 1, 0, 0, 0, 0);
HarmonicBondForce* bonds = new HarmonicBondForce(1);
bonds->setBondParameters(0, 2, 3, 2.0, 0.5);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(numAtoms, 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);
nonbonded->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);
system.addForce(nonbonded);
CMMotionRemover* remover = new CMMotionRemover();
system.addForce(remover);
OpenMMContext context(system, integrator, platform);
......
platforms/reference/tests/TestReferenceHarmonicAngleForce.cpp 0 → 100644
View file @ ad75a390
/* -------------------------------------------------------------------------- *
* 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 reference implementation of HarmonicAngleForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "ReferencePlatform.h"
#include "HarmonicAngleForce.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 testAngles() {
ReferencePlatform platform;
System system(4, 0);
VerletIntegrator integrator(0.01);
HarmonicAngleForce* forceField = new HarmonicAngleForce(2);
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);
}
int main() {
try {
testAngles();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/reference/tests/TestReferenceHarmonicBondForce.cpp 0 → 100644
View file @ ad75a390
/* -------------------------------------------------------------------------- *
* 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 reference implementation of HarmonicBondForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "ReferencePlatform.h"
#include "HarmonicBondForce.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() {
ReferencePlatform platform;
System system(3, 0);
VerletIntegrator integrator(0.01);
HarmonicBondForce* forceField = new HarmonicBondForce(2);
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);
}
int main() {
try {
testBonds();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/reference/tests/TestReferenceLangevinIntegrator.cpp
View file @ ad75a390
......@@ -36,7 +36,8 @@
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "ReferencePlatform.h"
#include "StandardMMForceField.h"
#include "HarmonicBondForce.h"
#include "NonbondedForce.h"
#include "System.h"
#include "LangevinIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
......@@ -55,7 +56,7 @@ void testSingleBond() {
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, 0);
HarmonicBondForce* forceField = new HarmonicBondForce(1);
forceField->setBondParameters(0, 0, 1, 1.5, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
......@@ -99,7 +100,7 @@ void testTemperature() {
ReferencePlatform platform;
System system(numAtoms, 0);
LangevinIntegrator integrator(temp, 2.0, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0, 0);
NonbondedForce* forceField = new NonbondedForce(numAtoms, 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);
......@@ -134,7 +135,7 @@ void testConstraints() {
ReferencePlatform platform;
System system(numAtoms, numAtoms-1);
LangevinIntegrator integrator(temp, 2.0, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0, 0);
NonbondedForce* forceField = new NonbondedForce(numAtoms, 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);
......
platforms/reference/tests/TestReferenceStandardMMForceField.cpp → platforms/reference/tests/TestReferenceNonbondedForce.cpp
View file @ ad75a390
......@@ -30,16 +30,17 @@
* -------------------------------------------------------------------------- */
/**
* This tests all the different force terms in the reference implementation of StandardMMForceField.
* This tests all the different force terms in the reference implementation of NonbondedForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "ReferencePlatform.h"
#include "StandardMMForceField.h"
#include "NonbondedForce.h"
#include "System.h"
#include "VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "HarmonicBondForce.h"
#include <iostream>
#include <vector>
......@@ -48,114 +49,11 @@ using namespace std;
const double TOL = 1e-5;
void testBonds() {
ReferencePlatform platform;
System system(3, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 2, 0, 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() {
ReferencePlatform platform;
System system(4, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 2, 0, 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() {
ReferencePlatform platform;
System system(4, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 1, 0, 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() {
ReferencePlatform platform;
System system(4, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 0, 1, 0);
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() {
ReferencePlatform platform;
System system(2, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0, 0);
NonbondedForce* forceField = new NonbondedForce(2, 0);
forceField->setAtomParameters(0, 0.5, 1, 0);
forceField->setAtomParameters(1, -1.5, 1, 0);
system.addForce(forceField);
......@@ -176,7 +74,7 @@ void testLJ() {
ReferencePlatform platform;
System system(2, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0, 0);
NonbondedForce* forceField = new NonbondedForce(2, 0);
forceField->setAtomParameters(0, 0, 1.2, 1);
forceField->setAtomParameters(1, 0, 1.4, 2);
system.addForce(forceField);
......@@ -199,12 +97,14 @@ void testExclusionsAnd14() {
ReferencePlatform platform;
System system(5, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0, 2);
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);
HarmonicBondForce* bonds = new HarmonicBondForce(4);
bonds->setBondParameters(0, 0, 1, 1, 0);
bonds->setBondParameters(1, 1, 2, 1, 0);
bonds->setBondParameters(2, 2, 3, 1, 0);
bonds->setBondParameters(3, 3, 4, 1, 0);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(5, 2);
system.addForce(nonbonded);
OpenMMContext context(system, integrator, platform);
for (int i = 1; i < 5; ++i) {
......@@ -213,13 +113,13 @@ void testExclusionsAnd14() {
vector<Vec3> positions(5);
const double r = 1.0;
for (int j = 0; j < 5; ++j) {
forceField->setAtomParameters(j, 0, 1.5, 0);
nonbonded->setAtomParameters(j, 0, 1.5, 0);
positions[j] = Vec3(0, j, 0);
}
forceField->setAtomParameters(0, 0, 1.5, 1);
forceField->setAtomParameters(i, 0, 1.5, 1);
forceField->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
nonbonded->setAtomParameters(0, 0, 1.5, 1);
nonbonded->setAtomParameters(i, 0, 1.5, 1);
nonbonded->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
positions[i] = Vec3(r, 0, 0);
context.reinitialize();
context.setPositions(positions);
......@@ -243,10 +143,10 @@ void testExclusionsAnd14() {
// Test Coulomb forces
forceField->setAtomParameters(0, 2, 1.5, 0);
forceField->setAtomParameters(i, 2, 1.5, 0);
forceField->setNonbonded14Parameters(0, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
nonbonded->setAtomParameters(0, 2, 1.5, 0);
nonbonded->setAtomParameters(i, 2, 1.5, 0);
nonbonded->setNonbonded14Parameters(0, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
......@@ -271,11 +171,11 @@ void testCutoff() {
ReferencePlatform platform;
System system(3, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 0, 0, 0, 0, 0);
NonbondedForce* forceField = new NonbondedForce(3, 0);
forceField->setAtomParameters(0, 1.0, 1, 0);
forceField->setAtomParameters(1, 1.0, 1, 0);
forceField->setAtomParameters(2, 1.0, 1, 0);
forceField->setNonbondedMethod(StandardMMForceField::CutoffNonPeriodic);
forceField->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
const double cutoff = 2.9;
forceField->setCutoffDistance(cutoff);
system.addForce(forceField);
......@@ -304,15 +204,17 @@ void testCutoff14() {
ReferencePlatform platform;
System system(5, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0, 2);
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);
forceField->setNonbondedMethod(StandardMMForceField::CutoffNonPeriodic);
HarmonicBondForce* bonds = new HarmonicBondForce(4);
bonds->setBondParameters(0, 0, 1, 1, 0);
bonds->setBondParameters(1, 1, 2, 1, 0);
bonds->setBondParameters(2, 2, 3, 1, 0);
bonds->setBondParameters(3, 3, 4, 1, 0);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(5, 2);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
const double cutoff = 3.5;
forceField->setCutoffDistance(cutoff);
system.addForce(forceField);
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(5);
positions[0] = Vec3(0, 0, 0);
......@@ -324,12 +226,12 @@ void testCutoff14() {
// Test LJ forces
forceField->setAtomParameters(0, 0, 1.5, 1);
nonbonded->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);
forceField->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
nonbonded->setAtomParameters(j, 0, 1.5, 0);
nonbonded->setAtomParameters(i, 0, 1.5, 1);
nonbonded->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
context.reinitialize();
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
......@@ -354,10 +256,10 @@ void testCutoff14() {
// Test Coulomb forces
const double q = 0.7;
forceField->setAtomParameters(0, q, 1.5, 0);
forceField->setAtomParameters(i, q, 1.5, 0);
forceField->setNonbonded14Parameters(0, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
nonbonded->setAtomParameters(0, q, 1.5, 0);
nonbonded->setAtomParameters(i, q, 1.5, 0);
nonbonded->setNonbonded14Parameters(0, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
......@@ -385,16 +287,18 @@ void testPeriodic() {
ReferencePlatform platform;
System system(3, 0);
VerletIntegrator integrator(0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 1, 0, 0, 0, 0);
forceField->setAtomParameters(0, 1.0, 1, 0);
forceField->setAtomParameters(1, 1.0, 1, 0);
forceField->setAtomParameters(2, 1.0, 1, 0);
forceField->setBondParameters(0, 0, 1, 1.0, 0.0);
forceField->setNonbondedMethod(StandardMMForceField::CutoffPeriodic);
HarmonicBondForce* bonds = new HarmonicBondForce(1);
bonds->setBondParameters(0, 0, 1, 1, 0);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(3, 0);
nonbonded->setAtomParameters(0, 1.0, 1, 0);
nonbonded->setAtomParameters(1, 1.0, 1, 0);
nonbonded->setAtomParameters(2, 1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 2.0;
forceField->setCutoffDistance(cutoff);
forceField->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
system.addForce(forceField);
nonbonded->setCutoffDistance(cutoff);
nonbonded->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
system.addForce(nonbonded);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
......@@ -415,10 +319,6 @@ void testPeriodic() {
int main() {
try {
testBonds();
testAngles();
testPeriodicTorsions();
testRBTorsions();
testCoulomb();
testLJ();
testExclusionsAnd14();
......
platforms/reference/tests/TestReferencePeriodicTorsionForce.cpp 0 → 100644
View file @ ad75a390
/* -------------------------------------------------------------------------- *
* 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 reference implementation of PeriodicTorsionForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "OpenMMContext.h"
#include "ReferencePlatform.h"
#include "PeriodicTorsionForce.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 testPeriodicTorsions() {
ReferencePlatform platform;
System system(4, 0);
VerletIntegrator integrator(0.01);
PeriodicTorsionForce* forceField = new PeriodicTorsionForce(1);
forceField->setTorsionParameters(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);
}
int main() {
try {
testPeriodicTorsions();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/reference/tests/TestReferenceRBTorsionForce.cpp 0 → 100644
View file @ ad75a390
This diff is collapsed.
tests/TestFindExclusions.cpp
View file @ ad75a390
This diff is collapsed.
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