Merge branch 'master' of github.com:pandegroup/openmm into dpme

platforms/reference/src/SimTKReference/ReferenceCCMAAlgorithm.cpp

 
-/* Portions copyright (c) 2006-2015 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2017 Stanford University and Simbios.
  * Contributors: Peter Eastman, Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -38,42 +38,6 @@
 using namespace OpenMM;
 using namespace std;
 
-// This class extracts columns from the inverse matrix one at a time.  It is done in parallel,
-// since this can be very slow.
-
-class ExtractMatrixTask : public ThreadPool::Task {
-public:
-    ExtractMatrixTask(int numConstraints, vector<vector<pair<int, RealOpenMM> > >& transposedMatrix, const vector<RealOpenMM>& distance, RealOpenMM elementCutoff,
-                      const int* qRowStart, const int* qColIndex, const int* rRowStart, const int* rColIndex, const double* qValue, const double* rValue) :
-                numConstraints(numConstraints), transposedMatrix(transposedMatrix), distance(distance), elementCutoff(elementCutoff), qRowStart(qRowStart), qColIndex(qColIndex),
-                rRowStart(rRowStart), rColIndex(rColIndex), qValue(qValue), rValue(rValue) {
-    }
-
-    void execute(ThreadPool& pool, int threadIndex) {
-        vector<double> rhs(numConstraints);
-        for (int i = threadIndex; i < numConstraints; i += pool.getNumThreads()) {
-            // Extract column i of the inverse matrix.
-
-            for (int j = 0; j < numConstraints; j++)
-                rhs[j] = (i == j ? 1.0 : 0.0);
-            QUERN_multiply_with_q_transpose(numConstraints, qRowStart, qColIndex, qValue, &rhs[0]);
-            QUERN_solve_with_r(numConstraints, rRowStart, rColIndex, rValue, &rhs[0], &rhs[0]);
-            for (int j = 0; j < numConstraints; j++) {
-                double value = rhs[j]*distance[i]/distance[j];
-                if (FABS((RealOpenMM) value) > elementCutoff)
-                    transposedMatrix[i].push_back(pair<int, RealOpenMM>(j, (RealOpenMM) value));
-            }
-        }
-    }
-private:
-    int numConstraints;
-    vector<vector<pair<int, RealOpenMM> > >& transposedMatrix;
-    const vector<RealOpenMM>& distance;
-    RealOpenMM elementCutoff;
-    const int *qRowStart, *qColIndex, *rRowStart, *rColIndex;
-    const double *qValue, *rValue;
-};
-
 ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm(int numberOfAtoms,
                                                int numberOfConstraints,
                                                const vector<pair<int, int> >& atomIndices,
@@ -194,9 +158,27 @@ ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm(int numberOfAtoms,
                 &qRowStart, &qColIndex, &qValue, &rRowStart, &rColIndex, &rValue);
         vector<vector<pair<int, RealOpenMM> > > transposedMatrix(numberOfConstraints);
         _matrix.resize(numberOfConstraints);
+
+        // Extract columns from the inverse matrix one at a time.  It is done in parallel,
+        // since this can be very slow.
+        
         ThreadPool threads;
-        ExtractMatrixTask task(numberOfConstraints, transposedMatrix, _distance, _elementCutoff, qRowStart, qColIndex, rRowStart, rColIndex, qValue, rValue);
-        threads.execute(task);
+        threads.execute([&] (ThreadPool& pool, int threadIndex) {
+            vector<double> rhs(numberOfConstraints);
+            for (int i = threadIndex; i < numberOfConstraints; i += pool.getNumThreads()) {
+                // Extract column i of the inverse matrix.
+
+                for (int j = 0; j < numberOfConstraints; j++)
+                    rhs[j] = (i == j ? 1.0 : 0.0);
+                QUERN_multiply_with_q_transpose(numberOfConstraints, qRowStart, qColIndex, qValue, &rhs[0]);
+                QUERN_solve_with_r(numberOfConstraints, rRowStart, rColIndex, rValue, &rhs[0], &rhs[0]);
+                for (int j = 0; j < numberOfConstraints; j++) {
+                    double value = rhs[j]*distance[i]/distance[j];
+                    if (FABS((RealOpenMM) value) > elementCutoff)
+                        transposedMatrix[i].push_back(pair<int, RealOpenMM>(j, (RealOpenMM) value));
+                }
+            }
+        });
         threads.waitForThreads();
 
         // For purposes of thread safety we extracted the matrix in transposed form, so we need to transpose it again.

plugins/amoeba/platforms/cuda/src/kernels/multipoleElectrostatics.cu

@@ -336,7 +336,7 @@ __device__ void computeOneInteraction(AtomData& atom1, AtomData& atom2, bool has
     iEIY -= eCoef*(qiUinpI.y*qiUindJ.x + qiUindI.y*qiUinpJ.x);
     iEJY -= eCoef*(qiUinpJ.y*qiUindI.x + qiUindJ.y*qiUinpI.x);
     fIZ += dCoef*(qiUinpI.x*qiUindJ.x + qiUindI.x*qiUinpJ.x);
-    fIZ += dCoef*(qiUinpJ.x*qiUindI.x + qiUindJ.x*qiUinpI.x);
+    fJZ += dCoef*(qiUinpJ.x*qiUindI.x + qiUindJ.x*qiUinpI.x);
     // Uind-Uind terms (m=1)
     eCoef = 2*rInvVec[3]*thole_d1;
     dCoef = -3*rInvVec[4]*dthole_d1;
@@ -345,7 +345,7 @@ __device__ void computeOneInteraction(AtomData& atom1, AtomData& atom2, bool has
     iEIY += eCoef*(qiUinpI.x*qiUindJ.y + qiUindI.x*qiUinpJ.y);
     iEJY += eCoef*(qiUinpJ.x*qiUindI.y + qiUindJ.x*qiUinpI.y);
     fIZ += dCoef*(qiUinpI.y*qiUindJ.y + qiUindI.y*qiUinpJ.y + qiUinpI.z*qiUindJ.z + qiUindI.z*qiUinpJ.z);
-    fIZ += dCoef*(qiUinpJ.y*qiUindI.y + qiUindJ.y*qiUinpI.y + qiUinpJ.z*qiUindI.z + qiUindJ.z*qiUinpI.z);
+    fJZ += dCoef*(qiUinpJ.y*qiUindI.y + qiUindJ.y*qiUinpI.y + qiUinpJ.z*qiUindI.z + qiUindJ.z*qiUinpI.z);
 #endif
 
     // The quasi-internal frame forces and torques.  Note that the induced torque intermediates are

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaExtrapolatedPolarization.cpp

@@ -284,8 +284,11 @@ vector<Vec3> setupWaterDimer(System& system,  AmoebaMultipoleForce* amoebaMultip
 
 static void check_finite_differences(vector<Vec3> analytic_forces, Context &context, vector<Vec3> positions)
 {
-    // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
+    double tol = 1e-5;
+    // We allow more permissive testing for single precision.
+    if(Platform::getPlatformByName("CUDA").getPropertyValue(context, "Precision") != "double") tol = 5e-4;
 
+    // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
     double norm = 0.0;
     for (int i = 0; i < (int) analytic_forces.size(); ++i)
         norm += analytic_forces[i].dot(analytic_forces[i]);
@@ -303,7 +306,7 @@ static void check_finite_differences(vector<Vec3> analytic_forces, Context &cont
     State state2 = context.getState(State::Energy);
     context.setPositions(positions3);
     State state3 = context.getState(State::Energy);
-    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
+    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, tol);
 }
 
 
@@ -493,6 +496,8 @@ int main(int numberOfArguments, char* argv[]) {
 
     try {
         registerAmoebaCudaKernelFactories();
+        if (numberOfArguments > 1)
+            Platform::getPlatformByName("CUDA").setPropertyDefaultValue("Precision", std::string(argv[1]));
 
         /*
          * Water dimer energy / force tests under various conditions.

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp

@@ -1509,7 +1509,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPairIxn(const Mu
         iEIY -= eCoef*(qiUinpI[1]*qiUindJ[0] + qiUindI[1]*qiUinpJ[0]);
         iEJY -= eCoef*(qiUinpJ[1]*qiUindI[0] + qiUindJ[1]*qiUinpI[0]);
         fIZ += dCoef*(qiUinpI[0]*qiUindJ[0] + qiUindI[0]*qiUinpJ[0]);
-        fIZ += dCoef*(qiUinpJ[0]*qiUindI[0] + qiUindJ[0]*qiUinpI[0]);
+        fJZ += dCoef*(qiUinpJ[0]*qiUindI[0] + qiUindJ[0]*qiUinpI[0]);
         // Uind-Uind terms (m=1)
         eCoef = 2.0*rInvVec[3]*uScale*thole_d1;
         dCoef = -3.0*rInvVec[4]*uScale*dthole_d1;
@@ -1518,7 +1518,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPairIxn(const Mu
         iEIY += eCoef*(qiUinpI[0]*qiUindJ[1] + qiUindI[0]*qiUinpJ[1]);
         iEJY += eCoef*(qiUinpJ[0]*qiUindI[1] + qiUindJ[0]*qiUinpI[1]);
         fIZ += dCoef*(qiUinpI[1]*qiUindJ[1] + qiUindI[1]*qiUinpJ[1] + qiUinpI[2]*qiUindJ[2] + qiUindI[2]*qiUinpJ[2]);
-        fIZ += dCoef*(qiUinpJ[1]*qiUindI[1] + qiUindJ[1]*qiUinpI[1] + qiUinpJ[2]*qiUindI[2] + qiUindJ[2]*qiUinpI[2]);
+        fJZ += dCoef*(qiUinpJ[1]*qiUindI[1] + qiUindJ[1]*qiUinpI[1] + qiUinpJ[2]*qiUindI[2] + qiUindJ[2]*qiUinpI[2]);
     }
 
     // The quasi-internal frame forces and torques.  Note that the induced torque intermediates are
@@ -6802,7 +6802,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeDirectElectrostaticPair
         iEIY -= eCoef*(qiUinpI[1]*qiUindJ[0] + qiUindI[1]*qiUinpJ[0]);
         iEJY -= eCoef*(qiUinpJ[1]*qiUindI[0] + qiUindJ[1]*qiUinpI[0]);
         fIZ += dCoef*(qiUinpI[0]*qiUindJ[0] + qiUindI[0]*qiUinpJ[0]);
-        fIZ += dCoef*(qiUinpJ[0]*qiUindI[0] + qiUindJ[0]*qiUinpI[0]);
+        fJZ += dCoef*(qiUinpJ[0]*qiUindI[0] + qiUindJ[0]*qiUinpI[0]);
         // Uind-Uind terms (m=1)
         eCoef = 2.0*rInvVec[3]*(uScale*thole_d1 + bVec[3] - twoThirds*alphaRVec[3]*X);
         dCoef = -3.0*rInvVec[4]*(uScale*dthole_d1 + bVec[3]);
@@ -6811,7 +6811,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeDirectElectrostaticPair
         iEIY += eCoef*(qiUinpI[0]*qiUindJ[1] + qiUindI[0]*qiUinpJ[1]);
         iEJY += eCoef*(qiUinpJ[0]*qiUindI[1] + qiUindJ[0]*qiUinpI[1]);
         fIZ += dCoef*(qiUinpI[1]*qiUindJ[1] + qiUindI[1]*qiUinpJ[1] + qiUinpI[2]*qiUindJ[2] + qiUindI[2]*qiUinpJ[2]);
-        fIZ += dCoef*(qiUinpJ[1]*qiUindI[1] + qiUindJ[1]*qiUinpI[1] + qiUinpJ[2]*qiUindI[2] + qiUindJ[2]*qiUinpI[2]);
+        fJZ += dCoef*(qiUinpJ[1]*qiUindI[1] + qiUindJ[1]*qiUinpI[1] + qiUinpJ[2]*qiUindI[2] + qiUindJ[2]*qiUinpI[2]);
     }
 
     // The quasi-internal frame forces and torques.  Note that the induced torque intermediates are

plugins/cpupme/src/CpuPmeKernels.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2017 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -400,16 +400,6 @@ static void interpolateForces(float* posq, float* force, float* grid, int gridx,
     }
 }
 
-class CpuCalcPmeReciprocalForceKernel::ComputeTask : public ThreadPool::Task {
-public:
-    ComputeTask(CpuCalcPmeReciprocalForceKernel& owner) : owner(owner) {
-    }
-    void execute(ThreadPool& threads, int threadIndex) {
-        owner.runWorkerThread(threads, threadIndex);
-    }
-    CpuCalcPmeReciprocalForceKernel& owner;
-};
-
 static void* threadBody(void* args) {
     CpuCalcPmeReciprocalForceKernel& owner = *reinterpret_cast<CpuCalcPmeReciprocalForceKernel*>(args);
     owner.runMainThread();
@@ -549,9 +539,8 @@ void CpuCalcPmeReciprocalForceKernel::runMainThread() {
         if (isDeleted)
             break;
         posq = io->getPosq();
-        ComputeTask task(*this);
         gmx_atomic_set(&atomicCounter, 0);
-        threads.execute(task); // Signal threads to perform charge spreading.
+        threads.execute([&] (ThreadPool& threads, int threadIndex) { runWorkerThread(threads, threadIndex); }); // Signal threads to perform charge spreading.
         threads.waitForThreads();
         threads.resumeThreads(); // Signal threads to sum the charge grids.
         threads.waitForThreads();

plugins/cpupme/src/CpuPmeKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2017 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -102,7 +102,6 @@ public:
      */
     void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
 private:
