Bug fixes and additional test cases

platforms/cuda2/src/CudaSort.h

@@ -46,9 +46,9 @@ namespace OpenMM {
  *     int getKeySize() const {return 4;}
  *     const char* getDataType() const {return "float";}
  *     const char* getKeyType() const {return "float";}
- *     const char* getMinKey() const {return "-MAXFLOAT";}
- *     const char* getMaxKey() const {return "MAXFLOAT";}
- *     const char* getMaxValue() const {return "MAXFLOAT";}
+ *     const char* getMinKey() const {return "-3.40282e+38f";}
+ *     const char* getMaxKey() const {return "3.40282e+38f";}
+ *     const char* getMaxValue() const {return "3.40282e+38f";}
  *     const char* getSortKey() const {return "value";}
  * };
  *

platforms/cuda2/src/kernels/findInteractingBlocks.cu

@@ -218,7 +218,7 @@ extern "C" __global__ void findInteractionsWithinBlocks(real4 periodicBoxSize, r
             if (index == 0)
                 allFlags = flags[threadIdx.x]+flags[threadIdx.x+4]+flags[threadIdx.x+8]+flags[threadIdx.x+12]+flags[threadIdx.x+16]+flags[threadIdx.x+20]+flags[threadIdx.x+24]+flags[threadIdx.x+28];
 #else
-            unsigned int allFlags = __ballot(delta.x*delta.x+delta.y*delta.y+delta.z*delta.z > CUTOFF_SQUARED);
+            unsigned int allFlags = __ballot(delta.x*delta.x+delta.y*delta.y+delta.z*delta.z < CUTOFF_SQUARED);
 #endif
 
             // Sum the flags.

platforms/cuda2/tests/TestCudaLocalEnergyMinimizer.cpp

+
+/* -------------------------------------------------------------------------- *
+ *                                   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) 2010-2012 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 "CudaPlatform.h"
+#include "openmm/Context.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/LocalEnergyMinimizer.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/VerletIntegrator.h"
+#include "openmm/VirtualSite.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testHarmonicBonds() {
+    const int numParticles = 10;
+    System system;
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+
+    // Create a chain of particles connected by harmonic bonds.
+
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        positions[i] = Vec3(i, 0, 0);
+        if (i > 0)
+            bonds->addBond(i-1, i, 1+0.1*i, 1);
+    }
+
+    // Minimize it and check that all bonds are at their equilibrium distances.
+
+    VerletIntegrator integrator(0.01);
+    CudaPlatform platform;
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    LocalEnergyMinimizer::minimize(context, 1e-5);
+    State state = context.getState(State::Positions);
+    for (int i = 1; i < numParticles; i++) {
+        Vec3 delta = state.getPositions()[i]-state.getPositions()[i-1];
+        ASSERT_EQUAL_TOL(1+0.1*i, sqrt(delta.dot(delta)), 1e-4);
+    }
+}
+
+void testLargeSystem() {
+    const int numMolecules = 50;
+    const int numParticles = numMolecules*2;
+    const double cutoff = 2.0;
+    const double boxSize = 5.0;
+    const double tolerance = 5;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    system.addForce(nonbonded);
+
+    // Create a cloud of molecules.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        nonbonded->addParticle(-1.0, 0.2, 0.2);
+        nonbonded->addParticle(1.0, 0.2, 0.2);
+        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        system.addConstraint(2*i, 2*i+1, 1.0);
+    }
+
+    // Minimize it and verify that the energy has decreased.
+
+    CudaPlatform platform;
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State initialState = context.getState(State::Forces | State::Energy);
+    LocalEnergyMinimizer::minimize(context, tolerance);
+    State finalState = context.getState(State::Forces | State::Energy | State::Positions);
+    ASSERT(finalState.getPotentialEnergy() < initialState.getPotentialEnergy());
+
+    // Compute the force magnitude, substracting off any component parallel to a constraint, and
+    // check that it satisfies the requested tolerance.
+
+    double forceNorm = 0.0;
+    for (int i = 0; i < numParticles; i += 2) {
+        Vec3 dir = finalState.getPositions()[i+1]-finalState.getPositions()[i];
+        double distance = sqrt(dir.dot(dir));
+        dir *= 1.0/distance;
+        Vec3 f = finalState.getForces()[i];
+        f -= dir*dir.dot(f);
+        forceNorm += f.dot(f);
+        f = finalState.getForces()[i+1];
+        f -= dir*dir.dot(f);
+        forceNorm += f.dot(f);
+    }
+    forceNorm = sqrt(forceNorm/(4*numMolecules));
+    ASSERT(forceNorm < 3*tolerance);
+}
+
+void testVirtualSites() {
+    const int numMolecules = 50;
+    const int numParticles = numMolecules*3;
+    const double cutoff = 2.0;
+    const double boxSize = 5.0;
+    const double tolerance = 5;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    system.addForce(nonbonded);
+
+    // Create a cloud of molecules.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        system.addParticle(0.0);
+        nonbonded->addParticle(-1.0, 0.2, 0.2);
+        nonbonded->addParticle(0.5, 0.2, 0.2);
+        nonbonded->addParticle(0.5, 0.2, 0.2);
+        positions[3*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[3*i+1] = Vec3(positions[3*i][0]+1.0, positions[3*i][1], positions[3*i][2]);
+        positions[3*i+2] = Vec3();
+        system.addConstraint(3*i, 3*i+1, 1.0);
+        system.setVirtualSite(3*i+2, new TwoParticleAverageSite(3*i, 3*i+1, 0.5, 0.5));
+    }
+
+    // Minimize it and verify that the energy has decreased.
+    
+    CudaPlatform platform;
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State initialState = context.getState(State::Forces | State::Energy);
+    LocalEnergyMinimizer::minimize(context, tolerance);
+    State finalState = context.getState(State::Forces | State::Energy | State::Positions);
+    ASSERT(finalState.getPotentialEnergy() < initialState.getPotentialEnergy());
+
+    // Compute the force magnitude, subtracting off any component parallel to a constraint, and
+    // check that it satisfies the requested tolerance.
+
+    double forceNorm = 0.0;
+    for (int i = 0; i < numParticles; i += 3) {
+        Vec3 dir = finalState.getPositions()[i+1]-finalState.getPositions()[i];
+        double distance = sqrt(dir.dot(dir));
+        dir *= 1.0/distance;
+        Vec3 f = finalState.getForces()[i];
+        f -= dir*dir.dot(f);
+        forceNorm += f.dot(f);
+        f = finalState.getForces()[i+1];
+        f -= dir*dir.dot(f);
+        forceNorm += f.dot(f);
+        
+        // Check the virtual site location.
+        
+        ASSERT_EQUAL_VEC((finalState.getPositions()[i+1]+finalState.getPositions()[i])*0.5, finalState.getPositions()[i+2], 1e-5);
+    }
+    forceNorm = sqrt(forceNorm/(4*numMolecules));
+    ASSERT(forceNorm < 3*tolerance);
+}
+
+int main() {
+    try {
+        testHarmonicBonds();
+        testLargeSystem();
+        testVirtualSites();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cuda2/tests/TestCudaMultipleForces.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   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) 2011-2012 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 a system with multiple forces, to make sure CudaBondedUtilities is
+ * processing them correctly.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "ReferencePlatform.h"
+#include "openmm/HarmonicAngleForce.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/PeriodicTorsionForce.h"
+#include "openmm/RBTorsionForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-4;
+
+void testForces() {
+    const int numParticles = 100;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    for (int i = 0; i < numParticles-1; i++)
+        bonds->addBond(i, i+1, 1.0, 1.5);
+    system.addForce(bonds);
+    HarmonicAngleForce* angles = new HarmonicAngleForce();
+    for (int i = 0; i < numParticles-2; i++)
+        angles->addAngle(i, i+1, i+2, 2.0, 1.5);
+    system.addForce(angles);
+    PeriodicTorsionForce* periodic = new PeriodicTorsionForce();
+    for (int i = 0; i < numParticles-3; i++)
+        periodic->addTorsion(i, i+1, i+2, i+3, 2, PI_M/3, 1.1);
+    system.addForce(periodic);
+    RBTorsionForce* rb = new RBTorsionForce();
+    for (int i = 0; i < numParticles-3; i += 3)
+        rb->addTorsion(i, i+1, i+2, i+3, 1.0, 1.1, 1.2, 0.3, 0.4, 0.5);
+    system.addForce(rb);
+    ReferencePlatform ref;
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, ref);
+    CudaPlatform cl;
+    VerletIntegrator integrator2(0.01);
+    Context context2(system, integrator2, cl);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        positions[i] = Vec3(i, genrand_real2(sfmt), genrand_real2(sfmt));
+    context1.setPositions(positions);
+    context2.setPositions(positions);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces1 = state1.getForces();
+    const vector<Vec3>& forces2 = state2.getForces();
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(forces1[i], forces2[i], TOL);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), TOL);
+}
+
+int main() {
+    try {
+        testForces();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cuda2/tests/TestCudaSettle.cpp

+
+/* -------------------------------------------------------------------------- *
+ *                                   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-2012 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 the SETTLE algorithm.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testConstraints() {
+    const int numMolecules = 10;
+    const int numParticles = numMolecules*3;
+    const int numConstraints = numMolecules*3;
+    const double temp = 100.0;
+    CudaPlatform platform;
+    System system;
+    LangevinIntegrator integrator(temp, 2.0, 0.001);
+    integrator.setConstraintTolerance(1e-5);
+    NonbondedForce* forceField = new NonbondedForce();
+    for (int i = 0; i < numMolecules; ++i) {
+        system.addParticle(16.0);
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        forceField->addParticle(-0.82, 0.317, 0.65);
+        forceField->addParticle(0.41, 1.0, 0.0);
+        forceField->addParticle(0.41, 1.0, 0.0);
+        system.addConstraint(i*3, i*3+1, 0.1);
+        system.addConstraint(i*3, i*3+2, 0.1);
+        system.addConstraint(i*3+1, i*3+2, 0.163);
+    }
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numMolecules; ++i) {
+        positions[i*3] = Vec3((i%4)*0.4, (i/4)*0.4, 0);
+        positions[i*3+1] = positions[i*3]+Vec3(0.1, 0, 0);
+        positions[i*3+2] = positions[i*3]+Vec3(-0.03333, 0.09428, 0);
+        velocities[i*3] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
+        velocities[i*3+1] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
+        velocities[i*3+2] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
+    }
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+
+    // Simulate it and see whether the constraints remain satisfied.
+
+    for (int i = 0; i < 1000; ++i) {
+        integrator.step(1);
+        State state = context.getState(State::Positions | State::Forces);
+        for (int j = 0; j < numConstraints; ++j) {
+            int particle1, particle2;
+            double distance;
+            system.getConstraintParameters(j, particle1, particle2, distance);
+            Vec3 p1 = state.getPositions()[particle1];
+            Vec3 p2 = state.getPositions()[particle2];
+            double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
+            ASSERT_EQUAL_TOL(distance, dist, 1e-5);
+        }
+    }
+}
+
+int main() {
+    try {
+        testConstraints();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cuda2/tests/TestCudaSort.cpp

@@ -51,9 +51,9 @@ class SortTrait : public CudaSort::SortTrait {
     int getKeySize() const {return 4;}
     const char* getDataType() const {return "float";}
     const char* getKeyType() const {return "float";}
-    const char* getMinKey() const {return "-MAXFLOAT";}
-    const char* getMaxKey() const {return "MAXFLOAT";}
-    const char* getMaxValue() const {return "MAXFLOAT";}
+    const char* getMinKey() const {return "-3.40282e+38f";}
+    const char* getMaxKey() const {return "3.40282e+38f";}
+    const char* getMaxValue() const {return "3.40282e+38f";}
     const char* getSortKey() const {return "value";}
 };