Cleaned up constraint code and deleted obsolete constraint implementations

platforms/cuda/src/CudaKernels.cpp

@@ -399,7 +399,7 @@ static void initializeIntegration(const System& system, CudaPlatform::PlatformDa
         invMass1[i] = 1.0f/mass[particle1Index];
         invMass2[i] = 1.0f/mass[particle2Index];
     }
-    gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, (float)integrator.getConstraintTolerance(), 4);
+    gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, (float)integrator.getConstraintTolerance());
     
     // Finish initialization.
 
@@ -438,7 +438,7 @@ void CudaIntegrateVerletStepKernel::execute(OpenMMContextImpl& context, const Ve
     kVerletUpdatePart1(gpu);
     kApplyFirstShake(gpu);
     kApplyFirstSettle(gpu);
-    kApplyFirstCShake(gpu);
+    kApplyFirstCCMA(gpu);
     if (data.removeCM) {
         int step = (int) (context.getTime()/stepSize);
         if (step%data.cmMotionFrequency == 0)
@@ -477,7 +477,7 @@ void CudaIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, const
     kUpdatePart1(gpu);
     kApplyFirstShake(gpu);
     kApplyFirstSettle(gpu);
-    kApplyFirstCShake(gpu);
+    kApplyFirstCCMA(gpu);
     if (data.removeCM) {
         int step = (int) (context.getTime()/stepSize);
         if (step%data.cmMotionFrequency == 0)
@@ -486,7 +486,7 @@ void CudaIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, const
     kUpdatePart2(gpu);
     kApplySecondShake(gpu);
     kApplySecondSettle(gpu);
-    kApplySecondCShake(gpu);
+    kApplySecondCCMA(gpu);
 }
 
 CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
@@ -519,7 +519,7 @@ void CudaIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, const
     kBrownianUpdatePart1(gpu);
     kApplyFirstShake(gpu);
     kApplyFirstSettle(gpu);
-    kApplyFirstCShake(gpu);
+    kApplyFirstCCMA(gpu);
     if (data.removeCM) {
         int step = (int) (context.getTime()/stepSize);
         if (step%data.cmMotionFrequency == 0)

platforms/cuda/src/kernels/cudaKernels.h

@@ -42,14 +42,12 @@ extern void kCalculateAndersenThermostat(gpuContext gpu);
 extern void kReduceBornSumAndForces(gpuContext gpu);
 extern void kUpdatePart1(gpuContext gpu);
 extern void kApplyFirstShake(gpuContext gpu);
-extern void kApplyFirstCShake(gpuContext gpu);
+extern void kApplyFirstCCMA(gpuContext gpu);
 extern void kApplyFirstSettle(gpuContext gpu);
-extern void kApplyFirstLincs(gpuContext gpu);
 extern void kUpdatePart2(gpuContext gpu);
 extern void kApplySecondShake(gpuContext gpu);
-extern void kApplySecondCShake(gpuContext gpu);
+extern void kApplySecondCCMA(gpuContext gpu);
 extern void kApplySecondSettle(gpuContext gpu);
-extern void kApplySecondLincs(gpuContext gpu);
 extern void kVerletUpdatePart1(gpuContext gpu);
 extern void kVerletUpdatePart2(gpuContext gpu);
 extern void kBrownianUpdatePart1(gpuContext gpu);
@@ -78,10 +76,8 @@ extern void SetUpdateShakeHSim(gpuContext gpu);
 extern void GetUpdateShakeHSim(gpuContext gpu);
 extern void SetSettleSim(gpuContext gpu);
 extern void GetSettleSim(gpuContext gpu);
-extern void SetCShakeSim(gpuContext gpu);
-extern void GetCShakeSim(gpuContext gpu);
-extern void SetLincsSim(gpuContext gpu);
-extern void GetLincsSim(gpuContext gpu);
+extern void SetCCMASim(gpuContext gpu);
+extern void GetCCMASim(gpuContext gpu);
 extern void SetVerletUpdateSim(gpuContext gpu);
 extern void GetVerletUpdateSim(gpuContext gpu);
 extern void SetBrownianUpdateSim(gpuContext gpu);

platforms/cuda/src/kernels/cudatypes.h

@@ -258,7 +258,7 @@ struct cudaGmxSimulation {
     unsigned int    max_shake_threads_per_block;    // Maximum threads per block in shake kernel calls
     unsigned int    shake_threads_per_block;        // Threads per block in shake kernel calls
     unsigned int    settle_threads_per_block;       // Threads per block in SETTLE kernel calls
-    unsigned int    lincs_threads_per_block;        // Threads per block in LINCS kernel calls
+    unsigned int    ccma_threads_per_block;         // Threads per block in CCMA kernel calls
     unsigned int    nonshake_threads_per_block;     // Threads per block in nonshaking kernel call
     unsigned int    max_localForces_threads_per_block;  // Threads per block in local forces kernel calls
     unsigned int    localForces_threads_per_block;  // Threads per block in local forces kernel calls
@@ -360,7 +360,7 @@ struct cudaGmxSimulation {
     float           inverseTotalMass;               // Used in linear momentum removal
     unsigned int    ShakeConstraints;               // Total number of Shake constraints
     unsigned int    settleConstraints;              // Total number of Settle constraints
-    unsigned int    lincsConstraints;               // Total number of LINCS constraints.
+    unsigned int    ccmaConstraints;                // Total number of CCMA constraints.
     unsigned int    rigidClusters;                  // Total number of rigid clusters
     unsigned int    maxRigidClusterSize;            // The size of the largest rigid cluster
     unsigned int    clusterShakeBlockSize;          // The number of threads to process each rigid cluster
@@ -368,7 +368,6 @@ struct cudaGmxSimulation {
     unsigned int    maxShakeIterations;             // Maximum shake iterations
     unsigned int    degreesOfFreedom;               // Number of degrees of freedom in system
     float           shakeTolerance;                 // Shake tolerance
-    unsigned int    lincsTerms;                     // Number of terms in the matrix expansion for LINCS
     float           InvMassJ;                       // Shake inverse mass for hydrogens
     int*            pNonShakeID;                    // Not Shaking atoms
     int4*           pShakeID;                       // Shake atoms and phase
@@ -385,23 +384,15 @@ struct cudaGmxSimulation {
     int*            pAtomIndex;                     // The original index of each atom
     float4*         pGridBoundingBox;               // The size of each grid cell
     float4*         pGridCenter;                    // The center of each grid cell
-    int2*           pLincsAtoms;                    // The atoms connected by each LINCS constraint
-    float4*         pLincsDistance;                 // The displacement vector (x, y, z) and constraint distance (w) for each LINCS constraint
-    int*            pLincsConnections;              // The indices of constraints that other constraints are connected to
-    int*            pLincsNumConnections;           // The number of other constraints that each constraint is linked to
-    float*          pLincsS;                        // S matrix for LINCS
-    float*          pLincsCoupling;                 // Coupling matrix for LINCS
-    float*          pLincsRhs1;                     // Workspace for LINCS
-    float*          pLincsRhs2;                     // Workspace for LINCS
-    float*          pLincsSolution;                 // Workspace for LINCS
-    int*            pLincsAtomConstraints;          // The indices of constraints involving each atom
-    int*            pLincsNumAtomConstraints;       // The number of constraints involving each atom
+    int2*           pCcmaAtoms;                     // The atoms connected by each CCMA constraint
+    float4*         pCcmaDistance;                  // The displacement vector (x, y, z) and constraint distance (w) for each CCMA constraint
+    float*          pCcmaDelta1;                    // Workspace for CCMA
+    float*          pCcmaDelta2;                    // Workspace for CCMA
+    int*            pCcmaAtomConstraints;           // The indices of constraints involving each atom
+    int*            pCcmaNumAtomConstraints;        // The number of constraints involving each atom
     short*          pSyncCounter;                   // Used for global thread synchronization
-    unsigned int*   pRequiredIterations;            // Used by SHAKE to communicate whether iteration has converged
-    float*          pShakeReducedMass;              // The reduced mass for each SHAKE constraint
-    float*          pRigidClusterMatrix;            // The inverse constraint matrix for each rigid cluster
-    unsigned int*   pRigidClusterConstraintIndex;   // The index of each cluster in the stream containing cluster constraints.
-    unsigned int*   pRigidClusterMatrixIndex;       // The index of each cluster in the stream containing cluster matrices.
+    unsigned int*   pRequiredIterations;            // Used by CCMA to communicate whether iteration has converged
+    float*          pCcmaReducedMass;               // The reduced mass for each CCMA constraint
     unsigned int*   pConstraintMatrixColumn;        // The column of each element in the constraint matrix.
     float*          pConstraintMatrixValue;         // The value of each element in the constraint matrix.
 

platforms/cuda/src/kernels/gputypes.h

@@ -127,20 +127,15 @@ struct _gpuContext {
     CUDAStream<int>* psAtomIndex;           // The original index of each atom
     CUDAStream<float4>* psGridBoundingBox;  // The size of each grid cell
     CUDAStream<float4>* psGridCenter;       // The center and radius for each grid cell
-    CUDAStream<int2>* psLincsAtoms;         // The atoms connected by each LINCS constraint
-    CUDAStream<float4>* psLincsDistance;    // The displacement vector (x, y, z) and constraint distance (w) for each LINCS constraint
-    CUDAStream<int>* psLincsConnections;    // The indices of constraints that other constraints are connected to
-    CUDAStream<int>* psLincsNumConnections; // The number of other constraints that each constraint is linked to
-    CUDAStream<int>* psLincsAtomConstraints; // The indices of constraints involving each atom
-    CUDAStream<int>* psLincsNumAtomConstraints; // The number of constraints involving each atom
-    CUDAStream<float>* psLincsS;            // S matrix for LINCS
-    CUDAStream<float>* psLincsCoupling;     // Coupling matrix for LINCS
-    CUDAStream<float>* psLincsRhs1;         // Workspace for LINCS
-    CUDAStream<float>* psLincsRhs2;         // Workspace for LINCS
-    CUDAStream<float>* psLincsSolution;     // Workspace for LINCS
+    CUDAStream<int2>* psCcmaAtoms;          // The atoms connected by each CCMA constraint
+    CUDAStream<float4>* psCcmaDistance;     // The displacement vector (x, y, z) and constraint distance (w) for each CCMA constraint
+    CUDAStream<int>* psCcmaAtomConstraints; // The indices of constraints involving each atom
+    CUDAStream<int>* psCcmaNumAtomConstraints; // The number of constraints involving each atom
+    CUDAStream<float>* psCcmaDelta1;        // Workspace for CCMA
+    CUDAStream<float>* psCcmaDelta2;        // Workspace for CCMA
     CUDAStream<short>* psSyncCounter;       // Used for global thread synchronization
     CUDAStream<unsigned int>* psRequiredIterations; // Used by SHAKE to communicate whether iteration has converged
-    CUDAStream<float>* psShakeReducedMass;  // The reduced mass for each SHAKE constraint
+    CUDAStream<float>* psCcmaReducedMass;   // The reduced mass for each SHAKE constraint
     CUDAStream<float>* psRigidClusterMatrix;// The inverse constraint matrix for each rigid cluster
     CUDAStream<unsigned int>* psRigidClusterConstraintIndex; // The index of each cluster in the stream containing cluster constraints.
     CUDAStream<unsigned int>* psRigidClusterMatrixIndex; // The index of each cluster in the stream containing cluster matrices.
@@ -192,7 +187,7 @@ void gpuSetObcParameters(gpuContext gpu, float innerDielectric, float solventDie
 
 extern "C"
 void gpuSetConstraintParameters(gpuContext gpu, const std::vector<int>& atom1, const std::vector<int>& atom2, const std::vector<float>& distance,
-        const std::vector<float>& invMass1, const std::vector<float>& invMass2, float shakeTolerance, unsigned int lincsTerms);
+        const std::vector<float>& invMass1, const std::vector<float>& invMass2, float constraintTolerance);
 
 extern "C"
 int gpuAllocateInitialBuffers(gpuContext gpu);

platforms/cuda/src/kernels/kCShake.cu → platforms/cuda/src/kernels/kCCMA.cu

@@ -34,14 +34,14 @@ using namespace std;
 
 static __constant__ cudaGmxSimulation cSim;
 
-void SetCShakeSim(gpuContext gpu)
+void SetCCMASim(gpuContext gpu)
 {
     cudaError_t status;
     status = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));
     RTERROR(status, "cudaMemcpyToSymbol: SetSim copy to cSim failed");
 }
 
-void GetCShakeSim(gpuContext gpu)
+void GetCCMASim(gpuContext gpu)
 {
     cudaError_t status;
     status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));
@@ -66,7 +66,7 @@ __device__ void kSyncAllThreads_kernel(short* syncCounter, short newCount)
     __syncthreads();
 }
 
-__global__ void kApplyCShake_kernel(float4* atomPositions, bool addOldPosition)
+__global__ void kApplyCCMA_kernel(float4* atomPositions, bool addOldPosition)
 {
     // Initialize counters used for monitoring convergence and doing global thread synchronization.
 
@@ -82,16 +82,16 @@ __global__ void kApplyCShake_kernel(float4* atomPositions, bool addOldPosition)
     // Calculate the direction of each constraint.
 
     unsigned int pos = threadIdx.x + blockIdx.x * blockDim.x;
-    while (pos < cSim.lincsConstraints)
+    while (pos < cSim.ccmaConstraints)
     {
-        int2 atoms = cSim.pLincsAtoms[pos];
-        float4 dir = cSim.pLincsDistance[pos];
+        int2 atoms = cSim.pCcmaAtoms[pos];
+        float4 dir = cSim.pCcmaDistance[pos];
         float4 oldPos1 = cSim.pOldPosq[atoms.x];
         float4 oldPos2 = cSim.pOldPosq[atoms.y];
         dir.x = oldPos1.x-oldPos2.x;
         dir.y = oldPos1.y-oldPos2.y;
         dir.z = oldPos1.z-oldPos2.z;
-        cSim.pLincsDistance[pos] = dir;
+        cSim.pCcmaDistance[pos] = dir;
         pos += blockDim.x*gridDim.x;
     }
     __syncthreads();
@@ -106,12 +106,12 @@ __global__ void kApplyCShake_kernel(float4* atomPositions, bool addOldPosition)
         // Calculate the constraint force for each constraint.
 
         pos = threadIdx.x + blockIdx.x * blockDim.x;
-        while (pos < cSim.lincsConstraints)
+        while (pos < cSim.ccmaConstraints)
         {
-            int2 atoms = cSim.pLincsAtoms[pos];
+            int2 atoms = cSim.pCcmaAtoms[pos];
             float4 delta1 = atomPositions[atoms.x];
             float4 delta2 = atomPositions[atoms.y];
-            float4 dir = cSim.pLincsDistance[pos];
+            float4 dir = cSim.pCcmaDistance[pos];
             float3 rp_ij = make_float3(delta1.x-delta2.x, delta1.y-delta2.y, delta1.z-delta2.z);
             if (addOldPosition)
             {
@@ -124,8 +124,8 @@ __global__ void kApplyCShake_kernel(float4* atomPositions, bool addOldPosition)
             float diff = dist2 - rp2;
             float rrpr  = rp_ij.x*dir.x + rp_ij.y*dir.y + rp_ij.z*dir.z;
             float d_ij2  = dir.x*dir.x + dir.y*dir.y + dir.z*dir.z;
-            float reducedMass = cSim.pShakeReducedMass[pos];
-            cSim.pLincsRhs1[pos] = (rrpr > d_ij2*1e-6f ? reducedMass*diff/rrpr : 0.0f);
+            float reducedMass = cSim.pCcmaReducedMass[pos];
+            cSim.pCcmaDelta1[pos] = (rrpr > d_ij2*1e-6f ? reducedMass*diff/rrpr : 0.0f);
             if (requiredIterations == iteration && (rp2 < lowerTol*dist2 || rp2 > upperTol*dist2))
                 requiredIterations = iteration+1;
             pos += blockDim.x * gridDim.x;
@@ -140,18 +140,18 @@ __global__ void kApplyCShake_kernel(float4* atomPositions, bool addOldPosition)
         // Multiply by the inverse constraint matrix.
 
         pos = threadIdx.x + blockIdx.x * blockDim.x;
-        while (pos < cSim.lincsConstraints)
+        while (pos < cSim.ccmaConstraints)
         {
             float sum = 0.0f;
             for (unsigned int i = 0; ; i++)
             {
-                unsigned int index = pos+i*cSim.lincsConstraints;
+                unsigned int index = pos+i*cSim.ccmaConstraints;
                 unsigned int column = cSim.pConstraintMatrixColumn[index];
-                if (column >= cSim.lincsConstraints)
+                if (column >= cSim.ccmaConstraints)
                     break;
-                sum += cSim.pLincsRhs1[column]*cSim.pConstraintMatrixValue[index];
+                sum += cSim.pCcmaDelta1[column]*cSim.pConstraintMatrixValue[index];
             }
-            cSim.pLincsSolution[pos] = sum;
+            cSim.pCcmaDelta2[pos] = sum;
             pos += blockDim.x * gridDim.x;
         }
         kSyncAllThreads_kernel(&cSim.pSyncCounter[gridDim.x], iteration);
@@ -164,16 +164,16 @@ __global__ void kApplyCShake_kernel(float4* atomPositions, bool addOldPosition)
         {
             float4 atomPos = atomPositions[pos];
             float invMass = cSim.pVelm4[pos].w;
-            int num = cSim.pLincsNumAtomConstraints[pos];
+            int num = cSim.pCcmaNumAtomConstraints[pos];
             for (int i = 0; i < num; i++)
             {
                 int index = pos+i*cSim.atoms;
-                int constraint = cSim.pLincsAtomConstraints[index];
+                int constraint = cSim.pCcmaAtomConstraints[index];
                 bool forward = (constraint > 0);
                 constraint = (forward ? constraint-1 : -constraint-1);
-                float constraintForce = damping*invMass*cSim.pLincsSolution[constraint];
+                float constraintForce = damping*invMass*cSim.pCcmaDelta2[constraint];
                 constraintForce = (forward ? constraintForce : -constraintForce);
-                float4 dir = cSim.pLincsDistance[constraint];
+                float4 dir = cSim.pCcmaDistance[constraint];
                 atomPos.x += constraintForce*dir.x;
                 atomPos.y += constraintForce*dir.y;
                 atomPos.z += constraintForce*dir.z;
@@ -192,22 +192,22 @@ __global__ void kApplyCShake_kernel(float4* atomPositions, bool addOldPosition)
         cSim.pSyncCounter[blockIdx.x] = -1;
 }
 
-void kApplyFirstCShake(gpuContext gpu)
+void kApplyFirstCCMA(gpuContext gpu)
 {
-//    printf("kApplyFirstCShake\n");
-    if (gpu->sim.lincsConstraints > 0)
+//    printf("kApplyFirstCCMA\n");
+    if (gpu->sim.ccmaConstraints > 0)
     {
-        kApplyCShake_kernel<<<gpu->sim.blocks, gpu->sim.lincs_threads_per_block>>>(gpu->sim.pPosqP, true);
-        LAUNCHERROR("kApplyCShake");
+        kApplyCCMA_kernel<<<gpu->sim.blocks, gpu->sim.ccma_threads_per_block>>>(gpu->sim.pPosqP, true);
+        LAUNCHERROR("kApplyCCMA");
     }
 }
 
-void kApplySecondCShake(gpuContext gpu)
+void kApplySecondCCMA(gpuContext gpu)
 {
-//    printf("kApplySecondCShake\n");
-    if (gpu->sim.lincsConstraints > 0)
+//    printf("kApplySecondCCMA\n");
+    if (gpu->sim.ccmaConstraints > 0)
     {
-        kApplyCShake_kernel<<<gpu->sim.blocks, gpu->sim.lincs_threads_per_block>>>(gpu->sim.pPosq, false);
-        LAUNCHERROR("kApplyCShake");
+        kApplyCCMA_kernel<<<gpu->sim.blocks, gpu->sim.ccma_threads_per_block>>>(gpu->sim.pPosq, false);
+        LAUNCHERROR("kApplyCCMA");
     }
 }

platforms/cuda/src/kernels/kLincs.cu

-/* -------------------------------------------------------------------------- *
- *                                   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) 2009 Stanford University and the Authors.           *
- * Authors: Scott Le Grand, Peter Eastman                                     *
- * Contributors:                                                              *
- *                                                                            *
- * This program is free software: you can redistribute it and/or modify       *
- * it under the terms of the GNU Lesser General Public License as published   *
- * by the Free Software Foundation, either version 3 of the License, or       *
- * (at your option) any later version.                                        *
- *                                                                            *
- * This program is distributed in the hope that it will be useful,            *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
- * GNU Lesser General Public License for more details.                        *
- *                                                                            *
- * You should have received a copy of the GNU Lesser General Public License   *
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
- * -------------------------------------------------------------------------- */
-
-#include <cuda.h>
-#include <vector_functions.h>
-using namespace std;
-
-#include "gputypes.h"
-
-
-static __constant__ cudaGmxSimulation cSim;
-
-void SetLincsSim(gpuContext gpu)
-{
-    cudaError_t status;
-    status = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));
-    RTERROR(status, "cudaMemcpyToSymbol: SetSim copy to cSim failed");
-}
-
-void GetLincsSim(gpuContext gpu)
-{
-    cudaError_t status;
-    status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));
-    RTERROR(status, "cudaMemcpyFromSymbol: SetSim copy from cSim failed");
-}
-
-/**
- * Synchronize all threads across all blocks.
- */
-__device__ void kSyncAllThreads_kernel(short* syncCounter, short newCount)
-{
-//    short* syncCounter = &cSim.pSyncCounter[newCount%2 == 0 ? 0 : gridDim.x];
-    __syncthreads();
-    if (threadIdx.x == 0)
-        syncCounter[blockIdx.x] = newCount;
-    if (threadIdx.x < gridDim.x)
-    {
-        volatile short* counter = &syncCounter[threadIdx.x];
-        do
-        {
-        } while (*counter != newCount);
-    }
-    __syncthreads();
-}
-
-__global__ void kSolveLincsMatrix_kernel(float4* atomPositions)
-{
-    for (unsigned int iteration = 0; iteration < cSim.lincsTerms; iteration++) {
-        float* rhs1 = (iteration%2 == 0 ? cSim.pLincsRhs1 : cSim.pLincsRhs2);
-        float* rhs2 = (iteration%2 == 0 ? cSim.pLincsRhs2 : cSim.pLincsRhs1);
-        unsigned int pos = threadIdx.x + blockIdx.x * blockDim.x;
-        while (pos < cSim.lincsConstraints)
-        {
-            float rhs = 0.0f;
-            int num = cSim.pLincsNumConnections[pos];
-            for (int i = 0; i < num; i++)
-            {
-                int index = pos+i*cSim.lincsConstraints;
-                int otherConstraint = cSim.pLincsConnections[index];
-                rhs += cSim.pLincsCoupling[index]*rhs1[otherConstraint];
-            }
-            rhs2[pos] = rhs;
-            cSim.pLincsSolution[pos] += rhs;
-            pos += blockDim.x * gridDim.x;
-        }
-        kSyncAllThreads_kernel(&cSim.pSyncCounter[iteration%2 == 0 ? 0 : gridDim.x], iteration);
-    }
-
-    // Update the atom positions based on the solution to the matrix equations.
-
-    unsigned int pos = threadIdx.x + blockIdx.x * blockDim.x;
-    while (pos < cSim.atoms)
-    {
-        float4 atomPos = atomPositions[pos];
-        float invMass = cSim.pVelm4[pos].w;
-        int num = cSim.pLincsNumAtomConstraints[pos];
-        for (int i = 0; i < num; i++)
-        {
-            int index = pos+i*cSim.atoms;
-            int constraint = cSim.pLincsAtomConstraints[index];
-            bool forward = (constraint > 0);
-            constraint = (forward ? constraint-1 : -constraint-1);
-            float4 dir = cSim.pLincsDistance[constraint];
-            float c = invMass*cSim.pLincsS[constraint]*cSim.pLincsSolution[constraint];
-            c = (forward ? -c : c);
-            atomPos.x += c*dir.x;
-            atomPos.y += c*dir.y;
-            atomPos.z += c*dir.z;
-        }
-        atomPositions[pos] = atomPos;
-        pos += blockDim.x * gridDim.x;
-    }
-}
-
-__global__ void kApplyLincsPart1_kernel(float4* atomPositions, bool addOldPosition)
-{
-   // Calculate the direction of each constraint, along with the initial RHS and solution vectors.
-
-    unsigned int pos = threadIdx.x + blockIdx.x * blockDim.x;
-    while (pos < cSim.lincsConstraints)
-    {
-        int2 atoms = cSim.pLincsAtoms[pos];
-        float4 delta1 = atomPositions[atoms.x];
-        float4 delta2 = atomPositions[atoms.y];
-        float4 dir = cSim.pLincsDistance[pos];
-        if (addOldPosition)
-        {
-            float4 oldPos1 = cSim.pOldPosq[atoms.x];
-            float4 oldPos2 = cSim.pOldPosq[atoms.y];
-            dir.x = (oldPos1.x-oldPos2.x)+(delta1.x-delta2.x);
-            dir.y = (oldPos1.y-oldPos2.y)+(delta1.y-delta2.y);
-            dir.z = (oldPos1.z-oldPos2.z)+(delta1.z-delta2.z);
-        }
-        else
-        {
-            dir.x = delta1.x-delta2.x;
-            dir.y = delta1.y-delta2.y;
-            dir.z = delta1.z-delta2.z;
-        }
-        float invLength = 1.0f/sqrt(dir.x*dir.x + dir.y*dir.y + dir.z*dir.z);
-        dir.x *= invLength;
-        dir.y *= invLength;
-        dir.z *= invLength;
-        cSim.pLincsDistance[pos] = dir;
-        float diff = cSim.pLincsS[pos]*(1.0f/invLength-dir.w);
-        cSim.pLincsRhs1[pos] = diff;
-        cSim.pLincsSolution[pos] = diff;
-        pos += blockDim.x * gridDim.x;
-    }
-    kSyncAllThreads_kernel(cSim.pSyncCounter, cSim.lincsTerms+1);
-
-    // Build the coupling matrix.
-
-    pos = threadIdx.x + blockIdx.x * blockDim.x;
-    while (pos < cSim.lincsConstraints)
-    {
-        float4 dir1 = cSim.pLincsDistance[pos];
-        int2 atoms1 = cSim.pLincsAtoms[pos];
-        int num = cSim.pLincsNumConnections[pos];
-        float s = cSim.pLincsS[pos];
-        float invMass = cSim.pVelm4[atoms1.x].w;
-        for (int i = 0; i < num; i++)
-        {
-            int index = pos+i*cSim.lincsConstraints;
-            int otherConstraint = cSim.pLincsConnections[index];
-            float4 dir2 = cSim.pLincsDistance[otherConstraint];
-            int2 atoms2 = cSim.pLincsAtoms[otherConstraint];
-            float signedMass = (atoms1.x == atoms2.x || atoms1.y == atoms2.y ? -invMass : cSim.pVelm4[atoms1.y].w);
-            cSim.pLincsCoupling[index] = signedMass*s*(dir1.x*dir2.x+dir1.y*dir2.y+dir1.z*dir2.z)*cSim.pLincsS[otherConstraint];
-        }
-        pos += blockDim.x * gridDim.x;
-    }
-}
-
-__global__ void kApplyLincsPart2_kernel(float4* atomPositions, bool addOldPosition)
-{
-    // Correct for rotational lengthening.
-
-    unsigned int pos = threadIdx.x + blockIdx.x * blockDim.x;
-    while (pos < cSim.lincsConstraints)
-    {
-        int2 atoms = cSim.pLincsAtoms[pos];
-        float4 delta1 = atomPositions[atoms.x];
-        float4 delta2 = atomPositions[atoms.y];
-        float3 delta;
-        if (addOldPosition)
-        {
-            float4 oldPos1 = cSim.pOldPosq[atoms.x];
-            float4 oldPos2 = cSim.pOldPosq[atoms.y];
-            delta = make_float3((oldPos1.x-oldPos2.x)+(delta1.x-delta2.x),
-                                (oldPos1.y-oldPos2.y)+(delta1.y-delta2.y),
-                                (oldPos1.z-oldPos2.z)+(delta1.z-delta2.z));
-        }
-        else
-        {
-            delta = make_float3(delta1.x-delta2.x, delta1.y-delta2.y, delta1.z-delta2.z);
-        }
-        float distance = cSim.pLincsDistance[pos].w;
-        float p2 = 2.0f*distance*distance-(delta.x*delta.x+delta.y*delta.y+delta.z*delta.z);
-        p2 = (p2 < 0.0f ? 0.0f : p2);
-        float diff = cSim.pLincsS[pos]*(distance-sqrt(p2));
-        cSim.pLincsRhs1[pos] = diff;
-        cSim.pLincsSolution[pos] = diff;
-        pos += blockDim.x * gridDim.x;
-    }
-}
-
-static void kApplyLincs(gpuContext gpu, float4* atomPositions, bool addOldPosition)
-{
-    kApplyLincsPart1_kernel<<<gpu->sim.blocks, gpu->sim.lincs_threads_per_block>>>(atomPositions, addOldPosition);
-    LAUNCHERROR("kApplyLincsPart1");
-    kSolveLincsMatrix_kernel<<<gpu->sim.blocks, gpu->sim.lincs_threads_per_block>>>(atomPositions);
-    LAUNCHERROR("kSolveLincsMatrix_kernel");
-    kApplyLincsPart2_kernel<<<gpu->sim.blocks, gpu->sim.lincs_threads_per_block>>>(atomPositions, addOldPosition);
-    LAUNCHERROR("kApplyLincsPart2");
-    kSolveLincsMatrix_kernel<<<gpu->sim.blocks, gpu->sim.lincs_threads_per_block>>>(atomPositions);
-    LAUNCHERROR("kSolveLincsMatrix_kernel");
-}
-
-void kApplyFirstLincs(gpuContext gpu)
-{
-//    printf("kApplyFirstLincs\n");
-    if (gpu->sim.lincsConstraints > 0)
-        kApplyLincs(gpu, gpu->sim.pPosqP, true);
-}
-
-void kApplySecondLincs(gpuContext gpu)
-{
-//    printf("kApplySecondLincs\n");
-    if (gpu->sim.lincsConstraints > 0)
-        kApplyLincs(gpu, gpu->sim.pPosq, false);
-}

platforms/reference/src/ReferenceKernels.cpp

@@ -38,14 +38,12 @@
 #include "SimTKReference/ReferenceBondForce.h"
 #include "SimTKReference/ReferenceBrownianDynamics.h"
 #include "SimTKReference/ReferenceHarmonicBondIxn.h"
-#include "SimTKReference/ReferenceLincsAlgorithm.h"
 #include "SimTKReference/ReferenceLJCoulomb14.h"
 #include "SimTKReference/ReferenceLJCoulombIxn.h"
 #include "SimTKReference/ReferenceProperDihedralBond.h"
 #include "SimTKReference/ReferenceRbDihedralBond.h"
 #include "SimTKReference/ReferenceStochasticDynamics.h"
-#include "SimTKReference/ReferenceRigidShakeAlgorithm.h"
-#include "SimTKReference/ReferenceShakeAlgorithm.h"
+#include "SimTKReference/ReferenceCCMAAlgorithm.h"
 #include "SimTKReference/ReferenceVerletDynamics.h"
 #include "openmm/CMMotionRemover.h"
 #include "openmm/System.h"
@@ -108,7 +106,7 @@ static void disposeRealArray(RealOpenMM** array, int size) {
     }
 }
 
-static void findAnglesForShake(const System& system, vector<ReferenceRigidShakeAlgorithm::AngleInfo>& angles) {
+static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorithm::AngleInfo>& angles) {
     for (int i = 0; i < system.getNumForces(); i++) {
         const HarmonicAngleForce* force = dynamic_cast<const HarmonicAngleForce*>(&system.getForce(i));
         if (force != NULL) {
@@ -116,7 +114,7 @@ static void findAnglesForShake(const System& system, vector<ReferenceRigidShakeA
                 int atom1, atom2, atom3;
                 double angle, k;
                 force->getAngleParameters(j, atom1, atom2, atom3, angle, k);
-                angles.push_back(ReferenceRigidShakeAlgorithm::AngleInfo(atom1, atom2, atom3, angle));
+                angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, angle));
             }
         }
     }
@@ -623,9 +621,9 @@ void ReferenceIntegrateVerletStepKernel::execute(OpenMMContextImpl& context, con
             delete constraints;
         }
         dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) );
