Skip to content

GitLab

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

Continuing to refactor tests

parent a5a52dd1
Expand all Show whitespace changes
Inline Side-by-side
Showing with 1839 additions and 1928 deletions
platforms/reference/tests/TestReferenceCustomAngleForce.cpp
View file @ ccd811da
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2010 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,122 +29,9 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of CustomAngleForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomAngleForce.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
#include "ReferenceTests.h"
#include "TestCustomAngleForce.h"
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testAngles() {
// Create a system using a CustomAngleForce.
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomAngleForce* custom = new CustomAngleForce("scale*k*(theta-theta0)^2");
custom->addPerAngleParameter("theta0");
custom->addPerAngleParameter("k");
custom->addGlobalParameter("scale", 0.5);
vector<double> parameters(2);
parameters[0] = 1.5;
parameters[1] = 0.8;
custom->addAngle(0, 1, 2, parameters);
parameters[0] = 2.0;
parameters[1] = 0.5;
custom->addAngle(1, 2, 3, parameters);
customSystem.addForce(custom);
ASSERT(!custom->usesPeriodicBoundaryConditions());
ASSERT(!customSystem.usesPeriodicBoundaryConditions());
// Create an identical system using a HarmonicAngleForce.
System harmonicSystem;
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
HarmonicAngleForce* harmonic = new HarmonicAngleForce();
harmonic->addAngle(0, 1, 2, 1.5, 0.8);
harmonic->addAngle(1, 2, 3, 2.0, 0.5);
harmonicSystem.addForce(harmonic);
// Set the atoms in various positions, and verify that both systems give identical forces and energy.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(4);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context c1(customSystem, integrator1, platform);
Context c2(harmonicSystem, integrator2, platform);
for (int i = 0; i < 10; i++) {
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = s1.getForces();
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
// Try changing the angle parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[1] = 0.9;
custom->setAngleParameters(0, 0, 1, 2, parameters);
parameters[0] = 2.1;
parameters[1] = 0.6;
custom->setAngleParameters(1, 1, 2, 3, parameters);
custom->updateParametersInContext(c1);
harmonic->setAngleParameters(0, 0, 1, 2, 1.6, 0.9);
harmonic->setAngleParameters(1, 1, 2, 3, 2.1, 0.6);
harmonic->updateParametersInContext(c2);
{
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = s1.getForces();
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
void runPlatformTests() {
}
int main() {
try {
testAngles();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/reference/tests/TestReferenceCustomBondForce.cpp
View file @ ccd811da
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,90 +29,9 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of CustomBondForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomBondForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
#include "ReferenceTests.h"
#include "TestCustomBondForce.h"
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testBonds() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomBondForce* forceField = new CustomBondForce("scale*k*(r-r0)^2");
forceField->addPerBondParameter("r0");
forceField->addPerBondParameter("k");
forceField->addGlobalParameter("scale", 0.5);
vector<double> parameters(2);
parameters[0] = 1.5;
parameters[1] = 0.8;
forceField->addBond(0, 1, parameters);
parameters[0] = 1.2;
parameters[1] = 0.7;
forceField->addBond(1, 2, parameters);
system.addForce(forceField);
ASSERT(!forceField->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context 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);
}
// Try changing the bond parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[1] = 0.9;
forceField->setBondParameters(0, 0, 1, parameters);
parameters[0] = 1.3;
parameters[1] = 0.8;
forceField->setBondParameters(1, 1, 2, parameters);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, -0.9*0.4, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.8*0.3, 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.9*0.4*0.4 + 0.5*0.8*0.3*0.3, state.getPotentialEnergy(), TOL);
}
}
int main() {
try {
testBonds();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/reference/tests/TestReferenceCustomCentroidBondForce.cpp
View file @ ccd811da
......@@ -29,246 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of CustomCompoundBondForce.
*/
#include "ReferenceTests.h"
#include "TestCustomCentroidBondForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomCentroidBondForce.h"
#include "openmm/CustomCompoundBondForce.h"
#include "openmm/System.h"
#include "openmm/TabulatedFunction.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testHarmonicBond() {
System system;
system.addParticle(1.0);
system.addParticle(2.0);
system.addParticle(3.0);
system.addParticle(4.0);
system.addParticle(5.0);
CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "k*distance(g1,g2)^2");
force->addPerBondParameter("k");
vector<int> particles1;
particles1.push_back(0);
particles1.push_back(1);
vector<int> particles2;
particles2.push_back(2);
particles2.push_back(3);
particles2.push_back(4);
force->addGroup(particles1);
force->addGroup(particles2);
vector<int> groups;
groups.push_back(0);
groups.push_back(1);
vector<double> parameters;
parameters.push_back(1.0);
force->addBond(groups, parameters);
system.addForce(force);
ASSERT(!system.usesPeriodicBoundaryConditions());
// The center of mass of group 0 is (1.5, 0, 0).
vector<Vec3> positions(5);
positions[0] = Vec3(2.5, 0, 0);
positions[1] = Vec3(1, 0, 0);
// The center of mass of group 1 is (-1, 0, 0).
positions[2] = Vec3(-6, 0, 0);
positions[3] = Vec3(-1, 0, 0);
positions[4] = Vec3(2, 0, 0);
// Check the forces and energy.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(2.5*2.5, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(Vec3(-2*2.5*(1.0/3.0), 0, 0), state.getForces()[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-2*2.5*(2.0/3.0), 0, 0), state.getForces()[1], TOL);
ASSERT_EQUAL_VEC(Vec3(2*2.5*(3.0/12.0), 0, 0), state.getForces()[2], TOL);
ASSERT_EQUAL_VEC(Vec3(2*2.5*(4.0/12.0), 0, 0), state.getForces()[3], TOL);
ASSERT_EQUAL_VEC(Vec3(2*2.5*(5.0/12.0), 0, 0), state.getForces()[4], TOL);
// Update the per-bond parameter and see if the results change.
parameters[0] = 2.0;
force->setBondParameters(0, groups, parameters);
force->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(2*2.5*2.5, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(Vec3(-4*2.5*(1.0/3.0), 0, 0), state.getForces()[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-4*2.5*(2.0/3.0), 0, 0), state.getForces()[1], TOL);
ASSERT_EQUAL_VEC(Vec3(4*2.5*(3.0/12.0), 0, 0), state.getForces()[2], TOL);
ASSERT_EQUAL_VEC(Vec3(4*2.5*(4.0/12.0), 0, 0), state.getForces()[3], TOL);
ASSERT_EQUAL_VEC(Vec3(4*2.5*(5.0/12.0), 0, 0), state.getForces()[4], TOL);
// All the particles should be treated as a single molecule.
vector<std::vector<int> > molecules = context.getMolecules();
ASSERT_EQUAL(1, molecules.size());
ASSERT_EQUAL(5, molecules[0].size());
}
void testComplexFunction() {
int numParticles = 5;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(2.0);
vector<double> table(20);
for (int i = 0; i < 20; i++)
table[i] = sin(0.11*i);
// When every group contains only one particle, a CustomCentroidBondForce is identical to a
// CustomCompoundBondForce. Use that to test a complicated energy function with lots of terms.
CustomCompoundBondForce* compound = new CustomCompoundBondForce(4, "x1+y2+z4+fn(distance(p1,p2))*angle(p3,p2,p4)+scale*dihedral(p2,p1,p4,p3)");
CustomCentroidBondForce* centroid = new CustomCentroidBondForce(4, "x1+y2+z4+fn(distance(g1,g2))*angle(g3,g2,g4)+scale*dihedral(g2,g1,g4,g3)");
compound->addGlobalParameter("scale", 0.5);
centroid->addGlobalParameter("scale", 0.5);
compound->addTabulatedFunction("fn", new Continuous1DFunction(table, -1, 10));
centroid->addTabulatedFunction("fn", new Continuous1DFunction(table, -1, 10));
// Add two bonds to the CustomCompoundBondForce.
vector<int> particles(4);
vector<double> parameters;
particles[0] = 0;
particles[1] = 1;
particles[2] = 2;
particles[3] = 3;
compound->addBond(particles, parameters);
particles[0] = 2;
particles[1] = 4;
particles[2] = 3;
particles[3] = 1;
compound->addBond(particles, parameters);
// Add identical bonds to the CustomCentroidBondForce. As a stronger test, make sure that
// group number is different from particle number.
vector<int> groupMembers(1);
groupMembers[0] = 3;
centroid->addGroup(groupMembers);
groupMembers[0] = 0;
centroid->addGroup(groupMembers);
groupMembers[0] = 1;
centroid->addGroup(groupMembers);
groupMembers[0] = 2;
centroid->addGroup(groupMembers);
groupMembers[0] = 4;
centroid->addGroup(groupMembers);
vector<int> groups(4);
groups[0] = 1;
groups[1] = 2;
groups[2] = 3;
groups[3] = 0;
centroid->addBond(groups, parameters);
groups[0] = 3;
groups[1] = 4;
groups[2] = 0;
groups[3] = 2;
centroid->addBond(groups, parameters);
// Add both forces as different force groups, and create a context.
centroid->setForceGroup(1);
system.addForce(compound);
system.addForce(centroid);
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
// Evaluate the force and energy for various positions and see if they match.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < 10; i++) {
for (int j = 0; j < numParticles; j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
context.setPositions(positions);
State state1 = context.getState(State::Forces | State::Energy, false, 1<<0);
State state2 = context.getState(State::Forces | State::Energy, false, 1<<1);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), TOL);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], TOL);
}
}
void testCustomWeights() {
System system;
system.addParticle(1.0);
system.addParticle(2.0);
system.addParticle(3.0);
system.addParticle(4.0);
CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "distance(g1,g2)^2");
vector<int> particles(2);
vector<double> weights(2);
particles[0] = 0;
particles[1] = 1;
weights[0] = 0.5;
weights[1] = 1.5;
force->addGroup(particles, weights);
particles[0] = 2;
particles[1] = 3;
weights[0] = 2.0;
weights[1] = 1.0;
force->addGroup(particles, weights);
vector<int> groups;
groups.push_back(0);
groups.push_back(1);
vector<double> parameters;
force->addBond(groups, parameters);
system.addForce(force);
// The center of mass of group 0 is (0, 1, 0).
vector<Vec3> positions(4);
positions[0] = Vec3(0, 4, 0);
positions[1] = Vec3(0, 0, 0);
// The center of mass of group 1 is (0, 10, 0).
positions[2] = Vec3(0, 9, 0);
positions[3] = Vec3(0, 12, 0);
// Check the forces and energy.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(9.0*9.0, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(Vec3(0, 2*9*(0.5/2.0), 0), state.getForces()[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 2*9*(1.5/2.0), 0), state.getForces()[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, -2*9*(2.0/3.0), 0), state.getForces()[2], TOL);
ASSERT_EQUAL_VEC(Vec3(0, -2*9*(1.0/3.0), 0), state.getForces()[3], TOL);
}
int main() {
try {
testHarmonicBond();
testComplexFunction();
testCustomWeights();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/reference/tests/TestReferenceCustomCompoundBondForce.cpp
View file @ ccd811da
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2012-2015 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,306 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of CustomCompoundBondForce.
*/
#include "ReferenceTests.h"
#include "TestCustomCompoundBondForce.h"
#ifdef WIN32
#define _USE_MATH_DEFINES // Needed to get M_PI
#endif
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomCompoundBondForce.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testBond() {
// Create a system using a CustomCompoundBondForce.
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(4, "0.5*kb*((distance(p1,p2)-b0)^2+(distance(p2,p3)-b0)^2)+0.5*ka*(angle(p2,p3,p4)-a0)^2+kt*(1+cos(dihedral(p1,p2,p3,p4)-t0))");
custom->addPerBondParameter("kb");
custom->addPerBondParameter("ka");
custom->addPerBondParameter("kt");
custom->addPerBondParameter("b0");
custom->addPerBondParameter("a0");
custom->addPerBondParameter("t0");
vector<int> particles(4);
particles[0] = 0;
particles[1] = 1;
particles[2] = 3;
particles[3] = 2;
vector<double> parameters(6);
parameters[0] = 1.5;
parameters[1] = 0.8;
parameters[2] = 0.6;
parameters[3] = 1.1;
parameters[4] = 2.9;
parameters[5] = 1.3;
custom->addBond(particles, parameters);
customSystem.addForce(custom);
ASSERT(!custom->usesPeriodicBoundaryConditions());
ASSERT(!customSystem.usesPeriodicBoundaryConditions());
// Create an identical system using standard forces.
System standardSystem;
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
HarmonicBondForce* bonds = new HarmonicBondForce();
bonds->addBond(0, 1, 1.1, 1.5);
bonds->addBond(1, 3, 1.1, 1.5);
standardSystem.addForce(bonds);
HarmonicAngleForce* angles = new HarmonicAngleForce();
angles->addAngle(1, 3, 2, 2.9, 0.8);
standardSystem.addForce(angles);
PeriodicTorsionForce* torsions = new PeriodicTorsionForce();
torsions->addTorsion(0, 1, 3, 2, 1, 1.3, 0.6);
standardSystem.addForce(torsions);
// Set the atoms in various positions, and verify that both systems give identical forces and energy.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context c1(customSystem, integrator1, platform);
Context c2(standardSystem, integrator2, platform);
vector<Vec3> positions(4);
for (int i = 0; i < 10; i++) {
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
// Try changing the bond parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[3] = 1.3;
custom->setBondParameters(0, particles, parameters);
custom->updateParametersInContext(c1);
bonds->setBondParameters(0, 0, 1, 1.3, 1.6);
bonds->setBondParameters(1, 1, 3, 1.3, 1.6);
bonds->updateParametersInContext(c2);
{
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = s1.getForces();
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
}
void testPositionDependence() {
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(2, "scale1*distance(p1,p2)+scale2*x1+2*y2");
custom->addGlobalParameter("scale1", 0.3);
custom->addGlobalParameter("scale2", 0.2);
vector<int> particles(2);
particles[0] = 1;
particles[1] = 0;
vector<double> parameters;
custom->addBond(particles, parameters);
customSystem.addForce(custom);
vector<Vec3> positions(2);
positions[0] = Vec3(1.5, 1, 0);
positions[1] = Vec3(0.5, 1, 0);
VerletIntegrator integrator(0.01);
Context context(customSystem, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(0.3*1.0+0.2*0.5+2*1, state.getPotentialEnergy(), 1e-5);
ASSERT_EQUAL_VEC(Vec3(-0.3, -2, 0), state.getForces()[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3-0.2, 0, 0), state.getForces()[1], 1e-5);
}
void testContinuous2DFunction() {
const int xsize = 10;
const int ysize = 11;
const double xmin = 0.4;
const double xmax = 1.1;
const double ymin = 0.0;
const double ymax = 0.9;
System system;
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomCompoundBondForce* forceField = new CustomCompoundBondForce(1, "fn(x1,y1)+1");
vector<int> particles(1, 0);
forceField->addBond(particles, vector<double>());
vector<double> table(xsize*ysize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
double x = xmin + i*(xmax-xmin)/xsize;
double y = ymin + j*(ymax-ymin)/ysize;
table[i+xsize*j] = sin(0.25*x)*cos(0.33*y);
}
}
forceField->addTabulatedFunction("fn", new Continuous2DFunction(xsize, ysize, table, xmin, xmax, ymin, ymax));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
positions[0] = Vec3(x, y, 1.5);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double energy = 1;
Vec3 force(0, 0, 0);
if (x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
energy = sin(0.25*x)*cos(0.33*y)+1;
force[0] = -0.25*cos(0.25*x)*cos(0.33*y);
force[1] = 0.3*sin(0.25*x)*sin(0.33*y);
}
ASSERT_EQUAL_VEC(force, forces[0], 0.1);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.05);
}
}
}
void testContinuous3DFunction() {
const int xsize = 10;
const int ysize = 11;
const int zsize = 12;
const double xmin = 0.4;
const double xmax = 1.1;
const double ymin = 0.0;
const double ymax = 0.9;
const double zmin = 0.2;
const double zmax = 1.3;
System system;
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomCompoundBondForce* forceField = new CustomCompoundBondForce(1, "fn(x1,y1,z1)+1");
vector<int> particles(1, 0);
forceField->addBond(particles, vector<double>());
vector<double> table(xsize*ysize*zsize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
for (int k = 0; k < zsize; k++) {
double x = xmin + i*(xmax-xmin)/xsize;
double y = ymin + j*(ymax-ymin)/ysize;
double z = zmin + k*(zmax-zmin)/zsize;
table[i+xsize*j+xsize*ysize*k] = sin(0.25*x)*cos(0.33*y)*(1+z);
}
}
}
forceField->addTabulatedFunction("fn", new Continuous3DFunction(xsize, ysize, zsize, table, xmin, xmax, ymin, ymax, zmin, zmax));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
for (double x = xmin-0.15; x < xmax+0.2; x += 0.1) {
for (double y = ymin-0.15; y < ymax+0.2; y += 0.1) {
for (double z = zmin-0.15; z < zmax+0.2; z += 0.1) {
positions[0] = Vec3(x, y, z);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double energy = 1;
Vec3 force(0, 0, 0);
if (x >= xmin && x <= xmax && y >= ymin && y <= ymax && z >= zmin && z <= zmax) {
energy = sin(0.25*x)*cos(0.33*y)*(1.0+z)+1;
force[0] = -0.25*cos(0.25*x)*cos(0.33*y)*(1.0+z);
force[1] = 0.3*sin(0.25*x)*sin(0.33*y)*(1.0+z);
force[2] = -sin(0.25*x)*cos(0.33*y);
}
ASSERT_EQUAL_VEC(force, forces[0], 0.1);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.05);
}
}
}
}
void testMultipleBonds() {
// Two compound bonds using Urey-Bradley example from API doc
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(3,
"0.5*(kangle*(angle(p1,p2,p3)-theta0)^2+kbond*(distance(p1,p3)-r0)^2)");
custom->addPerBondParameter("kangle");
custom->addPerBondParameter("kbond");
custom->addPerBondParameter("theta0");
custom->addPerBondParameter("r0");
vector<double> parameters(4);
parameters[0] = 1.0;
parameters[1] = 1.0;
parameters[2] = 2 * M_PI / 3;
parameters[3] = sqrt(3.0) / 2;
vector<int> particles0(3);
particles0[0] = 0;
particles0[1] = 1;
particles0[2] = 2;
vector<int> particles1(3);
particles1[0] = 1;
particles1[1] = 2;
particles1[2] = 3;
custom->addBond(particles0, parameters);
custom->addBond(particles1, parameters);
customSystem.addForce(custom);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0.5, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0.5, 0, 0);
positions[3] = Vec3(0.6, 0, 0.4);
VerletIntegrator integrator(0.01);
Context context(customSystem, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(0.199, state.getPotentialEnergy(), 1e-3);
vector<Vec3> forces(state.getForces());
ASSERT_EQUAL_VEC(Vec3(-1.160, 0.112, 0.0), forces[0], 1e-3);
ASSERT_EQUAL_VEC(Vec3(0.927, 1.047, -0.638), forces[1], 1e-3);
ASSERT_EQUAL_VEC(Vec3(-0.543, -1.160, 0.721), forces[2], 1e-3);
ASSERT_EQUAL_VEC(Vec3(0.776, 0.0, -0.084), forces[3], 1e-3);
}
int main() {
try {
testBond();
testPositionDependence();
testContinuous2DFunction();
testContinuous3DFunction();
testMultipleBonds();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/reference/tests/TestReferenceCustomExternalForce.cpp
View file @ ccd811da
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2015 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,132 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of CustomExternalForce.
*/
#include "ReferenceTests.h"
#include "TestCustomExternalForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testForce() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomExternalForce* forceField = new CustomExternalForce("scale*(x+yscale*(y-y0)^2)");
forceField->addPerParticleParameter("y0");
forceField->addPerParticleParameter("yscale");
forceField->addGlobalParameter("scale", 0.5);
vector<double> parameters(2);
parameters[0] = 0.5;
parameters[1] = 2.0;
forceField->addParticle(0, parameters);
parameters[0] = 1.5;
parameters[1] = 3.0;
forceField->addParticle(2, parameters);
system.addForce(forceField);
ASSERT(!forceField->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 2, 0);
positions[1] = Vec3(0, 0, 1);
positions[2] = Vec3(1, 0, 1);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-0.5, -0.5*2.0*2.0*1.5, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(-0.5, 0.5*3.0*2.0*1.5, 0), forces[2], TOL);
ASSERT_EQUAL_TOL(0.5*(1.0 + 2.0*1.5*1.5 + 3.0*1.5*1.5), state.getPotentialEnergy(), TOL);
}
// Try changing the parameters and make sure it's still correct.
parameters[0] = 1.4;
parameters[1] = 3.5;
forceField->setParticleParameters(1, 2, parameters);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-0.5, -0.5*2.0*2.0*1.5, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(-0.5, 0.5*3.5*2.0*1.4, 0), forces[2], TOL);
ASSERT_EQUAL_TOL(0.5*(1.0 + 2.0*1.5*1.5 + 3.5*1.4*1.4), state.getPotentialEnergy(), TOL);
}
}
void testPeriodic() {
Vec3 vx(5, 0, 0);
Vec3 vy(0, 6, 0);
Vec3 vz(1, 2, 7);
double x0 = 51, y0 = -17, z0 = 11.2;
System system;
system.setDefaultPeriodicBoxVectors(vx, vy, vz);
system.addParticle(1.0);
CustomExternalForce* force = new CustomExternalForce("periodicdistance(x, y, z, x0, y0, z0)^2");
force->addPerParticleParameter("x0");
force->addPerParticleParameter("y0");
force->addPerParticleParameter("z0");
vector<double> params(3);
params[0] = x0;
params[1] = y0;
params[2] = z0;
force->addParticle(0, params);
system.addForce(force);
ASSERT(force->usesPeriodicBoundaryConditions());
ASSERT(system.usesPeriodicBoundaryConditions());
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, 2, 0);
context.setPositions(positions);
for (int i = 0; i < 100; i++) {
State state = context.getState(State::Positions | State::Forces | State::Energy);
// Apply periodic boundary conditions to the difference between the two positions.
Vec3 delta = Vec3(x0, y0, z0)-state.getPositions()[0];
delta -= vz*floor(delta[2]/vz[2]+0.5);
delta -= vy*floor(delta[1]/vy[1]+0.5);
delta -= vx*floor(delta[0]/vx[0]+0.5);
// Verify that the force and energy are correct.
ASSERT_EQUAL_VEC(delta*2, state.getForces()[0], 1e-6);
ASSERT_EQUAL_TOL(delta.dot(delta), state.getPotentialEnergy(), 1e-6);
integrator.step(1);
}
void runPlatformTests() {
}
int main() {
try {
testForce();
testPeriodic();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/reference/tests/TestReferencePeriodicTorsionForce.cpp
View file @ ccd811da
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,80 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests all the different force terms in the reference implementation of PeriodicTorsionForce.
*/
#include "ReferenceTests.h"
#include "TestPeriodicTorsionForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testPeriodicTorsions() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
PeriodicTorsionForce* forceField = new PeriodicTorsionForce();
forceField->addTorsion(0, 1, 2, 3, 2, PI_M/3, 1.1);
system.addForce(forceField);
ASSERT(!forceField->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context 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);
}
// Try changing the torsion parameters and make sure it's still correct.
forceField->setTorsionParameters(0, 0, 1, 2, 3, 3, PI_M/3.2, 1.3);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double dtheta = (3*PI_M/2)-(PI_M/3.2);
double torque = -3*1.3*std::sin(dtheta);
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.3*(1+std::cos(dtheta)), state.getPotentialEnergy(), TOL);
}
}
int main() {
try {
testPeriodicTorsions();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/reference/tests/TestReferenceRBTorsionForce.cpp
View file @ ccd811da
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,99 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests all the different force terms in the reference implementation of RBTorsionForce.
*/
#include "ReferenceTests.h"
#include "TestRBTorsionForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/RBTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testRBTorsions() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
RBTorsionForce* forceField = new RBTorsionForce();
forceField->addTorsion(0, 1, 2, 3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6);
system.addForce(forceField);
ASSERT(!forceField->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context 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);
}
// Try changing the torsion parameters and make sure it's still correct.
forceField->setTorsionParameters(0, 0, 1, 2, 3, 0.11, 0.22, 0.33, 0.44, 0.55, 0.66);
forceField->updateParametersInContext(context);
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.11*(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.11*(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;
void runPlatformTests() {
}
tests/TestCustomAngleForce.h 0 → 100644
View file @ ccd811da
/* -------------------------------------------------------------------------- *
* 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-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/CustomAngleForce.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testAngles() {
// Create a system using a CustomAngleForce.
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomAngleForce* custom = new CustomAngleForce("scale*k*(theta-theta0)^2");
custom->addPerAngleParameter("theta0");
custom->addPerAngleParameter("k");
custom->addGlobalParameter("scale", 0.5);
vector<double> parameters(2);
parameters[0] = 1.5;
parameters[1] = 0.8;
custom->addAngle(0, 1, 2, parameters);
parameters[0] = 2.0;
parameters[1] = 0.5;
custom->addAngle(1, 2, 3, parameters);
customSystem.addForce(custom);
ASSERT(!custom->usesPeriodicBoundaryConditions());
ASSERT(!customSystem.usesPeriodicBoundaryConditions());
// Create an identical system using a HarmonicAngleForce.
System harmonicSystem;
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
HarmonicAngleForce* harmonic = new HarmonicAngleForce();
harmonic->addAngle(0, 1, 2, 1.5, 0.8);
harmonic->addAngle(1, 2, 3, 2.0, 0.5);
harmonicSystem.addForce(harmonic);
// Set the atoms in various positions, and verify that both systems give identical forces and energy.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(4);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context c1(customSystem, integrator1, platform);
Context c2(harmonicSystem, integrator2, platform);
for (int i = 0; i < 10; i++) {
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
// Try changing the angle parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[1] = 0.9;
custom->setAngleParameters(0, 0, 1, 2, parameters);
parameters[0] = 2.1;
parameters[1] = 0.6;
custom->setAngleParameters(1, 1, 2, 3, parameters);
custom->updateParametersInContext(c1);
harmonic->setAngleParameters(0, 0, 1, 2, 1.6, 0.9);
harmonic->setAngleParameters(1, 1, 2, 3, 2.1, 0.6);
harmonic->updateParametersInContext(c2);
{
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
}
void runPlatformTests();
int main(int argc, char* argv[]) {
try {
initializeTests(argc, argv);
testAngles();
runPlatformTests();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
tests/TestCustomBondForce.h 0 → 100644
View file @ ccd811da
/* -------------------------------------------------------------------------- *
* 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-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/CustomBondForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testBonds() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomBondForce* forceField = new CustomBondForce("scale*k*(r-r0)^2");
forceField->addPerBondParameter("r0");
forceField->addPerBondParameter("k");
forceField->addGlobalParameter("scale", 0.5);
vector<double> parameters(2);
parameters[0] = 1.5;
parameters[1] = 0.8;
forceField->addBond(0, 1, parameters);
parameters[0] = 1.2;
parameters[1] = 0.7;
forceField->addBond(1, 2, parameters);
system.addForce(forceField);
ASSERT(!forceField->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context 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);
}
// Try changing the bond parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[1] = 0.9;
forceField->setBondParameters(0, 0, 1, parameters);
parameters[0] = 1.3;
parameters[1] = 0.8;
forceField->setBondParameters(1, 1, 2, parameters);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, -0.9*0.4, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.8*0.3, 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.9*0.4*0.4 + 0.5*0.8*0.3*0.3, state.getPotentialEnergy(), TOL);
}
}
void testManyParameters() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomBondForce* forceField = new CustomBondForce("(a+b+c+d+e+f+g+h+i)*r");
forceField->addPerBondParameter("a");
forceField->addPerBondParameter("b");
forceField->addPerBondParameter("c");
forceField->addPerBondParameter("d");
forceField->addPerBondParameter("e");
forceField->addPerBondParameter("f");
forceField->addPerBondParameter("g");
forceField->addPerBondParameter("h");
forceField->addPerBondParameter("i");
vector<double> parameters(forceField->getNumPerBondParameters());
for (int i = 0; i < parameters.size(); i++)
parameters[i] = i;
forceField->addBond(0, 1, parameters);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(0, 2.5, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double f = 1+2+3+4+5+6+7+8;
ASSERT_EQUAL_VEC(Vec3(0, f, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, -f, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(f*2.5, state.getPotentialEnergy(), TOL);
}
void runPlatformTests();
int main(int argc, char* argv[]) {
try {
initializeTests(argc, argv);
testBonds();
testManyParameters();
runPlatformTests();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
tests/TestCustomCentroidBondForce.h 0 → 100644
View file @ ccd811da
This diff is collapsed.
tests/TestCustomCompoundBondForce.h 0 → 100644
View file @ ccd811da
This diff is collapsed.
tests/TestCustomExternalForce.h 0 → 100644
View file @ ccd811da
This diff is collapsed.
tests/TestCustomGBForce.h 0 → 100644
View file @ ccd811da
This diff is collapsed.
tests/TestHarmonicAngleForce.h
View file @ ccd811da
......@@ -29,10 +29,6 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests all the different force terms in the reference implementation of HarmonicAngleForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/HarmonicAngleForce.h"
......
tests/TestPeriodicTorsionForce.h 0 → 100644
View file @ ccd811da
/* -------------------------------------------------------------------------- *
* 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-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testPeriodicTorsions() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
PeriodicTorsionForce* forceField = new PeriodicTorsionForce();
forceField->addTorsion(0, 1, 2, 3, 2, PI_M/3, 1.1);
system.addForce(forceField);
ASSERT(!forceField->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context 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);
}
// Try changing the torsion parameters and make sure it's still correct.
forceField->setTorsionParameters(0, 0, 1, 2, 3, 3, PI_M/3.2, 1.3);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double dtheta = (3*PI_M/2)-(PI_M/3.2);
double torque = -3*1.3*std::sin(dtheta);
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.3*(1+std::cos(dtheta)), state.getPotentialEnergy(), TOL);
}
}
void runPlatformTests();
int main(int argc, char* argv[]) {
try {
initializeTests(argc, argv);
testPeriodicTorsions();
runPlatformTests();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
tests/TestRBTorsionForce.h 0 → 100644
View file @ ccd811da
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