-    class ComputeTask;
     /**
      * Select a size for one grid dimension that FFTW can handle efficiently.
      */

wrappers/python/simtk/openmm/mtsintegrator.py

@@ -28,7 +28,7 @@ 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.
 """
-from __future__ import absolute_import
+from __future__ import absolute_import, division
 __author__ = "Peter Eastman"
 __version__ = "1.0"
 
@@ -90,9 +90,10 @@ class MTSIntegrator(CustomIntegrator):
 
     def _createSubsteps(self, parentSubsteps, groups):
         group, substeps = groups[0]
-        stepsPerParentStep = substeps//parentSubsteps
+        stepsPerParentStep = substeps / parentSubsteps
         if stepsPerParentStep < 1 or stepsPerParentStep != int(stepsPerParentStep):
             raise ValueError("The number for substeps for each group must be a multiple of the number for the previous group")
+        stepsPerParentStep = int(stepsPerParentStep)
         if group < 0 or group > 31:
             raise ValueError("Force group must be between 0 and 31")
         for i in range(stepsPerParentStep):

wrappers/python/tests/TestIntegrators.py

@@ -96,5 +96,29 @@ class TestIntegrators(unittest.TestCase):
             # Take a step
             integrator.step(1)
 
+    def testBadGroups(self):
+        """Test the MTS integrator with bad force group substeps."""
+        # Create a periodic solvated system with PME
+        pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')
+        ff = ForceField('amber99sbildn.xml', 'tip3p.xml')
+        system = ff.createSystem(pdb.topology, cutoffMethod=PME)
+
+        # Split forces into groups
+        for force in system.getForces():
+            if force.__class__.__name__ == 'NonbondedForce':
+                force.setForceGroup(1)
+                force.setReciprocalSpaceForceGroup(2)
+            else:
+                force.setForceGroup(0)
+
+        with self.assertRaises(ValueError):
+            # Create an integrator
+            integrator = MTSIntegrator(4*femtoseconds, [(2,1), (1,3), (0,8)])
+
+            # Run a few steps of dynamics
+            context = Context(system, integrator)
+            context.setPositions(pdb.positions)
+            integrator.step(10)
+
 if __name__ == '__main__':
     unittest.main()