-        vector<ReferenceRigidShakeAlgorithm::AngleInfo> angles;
-        findAnglesForShake(context.getSystem(), angles);
-        constraints = new ReferenceRigidShakeAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
+        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
+        findAnglesForCCMA(context.getSystem(), angles);
+        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
         dynamics->setReferenceConstraintAlgorithm(constraints);
         prevStepSize = stepSize;
     }
@@ -684,9 +682,9 @@ void ReferenceIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, c
 				static_cast<RealOpenMM>(stepSize), 
 				static_cast<RealOpenMM>(tau), 
 				static_cast<RealOpenMM>(temperature) );
-        vector<ReferenceRigidShakeAlgorithm::AngleInfo> angles;
-        findAnglesForShake(context.getSystem(), angles);
-        constraints = new ReferenceRigidShakeAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
+        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
+        findAnglesForCCMA(context.getSystem(), angles);
+        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
         dynamics->setReferenceConstraintAlgorithm(constraints);
         prevTemp = temperature;
         prevFriction = friction;
@@ -746,9 +744,9 @@ void ReferenceIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, c
 				static_cast<RealOpenMM>(stepSize), 
 				static_cast<RealOpenMM>(friction), 
 				static_cast<RealOpenMM>(temperature) );
-        vector<ReferenceRigidShakeAlgorithm::AngleInfo> angles;
-        findAnglesForShake(context.getSystem(), angles);
-        constraints = new ReferenceRigidShakeAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
+        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
+        findAnglesForCCMA(context.getSystem(), angles);
+        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
         dynamics->setReferenceConstraintAlgorithm(constraints);
         prevTemp = temperature;
         prevFriction = friction;

