Finished CUDA implementation of periodicdistance()

platforms/cuda/src/CudaBondedUtilities.cpp

@@ -99,7 +99,7 @@ void CudaBondedUtilities::initialize(const System& system) {
     s<<CudaKernelSources::vectorOps;
     for (int i = 0; i < (int) prefixCode.size(); i++)
         s<<prefixCode[i];
-    s<<"extern \"C\" __global__ void computeBondedForces(unsigned long long* __restrict__ forceBuffer, real* __restrict__ energyBuffer, const real4* __restrict__ posq, int groups";
+    s<<"extern \"C\" __global__ void computeBondedForces(unsigned long long* __restrict__ forceBuffer, real* __restrict__ energyBuffer, const real4* __restrict__ posq, int groups, real4 periodicBoxSize, real4 invPeriodicBoxSize, real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ";
     for (int force = 0; force < numForces; force++) {
         for (int i = 0; i < (int) atomIndices[force].size(); i++) {
             int indexWidth = atomIndices[force][i]->getElementSize()/4;
@@ -161,6 +161,11 @@ void CudaBondedUtilities::computeInteractions(int groups) {
         kernelArgs.push_back(&context.getEnergyBuffer().getDevicePointer());
         kernelArgs.push_back(&context.getPosq().getDevicePointer());
         kernelArgs.push_back(NULL);
+        kernelArgs.push_back(context.getPeriodicBoxSizePointer());
+        kernelArgs.push_back(context.getInvPeriodicBoxSizePointer());
+        kernelArgs.push_back(context.getPeriodicBoxVecXPointer());
+        kernelArgs.push_back(context.getPeriodicBoxVecYPointer());
+        kernelArgs.push_back(context.getPeriodicBoxVecZPointer());
         for (int i = 0; i < (int) atomIndices.size(); i++)
             for (int j = 0; j < (int) atomIndices[i].size(); j++)
                 kernelArgs.push_back(&atomIndices[i][j]->getDevicePointer());

platforms/cuda/src/CudaExpressionUtilities.cpp

@@ -101,15 +101,15 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express
             if (node.getOperation().getName() == "periodicdistance") {
                 // This is the periodicdistance() function.
 
-                out << tempType << " periodicDistance_delta = make_real3(";
+                out << tempType << "3 periodicDistance_delta = make_real3(";
                 for (int i = 0; i < 3; i++) {
                     if (i > 0)
                         out << ", ";
                     out << getTempName(node.getChildren()[i], temps) << "-" << getTempName(node.getChildren()[i+3], temps);
-                    out << ");\n";
                 }
+                out << ");\n";
                 out << "APPLY_PERIODIC_TO_DELTA(periodicDistance_delta)\n";
-                out << tempType << " periodicDistance_r = SQRT(periodicDistance_delta.x*periodicDistance_delta.x + periodicDistance_delta.y*periodicDistance_delta.y + periodicDistance_delta.z*periodicDistance_delta.z);\n";
+                out << tempType << " periodicDistance_rinv = RSQRT(periodicDistance_delta.x*periodicDistance_delta.x + periodicDistance_delta.y*periodicDistance_delta.y + periodicDistance_delta.z*periodicDistance_delta.z);\n";
                 for (int j = 0; j < nodes.size(); j++) {
                     const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder();
                     int argIndex = -1;
@@ -121,19 +121,19 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express
                         }
                     }
                     if (argIndex == -1)
-                        out << nodeNames[j] << " = periodicDistance_r;\n";
+                        out << nodeNames[j] << " = RECIP(periodicDistance_rinv);\n";
                     else if (argIndex == 0)
-                        out << nodeNames[j] << " = periodicDistance_delta.x/periodicDistance_r;\n";
+                        out << nodeNames[j] << " = periodicDistance_delta.x*periodicDistance_rinv;\n";
                     else if (argIndex == 1)
-                        out << nodeNames[j] << " = periodicDistance_delta.y/periodicDistance_r;\n";
+                        out << nodeNames[j] << " = periodicDistance_delta.y*periodicDistance_rinv;\n";
                     else if (argIndex == 2)
-                        out << nodeNames[j] << " = periodicDistance_delta.z/periodicDistance_r;\n";
+                        out << nodeNames[j] << " = periodicDistance_delta.z*periodicDistance_rinv;\n";
                     else if (argIndex == 3)
-                        out << nodeNames[j] << " = -periodicDistance_delta.x/periodicDistance_r;\n";
+                        out << nodeNames[j] << " = -periodicDistance_delta.x*periodicDistance_rinv;\n";
                     else if (argIndex == 4)
-                        out << nodeNames[j] << " = -periodicDistance_delta.y/periodicDistance_r;\n";
+                        out << nodeNames[j] << " = -periodicDistance_delta.y*periodicDistance_rinv;\n";
                     else if (argIndex == 5)
-                        out << nodeNames[j] << " = -periodicDistance_delta.z/periodicDistance_r;\n";
+                        out << nodeNames[j] << " = -periodicDistance_delta.z*periodicDistance_rinv;\n";
                 }
             }
             else {
@@ -775,4 +775,4 @@ Lepton::CustomFunction* CudaExpressionUtilities::getFunctionPlaceholder(const Ta
 
 Lepton::CustomFunction* CudaExpressionUtilities::getPeriodicDistancePlaceholder() {
     return &periodicDistance;
-}
+}
\ No newline at end of file

platforms/cuda/tests/TestCudaCustomExternalForce.cpp

@@ -196,8 +196,8 @@ void testPeriodic() {
 
         // 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);
+        ASSERT_EQUAL_VEC(delta*2, state.getForces()[0], 1e-5);
+        ASSERT_EQUAL_TOL(delta.dot(delta), state.getPotentialEnergy(), 1e-5);
         integrator.step(1);
     }
 }