platforms/reference/src/SimTKReference/ReferenceRigidShakeAlgorithm.cpp → platforms/reference/src/SimTKReference/ReferenceCCMAAlgorithm.cpp

@@ -28,9 +28,8 @@
 #include "../SimTKUtilities/SimTKOpenMMCommon.h"
 #include "../SimTKUtilities/SimTKOpenMMUtilities.h"
 #include "../SimTKUtilities/SimTKOpenMMLog.h"
-#include "ReferenceRigidShakeAlgorithm.h"
+#include "ReferenceCCMAAlgorithm.h"
 #include "ReferenceDynamics.h"
-#include "jama_svd.h"
 #include "quern.h"
 #include "openmm/Vec3.h"
 #include <map>
@@ -40,22 +39,22 @@ using std::pair;
 using std::vector;
 using std::set;
 using OpenMM::Vec3;
-using TNT::Array2D;
-using JAMA::SVD;
 
 /**---------------------------------------------------------------------------------------
 
-   ReferenceQShakeAlgorithm constructor
+   ReferenceCCMAAlgorithm constructor
 
          @param numberOfAtoms    number of atoms
          @param numberOfConstraints      number of constraints
-   @param atomIndices              block of atom indices
-   @param shakeParameters          Shake parameters
+         @param atomIndices              atom indices for contraints
+         @param distance                 distances for constraints
+         @param masses                   atom masses
+         @param angles                   angle force field terms
          @param tolerance                constraint tolerance
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceRigidShakeAlgorithm::ReferenceRigidShakeAlgorithm( int numberOfAtoms,
+ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm( int numberOfAtoms,
                                                   int numberOfConstraints,
                                                   int** atomIndices,
                                                   RealOpenMM* distance,
@@ -65,7 +64,7 @@ ReferenceRigidShakeAlgorithm::ReferenceRigidShakeAlgorithm( int numberOfAtoms,
 
    // ---------------------------------------------------------------------------------------
 
-   // static const char* methodName = "\nReferenceQShakeAlgorithm::ReferenceQShakeAlgorithm";
+   // static const char* methodName = "\nReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm";
 
    static const RealOpenMM zero        =  0.0;
    static const RealOpenMM one         =  1.0;
@@ -201,181 +200,26 @@ ReferenceRigidShakeAlgorithm::ReferenceRigidShakeAlgorithm( int numberOfAtoms,
            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 (abs(value) > 0.02)
+               if (FABS(value) > 0.02)
                    _matrix[j].push_back(pair<int, RealOpenMM>(i, (RealOpenMM) value));
            }
        }
        QUERN_free_result(qRowStart, qColIndex, qValue);
        QUERN_free_result(rRowStart, rColIndex, rValue);
    }
-
-
-
-
-
-
-//   // Identify rigid clusters.
-//
-//   vector<map<int, int> > atomConstraints(numberOfAtoms);
-//   for (int i = 0; i < numberOfConstraints; i++) {
-//       atomConstraints[atomIndices[i][0]][atomIndices[i][1]] = i;
-//       atomConstraints[atomIndices[i][1]][atomIndices[i][0]] = i;
-//   }
-//   for (int i = 0; i < numberOfAtoms; i++) {
-//       if (atomConstraints[i].size() < 2)
-//           continue;
-//
-//       // Begin by looking for a triangle this atom is part of.
-//
-//       set<int> atoms;
-//       atoms.insert(i);
-//       for (map<int, int>::const_iterator atom1 = atomConstraints[i].begin(); atom1 != atomConstraints[i].end() && atoms.size() == 1; ++atom1) {
-//           for (map<int, int>::const_iterator atom2 = atomConstraints[atom1->first].begin(); atom2 != atomConstraints[atom1->first].end(); ++atom2) {
-//               if (atomConstraints[i].count(atom2->first) != 0) {
-//                   atoms.insert(atom1->first);
-//                   atoms.insert(atom2->first);
-//                   break;
-//               }
-//           }
-//       }
-//       if (atoms.size() == 1)
-//           continue;
-//
-//       // We have three atoms that are part of a cluster, so look for other atoms we can add.
-//
-//       bool done = false;
-//       while (!done) {
-//           done = true;
-//           for (set<int>::const_iterator atom1 = atoms.begin(); atom1 != atoms.end(); ++atom1) {
-//               for (map<int, int>::const_iterator atom2 = atomConstraints[*atom1].begin(); atom2 != atomConstraints[*atom1].end(); ++atom2) {
-//                   if (atoms.find(atom2->first) != atoms.end())
-//                       continue; // This atom is already in the cluster.
-//
-//                   // See if this atom is linked to three other atoms in the cluster.
-//
-//                   int linkCount = 0;
-//                   for (map<int, int>::const_iterator atom3 = atomConstraints[atom2->first].begin(); atom3 != atomConstraints[atom2->first].end(); ++atom3)
-//                       if (atoms.find(atom3->first) != atoms.end())
-//                           linkCount++;
-//                   if (linkCount > 2) {
-//                       atoms.insert(atom2->first);
-//                       done = false;
-//                   }
-//               }
-//           }
-//       }
-//
-//       // Record the cluster.
-//
-//       vector<int> constraints;
-//       for (set<int>::const_iterator atom1 = atoms.begin(); atom1 != atoms.end(); ++atom1) {
-//           for (map<int, int>::const_iterator atom2 = atomConstraints[*atom1].begin(); atom2 != atomConstraints[*atom1].end(); ++atom2) {
-//               if (*atom1 < atom2->first && atoms.find(atom2->first) != atoms.end())
-//                   constraints.push_back(atom2->second);
-//           }
-//       }
-//       _rigidClusters.push_back(constraints);
-//       for (set<int>::const_iterator atom1 = atoms.begin(); atom1 != atoms.end(); ++atom1) {
-//           for (map<int, int>::const_iterator atom2 = atomConstraints[*atom1].begin(); atom2 != atomConstraints[*atom1].end(); ++atom2)
-//               atomConstraints[atom2->first].erase(*atom1);
-//           atomConstraints[*atom1].clear();
-//       }
-//
-//       // Compute the constraint coupling matrix for this cluster.
-//
-//       const vector<int>& cluster = _rigidClusters[i];
-//       unsigned int size = cluster.size();
-//       Array2D<double> matrix(size, size);
-//       for (int j = 0; j < (int)size; j++) {
-//           for (int k = 0; k < (int)size; k++) {
-//               double scale;
-//               int atomj0 = _atomIndices[cluster[j]][0];
-//               int atomj1 = _atomIndices[cluster[j]][1];
-//               int atomk0 = _atomIndices[cluster[k]][0];
-//               int atomk1 = _atomIndices[cluster[k]][1];
-//               RealOpenMM invMass0 = one/masses[atomj0];
-//               RealOpenMM invMass1 = one/masses[atomj1];
-//               int atoma, atomb;
-//               if (atomj0 == atomk0) {
-//                   atoma = atomj1;
-//                   atomb = atomk1;
-//                   scale = invMass0/(invMass0+invMass1);
-//               }
-//               else if (atomj1 == atomk1) {
-//                   atoma = atomj0;
-//                   atomb = atomk0;
-//                   scale = invMass1/(invMass0+invMass1);
-//               }
-//               else if (atomj0 == atomk1) {
-//                   atoma = atomj1;
-//                   atomb = atomk0;
-//                   scale = invMass0/(invMass0+invMass1);
-//               }
-//               else if (atomj1 == atomk0) {
-//                   atoma = atomj0;
-//                   atomb = atomk1;
-//                   scale = invMass1/(invMass0+invMass1);
-//               }
-//               else {
-//                   matrix[j][k] = 0.0;
-//                   continue; // These constraints are not connected.
-//               }
-//
-//               // Find the third constraint forming a triangle with these two.
-//
-//               for (int m = 0; m < size; m++) {
-//                   int other = cluster[m];
-//                   if ((_atomIndices[other][0] == atoma && _atomIndices[other][1] == atomb) || (_atomIndices[other][0] == atomb && _atomIndices[other][1] == atoma)) {
-//                       double d1 = _distance[cluster[j]];
-//                       double d2 = _distance[cluster[k]];
-//                       double d3 = _distance[other];
-//                       matrix[j][k] = scale*(d1*d1+d2*d2-d3*d3)/(2.0*d1*d2);
-//                       break;
-//                   }
-//               }
-//           }
-//           matrix[j][j] = 1.0;
-//       }
-//
-//       // Invert it using SVD.
-//
-//       Array2D<double> u, v;
-//       Array1D<double> w;
-//       SVD<double> svd(matrix);
-//       svd.getU(u);
-//       svd.getV(v);
-//       svd.getSingularValues(w);
-//       double singularValueCutoff = 0.01*w[0];
-//       for (int j = 0; j < (int)size; j++)
-//           w[j] = (w[j] < singularValueCutoff ? 0.0 : 1.0/w[j]);
-//       for (int j = 0; j < (int)size; j++) {
-//           for (int k = 0; k < (int)size; k++) {
-//               matrix[j][k] = 0.0;
-//               for (int m = 0; m < (int)size; m++)
-//                   matrix[j][k] += v[j][m]*w[m]*u[k][m];
-//           }
-//       }
-//
-//       // Record the inverted matrix.
-//
-//       _matrices.push_back(SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray(constraints.size(), constraints.size(), NULL, false, 0.0, ""));
-//       for (int j = 0; j < (int)size; j++)
-//           for (int k = 0; k < (int)size; k++)
-//               _matrices[i][j][k] = (RealOpenMM)matrix[j][k]*_distance[cluster[k]]/_distance[cluster[j]];
-//   }
 }
 
 /**---------------------------------------------------------------------------------------
 
-   ReferenceQShakeAlgorithm destructor
+   ReferenceCCMAAlgorithm destructor
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceRigidShakeAlgorithm::~ReferenceRigidShakeAlgorithm( ){
+ReferenceCCMAAlgorithm::~ReferenceCCMAAlgorithm( ){
 
    // ---------------------------------------------------------------------------------------
 
-   // static const char* methodName = "\nReferenceQShakeAlgorithm::~ReferenceQShakeAlgorithm";
+   // static const char* methodName = "\nReferenceCCMAAlgorithm::~ReferenceCCMAAlgorithm";
 
    // ---------------------------------------------------------------------------------------
 
@@ -398,11 +242,11 @@ ReferenceRigidShakeAlgorithm::~ReferenceRigidShakeAlgorithm( ){
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceRigidShakeAlgorithm::getNumberOfConstraints( void ) const {
+int ReferenceCCMAAlgorithm::getNumberOfConstraints( void ) const {
 
    // ---------------------------------------------------------------------------------------
 
-   // static const char* methodName = "\nReferenceQShakeAlgorithm::getNumberOfConstraints";
+   // static const char* methodName = "\nReferenceCCMAAlgorithm::getNumberOfConstraints";
 
    // ---------------------------------------------------------------------------------------
 
@@ -417,11 +261,11 @@ int ReferenceRigidShakeAlgorithm::getNumberOfConstraints( void ) const {
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceRigidShakeAlgorithm::getMaximumNumberOfIterations( void ) const {
+int ReferenceCCMAAlgorithm::getMaximumNumberOfIterations( void ) const {
 
    // ---------------------------------------------------------------------------------------
 
-   // static const char* methodName = "\nReferenceQShakeAlgorithm::getMaximumNumberOfIterations";
+   // static const char* methodName = "\nReferenceCCMAAlgorithm::getMaximumNumberOfIterations";
 
    // ---------------------------------------------------------------------------------------
 
@@ -438,11 +282,11 @@ int ReferenceRigidShakeAlgorithm::getMaximumNumberOfIterations( void ) const {
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceRigidShakeAlgorithm::setMaximumNumberOfIterations( int maximumNumberOfIterations ){
+int ReferenceCCMAAlgorithm::setMaximumNumberOfIterations( int maximumNumberOfIterations ){
 
    // ---------------------------------------------------------------------------------------
 
-   // static const char* methodName = "\nReferenceQShakeAlgorithm::setMaximumNumberOfIterations";
+   // static const char* methodName = "\nReferenceCCMAAlgorithm::setMaximumNumberOfIterations";
 
    // ---------------------------------------------------------------------------------------
 
@@ -459,11 +303,11 @@ int ReferenceRigidShakeAlgorithm::setMaximumNumberOfIterations( int maximumNumbe
 
    --------------------------------------------------------------------------------------- */
 
-RealOpenMM ReferenceRigidShakeAlgorithm::getTolerance( void ) const {
+RealOpenMM ReferenceCCMAAlgorithm::getTolerance( void ) const {
 
    // ---------------------------------------------------------------------------------------
 
-   // static const char* methodName = "\nReferenceQShakeAlgorithm::getTolerance";
+   // static const char* methodName = "\nReferenceCCMAAlgorithm::getTolerance";
 
    // ---------------------------------------------------------------------------------------
 
@@ -480,12 +324,12 @@ RealOpenMM ReferenceRigidShakeAlgorithm::getTolerance( void ) const {
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceRigidShakeAlgorithm::setTolerance( RealOpenMM tolerance ){
+int ReferenceCCMAAlgorithm::setTolerance( RealOpenMM tolerance ){
 
 
    // ---------------------------------------------------------------------------------------
 
-   // static const char* methodName = "\nReferenceQShakeAlgorithm::setTolerance";
+   // static const char* methodName = "\nReferenceCCMAAlgorithm::setTolerance";
 
    // ---------------------------------------------------------------------------------------
 
@@ -497,7 +341,7 @@ int ReferenceRigidShakeAlgorithm::setTolerance( RealOpenMM tolerance ){
 
 /**---------------------------------------------------------------------------------------
 
-   Apply Shake algorithm
+   Apply CCMA algorithm
 
    @param numberOfAtoms    number of atoms
    @param atomCoordinates  atom coordinates
@@ -509,13 +353,13 @@ int ReferenceRigidShakeAlgorithm::setTolerance( RealOpenMM tolerance ){
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceRigidShakeAlgorithm::apply( int numberOfAtoms, RealOpenMM** atomCoordinates,
+int ReferenceCCMAAlgorithm::apply( int numberOfAtoms, RealOpenMM** atomCoordinates,
                                          RealOpenMM** atomCoordinatesP,
                                          RealOpenMM* inverseMasses ){
 
    // ---------------------------------------------------------------------------------------
 
-   static const char* methodName = "\nReferenceQShakeAlgorithm::apply";
+   static const char* methodName = "\nReferenceCCMAAlgorithm::apply";
 
    static const RealOpenMM zero        =  0.0;
    static const RealOpenMM one         =  1.0;
@@ -570,8 +414,8 @@ int ReferenceRigidShakeAlgorithm::apply( int numberOfAtoms, RealOpenMM** atomCoo
 
    int iterations           = 0;
    int numberConverged      = 0;
-   vector<RealOpenMM> constraintForce(_numberOfConstraints);
-   vector<RealOpenMM> tempForce(_numberOfConstraints);
+   vector<RealOpenMM> constraintDelta(_numberOfConstraints);
+   vector<RealOpenMM> tempDelta(_numberOfConstraints);
    while (iterations < getMaximumNumberOfIterations()) {
       numberConverged  = 0;
       for( int ii = 0; ii < _numberOfConstraints; ii++ ){
@@ -586,14 +430,14 @@ int ReferenceRigidShakeAlgorithm::apply( int numberOfAtoms, RealOpenMM** atomCoo
          RealOpenMM rp2  = DOT3( rp_ij, rp_ij );
          RealOpenMM dist2= _distance[ii]*_distance[ii];
          RealOpenMM diff = dist2 - rp2;
-         constraintForce[ii] = zero;
+         constraintDelta[ii] = zero;
          RealOpenMM rrpr  = DOT3(  rp_ij, r_ij[ii] );
          if( rrpr <  d_ij2[ii]*epsilon6 ){
              std::stringstream message;
              message << iterations <<" "<<atomI<<" "<<atomJ<< " Error: sign of rrpr < 0?\n";
              SimTKOpenMMLog::printMessage( message );
          } else {
-             constraintForce[ii] = reducedMasses[ii]*diff/rrpr;
+             constraintDelta[ii] = reducedMasses[ii]*diff/rrpr;
          }
          if (rp2 >= lowerTol*dist2 && rp2 <= upperTol*dist2) {
             numberConverged++;
@@ -608,35 +452,18 @@ int ReferenceRigidShakeAlgorithm::apply( int numberOfAtoms, RealOpenMM** atomCoo
               RealOpenMM sum = 0.0;
               for (int j = 0; j < (int) _matrix[i].size(); j++) {
                   pair<int, RealOpenMM> element = _matrix[i][j];
-                  sum += element.second*constraintForce[element.first];
-              }
-              tempForce[i] = sum;
-          }
-          constraintForce = tempForce;
-      }
-      
-//      for (unsigned int i = 0; i < _rigidClusters.size(); i++) {
-//          const vector<int>& cluster = _rigidClusters[i];
-//          RealOpenMM** matrix = _matrices[i];
-//          unsigned int size = cluster.size();
-//          if (size > tempForce.size())
-//              tempForce.resize(size);
-//          for (unsigned int j = 0; j < size; j++) {
-//              tempForce[j] = zero;
-//              for (unsigned int k = 0; k < size; k++)
-//                  tempForce[j] += matrix[j][k]*constraintForce[cluster[k]];
-//          }
-//          for (unsigned int j = 0; j < size; j++)
-//              constraintForce[cluster[j]] = tempForce[j];
-//
-//      }
-RealOpenMM damping = one;//(RealOpenMM) (iterations < 2 ? 0.5 : 1.0);
+                  sum += element.second*constraintDelta[element.first];
+              }
+              tempDelta[i] = sum;
+          }
+          constraintDelta = tempDelta;
+      }
       for( int ii = 0; ii < _numberOfConstraints; ii++ ){
 
          int atomI   = _atomIndices[ii][0];
          int atomJ   = _atomIndices[ii][1];
          for( int jj = 0; jj < 3; jj++ ){
-            RealOpenMM dr                = constraintForce[ii]*r_ij[ii][jj]*damping;
+            RealOpenMM dr                = constraintDelta[ii]*r_ij[ii][jj];
             atomCoordinatesP[atomI][jj] += inverseMasses[atomI]*dr;
             atomCoordinatesP[atomJ][jj] -= inverseMasses[atomJ]*dr;
          }
@@ -655,7 +482,7 @@ RealOpenMM damping = one;//(RealOpenMM) (iterations < 2 ? 0.5 : 1.0);
          message << " FAILED";
       }
       message << "\n";
-      int errors = reportShake( numberOfAtoms, atomCoordinatesP, message );
+      int errors = reportCCMA( numberOfAtoms, atomCoordinatesP, message );
       if( !errors ){
          message << "*** no errors recorded in explicit check ***";
       }
@@ -679,12 +506,12 @@ RealOpenMM damping = one;//(RealOpenMM) (iterations < 2 ? 0.5 : 1.0);
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceRigidShakeAlgorithm::reportShake( int numberOfAtoms, RealOpenMM** atomCoordinates,
+int ReferenceCCMAAlgorithm::reportCCMA( int numberOfAtoms, RealOpenMM** atomCoordinates,
                                           std::stringstream& message ){
 
    // ---------------------------------------------------------------------------------------
 
-   static const char* methodName = "\nReferenceQShakeAlgorithm::reportShake";
+   static const char* methodName = "\nReferenceCCMAAlgorithm::reportCCMA";
 
    static const RealOpenMM zero        =  0.0;
    static const RealOpenMM one         =  1.0;

platforms/reference/src/SimTKReference/ReferenceRigidShakeAlgorithm.h → platforms/reference/src/SimTKReference/ReferenceCCMAAlgorithm.h

@@ -22,8 +22,8 @@
  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  */
 
-#ifndef __ReferenceQShakeAlgorithm_H__
-#define __ReferenceQShakeAlgorithm_H__
+#ifndef __ReferenceCCMAAlgorithm_H__
+#define __ReferenceCCMAAlgorithm_H__
 
 #include "ReferenceConstraintAlgorithm.h"
 #include <utility>
@@ -32,7 +32,7 @@
 
 // ---------------------------------------------------------------------------------------
 
-class ReferenceRigidShakeAlgorithm : public ReferenceConstraintAlgorithm {
+class ReferenceCCMAAlgorithm : public ReferenceConstraintAlgorithm {
 
    protected:
 
@@ -48,8 +48,6 @@ class ReferenceRigidShakeAlgorithm : public ReferenceConstraintAlgorithm {
       RealOpenMM* _distanceTolerance;
       RealOpenMM* _reducedMasses;
       bool _hasInitializedMasses;
-//      std::vector<std::vector<int> > _rigidClusters;
-//      std::vector<RealOpenMM**> _matrices;
       std::vector<std::vector<std::pair<int, RealOpenMM> > > _matrix;
 
    public:
@@ -57,17 +55,19 @@ class ReferenceRigidShakeAlgorithm : public ReferenceConstraintAlgorithm {
 
       /**---------------------------------------------------------------------------------------
 
-         ReferenceQShakeAlgorithm constructor
+         ReferenceCCMAAlgorithm constructor
 
          @param numberOfAtoms    number of atoms
          @param numberOfConstraints      number of constraints
          @param atomIndices              atom indices for contraints
          @param distance                 distances for constraints
+         @param masses                   atom masses
+         @param angles                   angle force field terms
          @param tolerance                constraint tolerance
 
          --------------------------------------------------------------------------------------- */
 
-      ReferenceRigidShakeAlgorithm( int numberOfAtoms, int numberOfConstraints, int** atomIndices, RealOpenMM* distance, RealOpenMM* masses, std::vector<AngleInfo>& angles, RealOpenMM tolerance );
+      ReferenceCCMAAlgorithm( int numberOfAtoms, int numberOfConstraints, int** atomIndices, RealOpenMM* distance, RealOpenMM* masses, std::vector<AngleInfo>& angles, RealOpenMM tolerance );
 
       /**---------------------------------------------------------------------------------------
 
@@ -75,7 +75,7 @@ class ReferenceRigidShakeAlgorithm : public ReferenceConstraintAlgorithm {
 
          --------------------------------------------------------------------------------------- */
 
-       ~ReferenceRigidShakeAlgorithm( );
+       ~ReferenceCCMAAlgorithm( );
 
       /**---------------------------------------------------------------------------------------
 
@@ -132,19 +132,7 @@ class ReferenceRigidShakeAlgorithm : public ReferenceConstraintAlgorithm {
 
       /**---------------------------------------------------------------------------------------
 
-         Print parameters
-
-         @param message message
-
-         @return ReferenceQShakeAlgorithm::DefaultReturn
-
-         --------------------------------------------------------------------------------------- */
-
-      int printParameters( std::stringstream& message ) const;
-
-      /**---------------------------------------------------------------------------------------
-
-         Apply Shake algorithm
+         Apply CCMA algorithm
 
          @param numberOfAtoms    number of atoms
          @param atomCoordinates  atom coordinates
@@ -171,10 +159,10 @@ class ReferenceRigidShakeAlgorithm : public ReferenceConstraintAlgorithm {
 
          --------------------------------------------------------------------------------------- */
 
-      int reportShake( int numberOfAtoms, RealOpenMM** atomCoordinates, std::stringstream& message );
+      int reportCCMA( int numberOfAtoms, RealOpenMM** atomCoordinates, std::stringstream& message );
 };
 
-class ReferenceRigidShakeAlgorithm::AngleInfo
+class ReferenceCCMAAlgorithm::AngleInfo
 {
 public:
     int atom1, atom2, atom3;
@@ -187,4 +175,4 @@ public:
 
 // ---------------------------------------------------------------------------------------
 
-#endif // __ReferenceQShakeAlgorithm_H__
+#endif // __ReferenceCCMAAlgorithm_H__