Preliminary implementation of LINCS (both reference and Cuda)

platforms/cuda/cuda-cmake/FindCuda.cmake

@@ -117,7 +117,7 @@ ELSE(CUDA_BUILD_TYPE MATCHES "Emulation")
 ENDIF(CUDA_BUILD_TYPE MATCHES "Emulation")
 
 SET(CUDA_BUILD_CUBIN TRUE CACHE BOOL "Generate and parse .cubin files in Device mode.")
-SET(CUDA_NVCC_FLAGS "-maxrregcount=32;-use_fast_math;-O0" CACHE STRING "Semi-colon delimit multiple arguments.")
+SET(CUDA_NVCC_FLAGS "-maxrregcount=32;-use_fast_math;-O0;-arch=sm_11" CACHE STRING "Semi-colon delimit multiple arguments.")
 
 # Search for the cuda distribution.
 IF(NOT CUDA_INSTALL_PREFIX)

platforms/cuda/src/CudaKernels.cpp

@@ -419,6 +419,7 @@ void CudaIntegrateVerletStepKernel::execute(OpenMMContextImpl& context, const Ve
     kVerletUpdatePart1(gpu);
     kApplyFirstShake(gpu);
     kApplyFirstSettle(gpu);
+    kApplyFirstLincs(gpu);
     if (data.removeCM) {
         int step = (int) (context.getTime()/stepSize);
         if (step%data.cmMotionFrequency == 0)
@@ -457,6 +458,7 @@ void CudaIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, const
     kUpdatePart1(gpu);
     kApplyFirstShake(gpu);
     kApplyFirstSettle(gpu);
+    kApplyFirstLincs(gpu);
     if (data.removeCM) {
         int step = (int) (context.getTime()/stepSize);
         if (step%data.cmMotionFrequency == 0)
@@ -465,6 +467,7 @@ void CudaIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, const
     kUpdatePart2(gpu);
     kApplySecondShake(gpu);
     kApplySecondSettle(gpu);
+    kApplySecondLincs(gpu);
 }
 
 CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
@@ -497,6 +500,7 @@ void CudaIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, const
     kBrownianUpdatePart1(gpu);
     kApplyFirstShake(gpu);
     kApplyFirstSettle(gpu);
+    kApplyFirstLincs(gpu);
     if (data.removeCM) {
         int step = (int) (context.getTime()/stepSize);
         if (step%data.cmMotionFrequency == 0)

platforms/cuda/src/kernels/cudaKernels.h

@@ -48,9 +48,11 @@ extern void kReduceBornSumAndForces(gpuContext gpu);
 extern void kUpdatePart1(gpuContext gpu);
 extern void kApplyFirstShake(gpuContext gpu);
 extern void kApplyFirstSettle(gpuContext gpu);
+extern void kApplyFirstLincs(gpuContext gpu);
 extern void kUpdatePart2(gpuContext gpu);
 extern void kApplySecondShake(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);
@@ -79,6 +81,8 @@ extern void SetUpdateShakeHSim(gpuContext gpu);
 extern void GetUpdateShakeHSim(gpuContext gpu);
 extern void SetSettleSim(gpuContext gpu);
 extern void GetSettleSim(gpuContext gpu);
+extern void SetLincsSim(gpuContext gpu);
+extern void GetLincsSim(gpuContext gpu);
 extern void SetVerletUpdateSim(gpuContext gpu);
 extern void GetVerletUpdateSim(gpuContext gpu);
 extern void SetBrownianUpdateSim(gpuContext gpu);

platforms/cuda/src/kernels/cudatypes.h

@@ -254,6 +254,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    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
@@ -345,6 +346,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    NonShakeConstraints;            // Total number of NonShake atoms
     unsigned int    maxShakeIterations;             // Maximum shake iterations
     unsigned int    degreesOfFreedom;               // Number of degrees of freedom in system
@@ -365,6 +367,18 @@ 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*            pLincsConnectionsIndex;         // The index in pLincsConnections at which the connections for a particular constraint start
+    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*            pLincsAtomConstraintsIndex;     // The index in psLincsAtomConstraints at which the constraints for a particular atom start
+    unsigned int*   pSyncCounter;                   // Used for global thread synchronization
 
     // Mutable stuff
     float4*         pPosq;                          // Pointer to atom positions and charges

platforms/cuda/src/kernels/gpu.cpp

@@ -546,6 +546,106 @@ void gpuSetShakeParameters(gpuContext gpu, const vector<int>& atom1, const vecto
     if (gpu->sim.settle_threads_per_block < 1)
         gpu->sim.settle_threads_per_block = 1;
 
+    // Find connected constraints for LINCS.
+
+    vector<int> lincsConstraints;
+    for (unsigned int i = 0; i < atom1.size(); i++)
+        if (settleConstraints[atom1[i]].size() != 2)
+            lincsConstraints.push_back(i);
+    vector<vector<int> > atomConstraints(gpu->natoms);
+    for (unsigned int i = 0; i < lincsConstraints.size(); i++) {
+        atomConstraints[atom1[lincsConstraints[i]]].push_back(i);
+        atomConstraints[atom2[lincsConstraints[i]]].push_back(i);
+    }
+    vector<vector<int> > linkedConstraints(lincsConstraints.size());
+    for (unsigned int atom = 0; atom < atomConstraints.size(); atom++) {
+        for (unsigned int i = 0; i < atomConstraints[atom].size(); i++)
+            for (int j = 0; j < i; j++) {
+                int c1 = atomConstraints[atom][i];
+                int c2 = atomConstraints[atom][j];
+                linkedConstraints[c1].push_back(c2);
+                linkedConstraints[c2].push_back(c1);
+            }
+    }
+    int totalLinks = 0;
+    for (unsigned int i = 0; i < linkedConstraints.size(); i++)
+        totalLinks += linkedConstraints[i].size();
+
+    // Fill in the CUDA streams.
+
+    CUDAStream<int2>* psLincsAtoms = new CUDAStream<int2>((int) lincsConstraints.size(), 1);
+    gpu->psLincsAtoms              = psLincsAtoms;
+    gpu->sim.pLincsAtoms           = psLincsAtoms->_pDevData;
+    CUDAStream<float4>* psLincsDistance = new CUDAStream<float4>((int) lincsConstraints.size(), 1);
+    gpu->psLincsDistance                = psLincsDistance;
+    gpu->sim.pLincsDistance             = psLincsDistance->_pDevData;
+    CUDAStream<int>* psLincsConnections = new CUDAStream<int>(totalLinks, 1);
+    gpu->psLincsConnections             = psLincsConnections;
+    gpu->sim.pLincsConnections          = psLincsConnections->_pDevData;
+    CUDAStream<int>* psLincsConnectionsIndex = new CUDAStream<int>((int) lincsConstraints.size()+1, 1);
+    gpu->psLincsConnectionsIndex             = psLincsConnectionsIndex;
+    gpu->sim.pLincsConnectionsIndex          = psLincsConnectionsIndex->_pDevData;
+    CUDAStream<int>* psLincsAtomConstraints = new CUDAStream<int>((int) lincsConstraints.size()*2, 1);
+    gpu->psLincsAtomConstraints             = psLincsAtomConstraints;
+    gpu->sim.pLincsAtomConstraints          = psLincsAtomConstraints->_pDevData;
+    CUDAStream<int>* psLincsAtomConstraintsIndex = new CUDAStream<int>(gpu->natoms+1, 1);
+    gpu->psLincsAtomConstraintsIndex             = psLincsAtomConstraintsIndex;
+    gpu->sim.pLincsAtomConstraintsIndex          = psLincsAtomConstraintsIndex->_pDevData;
+    CUDAStream<float>* psLincsS = new CUDAStream<float>((int) lincsConstraints.size(), 1);
+    gpu->psLincsS             = psLincsS;
+    gpu->sim.pLincsS          = psLincsS->_pDevData;
+    CUDAStream<float>* psLincsCoupling = new CUDAStream<float>(totalLinks, 1);
+    gpu->psLincsCoupling               = psLincsCoupling;
+    gpu->sim.pLincsCoupling            = psLincsCoupling->_pDevData;
+    CUDAStream<float>* psLincsRhs1 = new CUDAStream<float>((int) lincsConstraints.size(), 1);
+    gpu->psLincsRhs1             = psLincsRhs1;
+    gpu->sim.pLincsRhs1          = psLincsRhs1->_pDevData;
+    CUDAStream<float>* psLincsRhs2 = new CUDAStream<float>((int) lincsConstraints.size(), 1);
+    gpu->psLincsRhs2             = psLincsRhs2;
+    gpu->sim.pLincsRhs2          = psLincsRhs2->_pDevData;
+    CUDAStream<float>* psLincsSolution = new CUDAStream<float>((int) lincsConstraints.size(), 1);
+    gpu->psLincsSolution             = psLincsSolution;
+    gpu->sim.pLincsSolution          = psLincsSolution->_pDevData;
+    CUDAStream<unsigned int>* psSyncCounter = new CUDAStream<unsigned int>(100, 1);
+    gpu->psSyncCounter                      = psSyncCounter;
+    gpu->sim.pSyncCounter                   = psSyncCounter->_pDevData;
+    gpu->sim.lincsConstraints = lincsConstraints.size();
+    int index = 0;
+    for (unsigned int i = 0; i < lincsConstraints.size(); i++) {
+        int c = lincsConstraints[i];
+        psLincsAtoms->_pSysData[i].x = atom1[c];
+        psLincsAtoms->_pSysData[i].y = atom2[c];
+        psLincsDistance->_pSysData[i].w = distance[c];
+        psLincsS->_pSysData[i] = 1.0f/sqrt(invMass1[c]+invMass2[c]);
+        psLincsConnectionsIndex->_pSysData[i] = index;
+        for (unsigned int j = 0; j < linkedConstraints[i].size(); j++)
+            psLincsConnections->_pSysData[index++] = linkedConstraints[i][j];
+    }
+    psLincsConnectionsIndex->_pSysData[lincsConstraints.size()] = index;
+    for (unsigned int i = 0; i < psSyncCounter->_length; i++)
+        psSyncCounter->_pSysData[i] = 0;
+    index = 0;
+    for (unsigned int i = 0; i < atomConstraints.size(); i++) {
+        psLincsAtomConstraintsIndex->_pSysData[i] = index;
+        for (unsigned int j = 0; j < atomConstraints[i].size(); j++)
+            psLincsAtomConstraints->_pSysData[index++] = atomConstraints[i][j];
+    }
+    psLincsAtomConstraintsIndex->_pSysData[atomConstraints.size()] = index;
+    psLincsAtoms->Upload();
+    psLincsDistance->Upload();
+    psLincsS->Upload();
+    psLincsConnections->Upload();
+    psLincsConnectionsIndex->Upload();
+    psLincsAtomConstraints->Upload();
+    psLincsAtomConstraintsIndex->Upload();
+    psSyncCounter->Upload();
+    gpu->sim.lincs_threads_per_block = (gpu->sim.lincsConstraints + gpu->sim.blocks - 1) / gpu->sim.blocks;
+    if (gpu->sim.lincs_threads_per_block > gpu->sim.max_shake_threads_per_block)
+        gpu->sim.lincs_threads_per_block = gpu->sim.max_shake_threads_per_block;
+    if (gpu->sim.lincs_threads_per_block < 1)
+        gpu->sim.lincs_threads_per_block = 1;
+
+/*
     // Find clusters consisting of a central atom with up to three peripheral atoms.
     
     map<int, ShakeCluster> clusters;
@@ -614,6 +714,11 @@ void gpuSetShakeParameters(gpuContext gpu, const vector<int>& atom1, const vecto
     psShakeID->Upload();
     psShakeParameter->Upload();
     gpu->sim.shakeTolerance = tolerance;
+*/
+
+
+    gpu->sim.ShakeConstraints = 0;
+
 
     gpu->sim.shake_threads_per_block     = (gpu->sim.ShakeConstraints + gpu->sim.blocks - 1) / gpu->sim.blocks; 
     if (gpu->sim.shake_threads_per_block > gpu->sim.max_shake_threads_per_block)
@@ -627,7 +732,7 @@ void gpuSetShakeParameters(gpuContext gpu, const vector<int>& atom1, const vecto
 
     int count = 0;
     for (int i = 0; i < gpu->natoms; i++)
-       if (constraintCount[i] == 0)
+//       if (constraintCount[i] == 0)
           count++;
 
     // Allocate NonShake parameters
@@ -651,9 +756,9 @@ void gpuSetShakeParameters(gpuContext gpu, const vector<int>& atom1, const vecto
 
        count = 0;
        for (int i = 0; i < gpu->natoms; i++){
-          if (constraintCount[i] == 0){
+//          if (constraintCount[i] == 0){
              psNonShakeID->_pSysStream[0][count++] = i;
-          }
+//          }
        }
        psNonShakeID->Upload();
 
@@ -978,7 +1083,18 @@ void* gpuInit(int numAtoms)
     gpu->psInteractionCount         = NULL;
     gpu->psGridBoundingBox          = NULL;
     gpu->psGridCenter               = NULL;
-
+    gpu->psLincsAtoms               = NULL;
+    gpu->psLincsDistance            = NULL;
+    gpu->psLincsConnections         = NULL;
+    gpu->psLincsConnectionsIndex    = NULL;
+    gpu->psLincsAtomConstraints     = NULL;
+    gpu->psLincsAtomConstraintsIndex= NULL;
+    gpu->psLincsS                   = NULL;
+    gpu->psLincsCoupling            = NULL;
+    gpu->psLincsRhs1                = NULL;
+    gpu->psLincsRhs2                = NULL;
+    gpu->psLincsSolution            = NULL;
+    gpu->psSyncCounter              = NULL;
 
     // Initialize output buffer before reading parameters
     gpu->pOutputBufferCounter       = new unsigned int[gpu->sim.paddedNumberOfAtoms];
@@ -1117,6 +1233,18 @@ void gpuShutDown(gpuContext gpu)
     delete gpu->psAtomIndex;
     delete gpu->psGridBoundingBox;
     delete gpu->psGridCenter;
+    delete gpu->psLincsAtoms;
+    delete gpu->psLincsDistance;
+    delete gpu->psLincsConnections;
+    delete gpu->psLincsConnectionsIndex;
+    delete gpu->psLincsAtomConstraints;
+    delete gpu->psLincsAtomConstraintsIndex;
+    delete gpu->psLincsS;
+    delete gpu->psLincsCoupling;
+    delete gpu->psLincsRhs1;
+    delete gpu->psLincsRhs2;
+    delete gpu->psLincsSolution;
+    delete gpu->psSyncCounter;
     if (gpu->cudpp != 0)
         cudppDestroyPlan(gpu->cudpp);
 
@@ -1424,6 +1552,7 @@ int gpuSetConstants(gpuContext gpu)
     SetVerletUpdateSim(gpu);
     SetBrownianUpdateSim(gpu);
     SetSettleSim(gpu);
+    SetLincsSim(gpu);
     SetRandomSim(gpu);
     return 1;
 }

platforms/cuda/src/kernels/gputypes.h

@@ -129,6 +129,18 @@ 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>* psLincsConnectionsIndex; // The index in psLincsConnections at which the connections for a particular constraint start
+    CUDAStream<int>* psLincsAtomConstraints; // The indices of constraints involving each atom
+    CUDAStream<int>* psLincsAtomConstraintsIndex; // The index in psLincsAtomConstraints at which the constraints for a particular atom start
+    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<unsigned int>* psSyncCounter;// Used for global thread synchronization
 };
 
 typedef struct _gpuContext *gpuContext;

platforms/cuda/tests/TestCudaLangevinIntegrator.cpp

@@ -131,7 +131,7 @@ void testTemperature() {
 
 void testConstraints() {
     const int numParticles = 8;
-    const int numConstraints = 4;
+    const int numConstraints = 7;
     const double temp = 100.0;
     CudaPlatform platform;
     System system(numParticles, numConstraints);
@@ -143,7 +143,7 @@ void testConstraints() {
         forceField->setParticleParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
     }
     for (int i = 0; i < numConstraints; ++i)
-        system.setConstraintParameters(i, 2*i, 2*i+1, 1.0);
+        system.setConstraintParameters(i, i, i+1, 1.0);
     system.addForce(forceField);
     OpenMMContext context(system, integrator, platform);
     vector<Vec3> positions(numParticles);
@@ -241,7 +241,7 @@ int main() {
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;
-        return 1;
+//        return 1;
     }
     cout << "Done" << endl;
     return 0;

platforms/reference/src/ReferenceKernels.cpp

@@ -37,6 +37,7 @@
 #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"
@@ -461,14 +462,14 @@ double ReferenceCalcGBSAOBCForceKernel::executeEnergy(OpenMMContextImpl& context
 ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
     if (dynamics)
         delete dynamics;
-    if (shake)
-        delete shake;
+    if (constraints)
+        delete constraints;
     if (masses)
         delete[] masses;
     if (constraintIndices)
         disposeIntArray(constraintIndices, numConstraints);
-    if (shakeParameters)
-        disposeRealArray(shakeParameters, numConstraints);
+    if (constraintDistances)
+        delete[] constraintDistances;
 }
 
 void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
@@ -478,14 +479,14 @@ void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
     numConstraints = system.getNumConstraints();
     constraintIndices = allocateIntArray(numConstraints, 2);
-    shakeParameters = allocateRealArray(numConstraints, 1);
+    constraintDistances = new RealOpenMM[numConstraints];
     for (int i = 0; i < numConstraints; ++i) {
         int particle1, particle2;
         double distance;
         system.getConstraintParameters(i, particle1, particle2, distance);
         constraintIndices[i][0] = particle1;
         constraintIndices[i][1] = particle2;
-        shakeParameters[i][0] = static_cast<RealOpenMM>(distance);
+        constraintDistances[i] = static_cast<RealOpenMM>(distance);
     }
 }
 
@@ -499,28 +500,28 @@ void ReferenceIntegrateVerletStepKernel::execute(OpenMMContextImpl& context, con
         
         if (dynamics) {
             delete dynamics;
-            delete shake;
+            delete constraints;
         }
         dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) );
-        shake = new ReferenceShakeAlgorithm(numConstraints, constraintIndices, shakeParameters);
-        dynamics->setReferenceShakeAlgorithm(shake);
+        constraints = new ReferenceShakeAlgorithm(numConstraints, constraintIndices, constraintDistances);
+        dynamics->setReferenceConstraintAlgorithm(constraints);
         prevStepSize = stepSize;
     }
-    shake->setTolerance(integrator.getConstraintTolerance());
+//    shake->setTolerance(integrator.getConstraintTolerance());
     dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
 }
 
 ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
     if (dynamics)
         delete dynamics;
-    if (shake)
-        delete shake;
+    if (constraints)
+        delete constraints;
     if (masses)
         delete[] masses;
     if (constraintIndices)
         disposeIntArray(constraintIndices, numConstraints);
-    if (shakeParameters)
-        disposeRealArray(shakeParameters, numConstraints);
+    if (constraintDistances)
+        delete[] constraintDistances;
 }
 
 void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
@@ -530,14 +531,14 @@ void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, cons
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
     numConstraints = system.getNumConstraints();
     constraintIndices = allocateIntArray(numConstraints, 2);
-    shakeParameters = allocateRealArray(numConstraints, 1);
+    constraintDistances = new RealOpenMM[numConstraints];
     for (int i = 0; i < numConstraints; ++i) {
         int particle1, particle2;
         double distance;
         system.getConstraintParameters(i, particle1, particle2, distance);
         constraintIndices[i][0] = particle1;
         constraintIndices[i][1] = particle2;
-        shakeParameters[i][0] = static_cast<RealOpenMM>(distance);
+        constraintDistances[i] = static_cast<RealOpenMM>(distance);
     }
     SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
 }
@@ -554,7 +555,7 @@ void ReferenceIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, c
         
         if (dynamics) {
             delete dynamics;
-            delete shake;
+            delete constraints;
         }
         RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
         dynamics = new ReferenceStochasticDynamics(
@@ -562,27 +563,28 @@ void ReferenceIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, c
 				static_cast<RealOpenMM>(stepSize), 
 				static_cast<RealOpenMM>(tau), 
 				static_cast<RealOpenMM>(temperature) );
-        shake = new ReferenceShakeAlgorithm(numConstraints, constraintIndices, shakeParameters);
-        dynamics->setReferenceShakeAlgorithm(shake);
+        constraints = new ReferenceLincsAlgorithm(numConstraints, constraintIndices, constraintDistances);
+//        shake = new ReferenceShakeAlgorithm(numConstraints, constraintIndices, shakeParameters);
+        dynamics->setReferenceConstraintAlgorithm(constraints);
         prevTemp = temperature;
         prevFriction = friction;
         prevStepSize = stepSize;
     }
-    shake->setTolerance(integrator.getConstraintTolerance());
+//    shake->setTolerance(integrator.getConstraintTolerance());
     dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
 }
 
 ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
     if (dynamics)
         delete dynamics;
-    if (shake)
-        delete shake;
+    if (constraints)
+        delete constraints;
     if (masses)
         delete[] masses;
     if (constraintIndices)
         disposeIntArray(constraintIndices, numConstraints);
-    if (shakeParameters)
-        disposeRealArray(shakeParameters, numConstraints);
+    if (constraintDistances)
+        delete[] constraintDistances;
 }
 
 void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
@@ -592,14 +594,14 @@ void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, cons
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
     numConstraints = system.getNumConstraints();
     constraintIndices = allocateIntArray(numConstraints, 2);
-    shakeParameters = allocateRealArray(numConstraints, 1);
+    constraintDistances = new RealOpenMM[numConstraints];
     for (int i = 0; i < numConstraints; ++i) {
         int particle1, particle2;
         double distance;
         system.getConstraintParameters(i, particle1, particle2, distance);
         constraintIndices[i][0] = particle1;
         constraintIndices[i][1] = particle2;
-        shakeParameters[i][0] = static_cast<RealOpenMM>(distance);
+        constraintDistances[i] = static_cast<RealOpenMM>(distance);
     }
     SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
 }
@@ -616,20 +618,20 @@ void ReferenceIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, c
         
         if (dynamics) {
             delete dynamics;
-            delete shake;
+            delete constraints;
         }
         dynamics = new ReferenceBrownianDynamics(
 				context.getSystem().getNumParticles(), 
 				static_cast<RealOpenMM>(stepSize), 
 				static_cast<RealOpenMM>(friction), 
 				static_cast<RealOpenMM>(temperature) );
-        shake = new ReferenceShakeAlgorithm(numConstraints, constraintIndices, shakeParameters);
-        dynamics->setReferenceShakeAlgorithm(shake);
+        constraints = new ReferenceShakeAlgorithm(numConstraints, constraintIndices, constraintDistances);
+        dynamics->setReferenceConstraintAlgorithm(constraints);
         prevTemp = temperature;
         prevFriction = friction;
         prevStepSize = stepSize;
     }
-    shake->setTolerance(integrator.getConstraintTolerance());
+//    shake->setTolerance(integrator.getConstraintTolerance());
     dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
 }
 

platforms/reference/src/ReferenceKernels.h

@@ -40,7 +40,7 @@ class CpuObc;
 class ReferenceAndersenThermostat;
 class ReferenceBrownianDynamics;
 class ReferenceStochasticDynamics;
-class ReferenceShakeAlgorithm;
+class ReferenceConstraintAlgorithm;
 class ReferenceVerletDynamics;
 
 namespace OpenMM {
@@ -282,7 +282,7 @@ private:
 class ReferenceIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
 public:
     ReferenceIntegrateVerletStepKernel(std::string name, const Platform& platform) : IntegrateVerletStepKernel(name, platform),
-        dynamics(0), shake(0), masses(0), shakeParameters(0), constraintIndices(0) {
+        dynamics(0), constraints(0), masses(0), constraintDistances(0), constraintIndices(0) {
     }
     ~ReferenceIntegrateVerletStepKernel();
     /**
@@ -301,9 +301,9 @@ public:
     void execute(OpenMMContextImpl& context, const VerletIntegrator& integrator);
 private:
     ReferenceVerletDynamics* dynamics;
-    ReferenceShakeAlgorithm* shake;
+    ReferenceConstraintAlgorithm* constraints;
     RealOpenMM* masses;
-    RealOpenMM** shakeParameters;
+    RealOpenMM* constraintDistances;
     int** constraintIndices;
     int numConstraints;
     double prevStepSize;
@@ -315,7 +315,7 @@ private:
 class ReferenceIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
 public:
     ReferenceIntegrateLangevinStepKernel(std::string name, const Platform& platform) : IntegrateLangevinStepKernel(name, platform),
-        dynamics(0), shake(0), masses(0), shakeParameters(0), constraintIndices(0) {
+        dynamics(0), constraints(0), masses(0), constraintDistances(0), constraintIndices(0) {
     }
     ~ReferenceIntegrateLangevinStepKernel();
     /**
@@ -334,9 +334,9 @@ public:
     void execute(OpenMMContextImpl& context, const LangevinIntegrator& integrator);
 private:
     ReferenceStochasticDynamics* dynamics;
-    ReferenceShakeAlgorithm* shake;
+    ReferenceConstraintAlgorithm* constraints;
     RealOpenMM* masses;
-    RealOpenMM** shakeParameters;
+    RealOpenMM* constraintDistances;
     int** constraintIndices;
     int numConstraints;
     double prevTemp, prevFriction, prevStepSize;
@@ -348,7 +348,7 @@ private:
 class ReferenceIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel {
 public:
     ReferenceIntegrateBrownianStepKernel(std::string name, const Platform& platform) : IntegrateBrownianStepKernel(name, platform),
-        dynamics(0), shake(0), masses(0), shakeParameters(0), constraintIndices(0) {
+        dynamics(0), constraints(0), masses(0), constraintDistances(0), constraintIndices(0) {
     }
     ~ReferenceIntegrateBrownianStepKernel();
     /**
@@ -367,9 +367,9 @@ public:
     void execute(OpenMMContextImpl& context, const BrownianIntegrator& integrator);
 private:
     ReferenceBrownianDynamics* dynamics;
-    ReferenceShakeAlgorithm* shake;
+    ReferenceConstraintAlgorithm* constraints;
     RealOpenMM* masses;
-    RealOpenMM** shakeParameters;
+    RealOpenMM* constraintDistances;
     int** constraintIndices;
     int numConstraints;
     double prevTemp, prevFriction, prevStepSize;

platforms/reference/src/SimTKReference/ReferenceBrownianDynamics.cpp

@@ -203,9 +203,9 @@ int ReferenceBrownianDynamics::update( int numberOfAtoms, RealOpenMM** atomCoord
            xPrime[i][j] = atomCoordinates[i][j] + forceScale*forces[i][j] + noiseAmplitude*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
        }
    }
-   ReferenceShakeAlgorithm* referenceShakeAlgorithm = getReferenceShakeAlgorithm();
-   if( referenceShakeAlgorithm )
-      referenceShakeAlgorithm->applyShake( numberOfAtoms, atomCoordinates, xPrime, inverseMasses );
+   ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
+   if( referenceConstraintAlgorithm )
+      referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime, inverseMasses );
    
    // Update the positions and velocities.
    

platforms/reference/src/SimTKReference/ReferenceConstraintAlgorithm.h

+
+/* Portions copyright (c) 2009 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * 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.
+ */
+
+#ifndef __ReferenceConstraintAlgorithm_H__
+#define __ReferenceConstraintAlgorithm_H__
+
+/**
+ * This abstract class defines the interface which constraint algorithms must implement.
+ */
+class ReferenceConstraintAlgorithm {
+public:
+
+    virtual ~ReferenceConstraintAlgorithm() {};
+
+      /**---------------------------------------------------------------------------------------
+
+         Apply constraint algorithm
+
+         @param numberOfAtoms    number of atoms
+         @param atomCoordinates  atom coordinates
+         @param atomCoordinatesP atom coordinates prime
+         @param inverseMasses    1/mass
+
+         @return ReferenceDynamics::DefaultReturn if converge; else
+          return ReferenceDynamics::ErrorReturn
+
+         --------------------------------------------------------------------------------------- */
+
+    virtual int apply(int numberOfAtoms, RealOpenMM** atomCoordinates,
+                     RealOpenMM** atomCoordinatesP, RealOpenMM* inverseMasses) = 0;
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // __ReferenceConstraintAlgorithm_H__

platforms/reference/src/SimTKReference/ReferenceDynamics.cpp

@@ -63,8 +63,8 @@ ReferenceDynamics::ReferenceDynamics( int numberOfAtoms,  RealOpenMM deltaT, Rea
    _oneDTempArrays       = 0;
    _oneDTempArrays       = NULL;
 
-   _ownReferenceShake    = false;
-   _referenceShake       = NULL;
+   _ownReferenceConstraint = false;
+   _referenceConstraint    = NULL;
 }
 
 /**---------------------------------------------------------------------------------------
@@ -84,8 +84,8 @@ ReferenceDynamics::~ReferenceDynamics( ){
    _freeTwoDArrays();
    _freeOneDArrays();
 
-   if( _ownReferenceShake ){
-      delete _referenceShake;
+   if( _ownReferenceConstraint ){
+      delete _referenceConstraint;
    }
 }
 
@@ -383,49 +383,49 @@ RealOpenMM ReferenceDynamics::getTemperature( void ) const {
 
 /**---------------------------------------------------------------------------------------
 
-   Get ReferenceShake
+   Get ReferenceConstraint
 
-   @return referenceShake  object
+   @return ReferenceConstraint  object
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceShakeAlgorithm* ReferenceDynamics::getReferenceShakeAlgorithm( void ) const {
+ReferenceConstraintAlgorithm* ReferenceDynamics::getReferenceConstraintAlgorithm( void ) const {
 
    // ---------------------------------------------------------------------------------------
 
-   // static const char* methodName  = "\nReferenceDynamics::getReferenceShake";
+   // static const char* methodName  = "\nReferenceDynamics::getReferenceConstraint";
 
    // ---------------------------------------------------------------------------------------
 
-   return _referenceShake;
+   return _referenceConstraint;
 }
 
 /**---------------------------------------------------------------------------------------
 
-   Set ReferenceShake
+   Set ReferenceConstraint
 
-   @param referenceShake  referenceShake object
+   @param referenceConstraint  ReferenceConstraint object
 
    @return ReferenceDynamics::DefaultReturn
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceDynamics::setReferenceShakeAlgorithm( ReferenceShakeAlgorithm* referenceShake ){
+int ReferenceDynamics::setReferenceConstraintAlgorithm( ReferenceConstraintAlgorithm* referenceConstraint ){
 
    // ---------------------------------------------------------------------------------------
 
-   // static const char* methodName  = "\nReferenceDynamics::setReferenceShake";
+   // static const char* methodName  = "\nReferenceDynamics::setReferenceConstraint";
 
    // ---------------------------------------------------------------------------------------
 
    // delete if own
 
-   if( _referenceShake && _ownReferenceShake ){
-      delete _referenceShake;
+   if( _referenceConstraint && _ownReferenceConstraint ){
+      delete _referenceConstraint;
    }
 
-   _referenceShake     = referenceShake;
-   _ownReferenceShake  = 0;
+   _referenceConstraint = referenceConstraint;
+   _ownReferenceConstraint = 0;
 
    return ReferenceDynamics::DefaultReturn;
 }

platforms/reference/src/SimTKReference/ReferenceDynamics.h

@@ -25,7 +25,7 @@
 #ifndef __ReferenceDynamics_H__
 #define __ReferenceDynamics_H__
 
-#include "ReferenceShakeAlgorithm.h"
+#include "ReferenceConstraintAlgorithm.h"
 #include "../SimTKUtilities/SimTKOpenMMCommon.h"
 #include <cstddef>
 
@@ -68,8 +68,8 @@ class ReferenceDynamics {
       int _numberOf1DTempArrays;
       RealOpenMM** _oneDTempArrays;
 
-      int _ownReferenceShake;
-      ReferenceShakeAlgorithm* _referenceShake;
+      int _ownReferenceConstraint;
+      ReferenceConstraintAlgorithm* _referenceConstraint;
 
       /**---------------------------------------------------------------------------------------
       
@@ -261,25 +261,25 @@ class ReferenceDynamics {
 
       /**---------------------------------------------------------------------------------------
       
-         Get ReferenceShake
+         Get ReferenceConstraint
       
-         @return referenceShake  object
+         @return referenceConstraint  object
       
          --------------------------------------------------------------------------------------- */
       
-      ReferenceShakeAlgorithm* getReferenceShakeAlgorithm( void ) const;
+      ReferenceConstraintAlgorithm* getReferenceConstraintAlgorithm( void ) const;
       
       /**---------------------------------------------------------------------------------------
       
-         Set ReferenceShake
+         Set ReferenceConstraint
       
-         @param referenceShake  referenceShake object
+         @param referenceConstraint  referenceConstraint object
       
          @return ReferenceDynamics::DefaultReturn
       
          --------------------------------------------------------------------------------------- */
       
-      int setReferenceShakeAlgorithm( ReferenceShakeAlgorithm* referenceShake );
+      int setReferenceConstraintAlgorithm( ReferenceConstraintAlgorithm* referenceConstraint );
       
       /**---------------------------------------------------------------------------------------
       

platforms/reference/src/SimTKReference/ReferenceLincsAlgorithm.cpp

+/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * 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 <string.h>
+#include <sstream>
+
+#include "../SimTKUtilities/SimTKOpenMMCommon.h"
+#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
+#include "../SimTKUtilities/SimTKOpenMMLog.h"
+#include "ReferenceLincsAlgorithm.h"
+#include "ReferenceDynamics.h"
+
+using std::vector;
+using OpenMM::Vec3;
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceLincsAlgorithm constructor
+
+   @param numberOfConstraints      number of constraints
+   @param atomIndices              block of atom indices
+   @param distance                 distances for constraints
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceLincsAlgorithm::ReferenceLincsAlgorithm( int numberOfConstraints,
+                                                  int** atomIndices,
+                                                  RealOpenMM* distance ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceLincsAlgorithm::ReferenceLincsAlgorithm";
+
+   // ---------------------------------------------------------------------------------------
+
+   _numberOfConstraints        = numberOfConstraints;
+   _atomIndices                = atomIndices;
+   _distance                   = distance;
+
+   _numTerms                   = 8;
+   _hasInitialized             = false;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get number of constraints
+
+   @return number of constraints
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceLincsAlgorithm::getNumberOfConstraints( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceLincsAlgorithm::getNumberOfConstraints";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _numberOfConstraints;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get the number of terms to use in the series expansion
+
+   @return terms
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceLincsAlgorithm::getNumTerms( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceLincsAlgorithm::getNumTerms";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _numTerms;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set the number of terms to use in the series expansion
+
+   --------------------------------------------------------------------------------------- */
+
+void ReferenceLincsAlgorithm::setNumTerms( int terms ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceLincsAlgorithm::setNumTerms";
+
+   // ---------------------------------------------------------------------------------------
+
+   _numTerms = terms;
+}
+
+
+/**---------------------------------------------------------------------------------------
+
+   Initialize internal data structures.
+
+   @param numberOfAtoms    number of atoms
+   @param inverseMasses    1/mass
+
+   --------------------------------------------------------------------------------------- */
+
+void ReferenceLincsAlgorithm::initialize(int numberOfAtoms, RealOpenMM* inverseMasses) {
+    static const RealOpenMM one = 1.0;
+    _hasInitialized = true;
+    vector<vector<int> > atomConstraints(numberOfAtoms);
+    for (int constraint = 0; constraint < _numberOfConstraints; constraint++) {
+        atomConstraints[_atomIndices[constraint][0]].push_back(constraint);
+        atomConstraints[_atomIndices[constraint][1]].push_back(constraint);
+    }
+    _linkedConstraints.resize(_numberOfConstraints);
+    for (int atom = 0; atom < numberOfAtoms; atom++) {
+        for (int i = 0; i < atomConstraints[atom].size(); i++)
+            for (int j = 0; j < i; j++) {
+                int c1 = atomConstraints[atom][i];
+                int c2 = atomConstraints[atom][j];
+                _linkedConstraints[c1].push_back(c2);
+                _linkedConstraints[c2].push_back(c1);
+            }
+    }
+    _sMatrix.resize(_numberOfConstraints);
+    for (int constraint = 0; constraint < _numberOfConstraints; constraint++)
+        _sMatrix[constraint] = one/SQRT(inverseMasses[_atomIndices[constraint][0]]+inverseMasses[_atomIndices[constraint][1]]);
+    _couplingMatrix.resize(_numberOfConstraints);
+    for (int constraint = 0; constraint < _numberOfConstraints; constraint++)
+        _couplingMatrix[constraint].resize(_linkedConstraints[constraint].size());
+    _constraintDir.resize(_numberOfConstraints);
+    _rhs1.resize(_numberOfConstraints);
+    _rhs2.resize(_numberOfConstraints);
+    _solution.resize(_numberOfConstraints);
+}
+
+/**---------------------------------------------------------------------------------------
+
+ Solve the matrix equation
+
+ --------------------------------------------------------------------------------------- */
+
+void ReferenceLincsAlgorithm::solveMatrix() {
+    static const RealOpenMM zero = 0.0;
+    for (int iteration = 0; iteration < _numTerms; iteration++) {
+        vector<RealOpenMM>& rhs1 = (iteration%2 == 0 ? _rhs1 : _rhs2);
+        vector<RealOpenMM>& rhs2 = (iteration%2 == 0 ? _rhs2 : _rhs1);
+        for (int c1 = 0; c1 < _numberOfConstraints; c1++) {
+            rhs2[c1] = zero;
+            for (int j = 0; j < _linkedConstraints[c1].size(); j++) {
+                int c2 = _linkedConstraints[c1][j];
+                rhs2[c1] += _couplingMatrix[c1][j]*rhs1[c2];
+            }
+            _solution[c1] += rhs2[c1];
+        }
+    }
+}
+
+/**---------------------------------------------------------------------------------------
+
+ Update the atom position based on the solution to the matrix equation.
+
+ @param numberOfAtoms    number of atoms
+ @param atomCoordinates  atom coordinates
+ @param inverseMasses    1/mass
+
+ --------------------------------------------------------------------------------------- */
+
+void ReferenceLincsAlgorithm::updateAtomPositions(int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM* inverseMasses) {
+    for (int i = 0; i < _numberOfConstraints; i++) {
+        Vec3 delta(_sMatrix[i]*_solution[i]*_constraintDir[i][0],
+                   _sMatrix[i]*_solution[i]*_constraintDir[i][1],
+                   _sMatrix[i]*_solution[i]*_constraintDir[i][2]);
+        int atom1 = _atomIndices[i][0];
+        int atom2 = _atomIndices[i][1];
+        atomCoordinates[atom1][0] -= inverseMasses[atom1]*delta[0];
+        atomCoordinates[atom1][1] -= inverseMasses[atom1]*delta[1];
+        atomCoordinates[atom1][2] -= inverseMasses[atom1]*delta[2];
+        atomCoordinates[atom2][0] += inverseMasses[atom2]*delta[0];
+        atomCoordinates[atom2][1] += inverseMasses[atom2]*delta[1];
+        atomCoordinates[atom2][2] += inverseMasses[atom2]*delta[2];
+    }
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Apply Lincs algorithm
+
+   @param numberOfAtoms    number of atoms
+   @param atomCoordinates  atom coordinates
+   @param atomCoordinatesP atom coordinates prime
+   @param inverseMasses    1/mass
+
+   @return ReferenceDynamics::DefaultReturn if converge; else
+    return ReferenceDynamics::ErrorReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceLincsAlgorithm::apply( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                         RealOpenMM** atomCoordinatesP,
+                                         RealOpenMM* inverseMasses ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName = "\nReferenceLincsAlgorithm::apply";
+
+   static const RealOpenMM zero        =  0.0;
+   static const RealOpenMM one         =  1.0;
+   static const RealOpenMM two         =  2.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   if (_numberOfConstraints == 0)
+       return ReferenceDynamics::DefaultReturn;
+
+   if( !_hasInitialized )
+       initialize(numberOfAtoms, inverseMasses);
+
+   // Calculate the direction of each constraint, along with the initial RHS and solution vectors.
+
+   for (int i = 0; i < _numberOfConstraints; i++) {
+       int atom1 = _atomIndices[i][0];
+       int atom2 = _atomIndices[i][1];
+       _constraintDir[i] = Vec3(atomCoordinatesP[atom1][0]-atomCoordinatesP[atom2][0],
+                                atomCoordinatesP[atom1][1]-atomCoordinatesP[atom2][1],
+                                atomCoordinatesP[atom1][2]-atomCoordinatesP[atom2][2]);
+       RealOpenMM invLength = one/SQRT(_constraintDir[i].dot(_constraintDir[i]));
+       _constraintDir[i][0] *= invLength;
+       _constraintDir[i][1] *= invLength;
+       _constraintDir[i][2] *= invLength;
+       _rhs1[i] = _solution[i] = _sMatrix[i]*(one/invLength-_distance[i]);
+   }
+
+   // Build the coupling matrix.
+
+   for (int c1 = 0; c1 < _couplingMatrix.size(); c1++) {
+       Vec3& dir1 = _constraintDir[c1];
+       for (int j = 0; j < _couplingMatrix[c1].size(); j++) {
+           int c2 = _linkedConstraints[c1][j];
+           Vec3& dir2 = _constraintDir[c2];
+           if (_atomIndices[c1][0] == _atomIndices[c2][0] || _atomIndices[c1][1] == _atomIndices[c2][1])
+               _couplingMatrix[c1][j] = -inverseMasses[_atomIndices[c1][0]]*_sMatrix[c1]*dir1.dot(dir2)*_sMatrix[c2];
+           else
+               _couplingMatrix[c1][j] = inverseMasses[_atomIndices[c1][1]]*_sMatrix[c1]*dir1.dot(dir2)*_sMatrix[c2];
+       }
+   }
+
+   // Solve the matrix equation and update the positions.
+
+   solveMatrix();
+   updateAtomPositions(numberOfAtoms, atomCoordinatesP, inverseMasses);
+
+   // Correct for rotational lengthening.
+
+   for (int i = 0; i < _numberOfConstraints; i++) {
+       int atom1 = _atomIndices[i][0];
+       int atom2 = _atomIndices[i][1];
+       Vec3 delta(atomCoordinatesP[atom1][0]-atomCoordinatesP[atom2][0],
+                  atomCoordinatesP[atom1][1]-atomCoordinatesP[atom2][1],
+                  atomCoordinatesP[atom1][2]-atomCoordinatesP[atom2][2]);
+       RealOpenMM p2 = two*_distance[i]*_distance[i]-delta.dot(delta);
+       if (p2 < zero)
+           p2 = zero;
+       _rhs1[i] = _solution[i] = _sMatrix[i]*(_distance[i]-SQRT(p2));
+   }
+   solveMatrix();
+   updateAtomPositions(numberOfAtoms, atomCoordinatesP, inverseMasses);
+
+   return ReferenceDynamics::DefaultReturn;
+
+}
+

platforms/reference/src/SimTKReference/ReferenceLincsAlgorithm.h

+/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * 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.
+ */
+
+#ifndef __ReferenceLincsAlgorithm_H__
+#define __ReferenceLincsAlgorithm_H__
+
+#include "ReferenceConstraintAlgorithm.h"
+#include "Vec3.h"
+#include <vector>
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceLincsAlgorithm : public ReferenceConstraintAlgorithm {
+
+   protected:
+
+      int _numTerms;
+      int _numberOfConstraints;
+      int** _atomIndices;
+      RealOpenMM* _distance;
+
+      bool _hasInitialized;
+      std::vector<std::vector<int> > _linkedConstraints;
+      std::vector<RealOpenMM> _sMatrix;
+      std::vector<RealOpenMM> _rhs1;
+      std::vector<RealOpenMM> _rhs2;
+      std::vector<RealOpenMM> _solution;
+      std::vector<std::vector<RealOpenMM> > _couplingMatrix;
+      std::vector<OpenMM::Vec3> _constraintDir;
+
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the number of terms to use in the series expansion
+
+         @param numberOfAtoms    number of atoms
+         @param inverseMasses    1/mass
+
+         --------------------------------------------------------------------------------------- */
+
+      void initialize(int numberOfAtoms, RealOpenMM* inverseMasses);
+
+      /**---------------------------------------------------------------------------------------
+
+         Solve the matrix equation
+
+         --------------------------------------------------------------------------------------- */
+
+      void solveMatrix();
+
+      /**---------------------------------------------------------------------------------------
+
+         Update the atom position based on the solution to the matrix equation.
+
+         @param numberOfAtoms    number of atoms
+         @param atomCoordinates  atom coordinates
+         @param inverseMasses    1/mass
+
+         --------------------------------------------------------------------------------------- */
+
+      void updateAtomPositions(int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM* inverseMasses);
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         ReferenceLincsAlgorithm constructor
+
+         @param numberOfConstraints      number of constraints
+         @param atomIndices              atom indices for contraints
+         @param distance                 distances for constraints
+
+         --------------------------------------------------------------------------------------- */
+
+      ReferenceLincsAlgorithm( int numberOfConstraints, int** atomIndices, RealOpenMM* distance );
+
+      /**---------------------------------------------------------------------------------------
+
+         Get number of constraints
+
+         @return number of constraints
+
+         --------------------------------------------------------------------------------------- */
+
+      int getNumberOfConstraints( void ) const;
+
+      /**---------------------------------------------------------------------------------------
+
+         Get the number of terms to use in the series expansion
+
+         @return terms
+
+         --------------------------------------------------------------------------------------- */
+
+      int getNumTerms( void ) const;
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the number of terms to use in the series expansion
+
+         --------------------------------------------------------------------------------------- */
+
+      void setNumTerms( int terms );
+
+      /**---------------------------------------------------------------------------------------
+
+         Apply Lincs algorithm
+
+         @param numberOfAtoms    number of atoms
+         @param atomCoordinates  atom coordinates
+         @param atomCoordinatesP atom coordinates prime
+         @param inverseMasses    1/mass
+
+         @return ReferenceDynamics::DefaultReturn if converge; else
+          return ReferenceDynamics::ErrorReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int apply( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                       RealOpenMM** atomCoordinatesP, RealOpenMM* inverseMasses );
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // __ReferenceLincsAlgorithm_H__

platforms/reference/src/SimTKReference/ReferenceShakeAlgorithm.cpp

@@ -43,7 +43,7 @@
 
 ReferenceShakeAlgorithm::ReferenceShakeAlgorithm( int numberOfConstraints,
                                                   int** atomIndices,
-                                                  RealOpenMM** shakeParameters ){
+                                                  RealOpenMM* distance ){
 
    // ---------------------------------------------------------------------------------------
 
@@ -61,7 +61,7 @@ ReferenceShakeAlgorithm::ReferenceShakeAlgorithm( int numberOfConstraints,
 
    _numberOfConstraints        = numberOfConstraints;
    _atomIndices                = atomIndices;
-   _shakeParameters            = shakeParameters;
+   _distance                   = distance;
 
    _maximumNumberOfIterations  = 15;
    _tolerance                  = (RealOpenMM) 1.0e-04;
@@ -221,13 +221,13 @@ int ReferenceShakeAlgorithm::setTolerance( RealOpenMM tolerance ){
 
    --------------------------------------------------------------------------------------- */
 
-int ReferenceShakeAlgorithm::applyShake( int numberOfAtoms, RealOpenMM** atomCoordinates,
+int ReferenceShakeAlgorithm::apply( int numberOfAtoms, RealOpenMM** atomCoordinates,
                                          RealOpenMM** atomCoordinatesP,
                                          RealOpenMM* inverseMasses ){
 
    // ---------------------------------------------------------------------------------------
 
-   static const char* methodName = "\nReferenceShakeAlgorithm::applyShake";
+   static const char* methodName = "\nReferenceShakeAlgorithm::apply";
 
    static const RealOpenMM zero        =  0.0;
    static const RealOpenMM one         =  1.0;
@@ -298,7 +298,7 @@ int ReferenceShakeAlgorithm::applyShake( int numberOfAtoms, RealOpenMM** atomCoo
          }
 
          RealOpenMM rp2  = DOT3( rp_ij, rp_ij );
-         RealOpenMM dist2= _shakeParameters[ii][0]*_shakeParameters[ii][0];
+         RealOpenMM dist2= _distance[ii]*_distance[ii];
          RealOpenMM diff = dist2 - rp2;
          int iconv       = (int) ( FABS( diff )*distanceTolerance[ii] );
          if( iconv ){
@@ -409,9 +409,9 @@ int ReferenceShakeAlgorithm::reportShake( int numberOfAtoms, RealOpenMM** atomCo
       for( int jj = 0; jj < 3; jj++ ){
          rp2 += (atomCoordinates[atomI][jj] - atomCoordinates[atomJ][jj])*(atomCoordinates[atomI][jj] - atomCoordinates[atomJ][jj]); 
       }
-      RealOpenMM diff = FABS( rp2 - (_shakeParameters[ii][0]*_shakeParameters[ii][0]) );
+      RealOpenMM diff = FABS( rp2 - (_distance[ii]*_distance[ii]) );
       if( diff > tolerance ){
-         message << ii << " constraint violated: " << atomI << " " << atomJ << "] d=" << SQRT( rp2 ) << " " << rp2 << " d0=" << _shakeParameters[ii][0];
+         message << ii << " constraint violated: " << atomI << " " << atomJ << "] d=" << SQRT( rp2 ) << " " << rp2 << " d0=" << _distance[ii];
          message << " diff=" << diff;
          message << " [" << atomCoordinates[atomI][0] << " " << atomCoordinates[atomI][1] << " " << atomCoordinates[atomI][2] << "] ";
          message << " [" << atomCoordinates[atomJ][0] << " " << atomCoordinates[atomJ][1] << " " << atomCoordinates[atomJ][2] << "] ";

platforms/reference/src/SimTKReference/ReferenceShakeAlgorithm.h

@@ -25,9 +25,11 @@
 #ifndef __ReferenceShakeAlgorithm_H__
 #define __ReferenceShakeAlgorithm_H__
 
+#include "ReferenceConstraintAlgorithm.h"
+
 // ---------------------------------------------------------------------------------------
 
-class ReferenceShakeAlgorithm {
+class ReferenceShakeAlgorithm : public ReferenceConstraintAlgorithm {
 
    protected:
 
@@ -36,7 +38,7 @@ class ReferenceShakeAlgorithm {
 
       int _numberOfConstraints;
       int** _atomIndices;
-      RealOpenMM** _shakeParameters;
+      RealOpenMM* _distance;
 
       RealOpenMM** _r_ij;
       RealOpenMM* _d_ij2;
@@ -52,11 +54,11 @@ class ReferenceShakeAlgorithm {
       
          @param numberOfConstraints      number of constraints 
          @param atomIndices              atom indices for contraints
-         @param shakeParameters          Shake parameters
+         @param distance                 distances for constraints
       
          --------------------------------------------------------------------------------------- */
       
-      ReferenceShakeAlgorithm( int numberOfConstraints, int** atomIndices, RealOpenMM** shakeParameters );
+      ReferenceShakeAlgorithm( int numberOfConstraints, int** atomIndices, RealOpenMM* distance );
 
       /**---------------------------------------------------------------------------------------
       
@@ -145,7 +147,7 @@ class ReferenceShakeAlgorithm {
       
          --------------------------------------------------------------------------------------- */
       
-      int applyShake( int numberOfAtoms, RealOpenMM** atomCoordinates,
+      int apply( int numberOfAtoms, RealOpenMM** atomCoordinates,
                        RealOpenMM** atomCoordinatesP, RealOpenMM* inverseMasses );
       
       /**---------------------------------------------------------------------------------------

platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp

@@ -498,17 +498,17 @@ int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoo
 
    //writeState( numberOfAtoms, atomCoordinates, velocities, forces, masses, -1 , "Sd1" );
 
-   ReferenceShakeAlgorithm* referenceShakeAlgorithm = getReferenceShakeAlgorithm();
-   if( referenceShakeAlgorithm ){
+   ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
+   if( referenceConstraintAlgorithm ){
 
 /*
       std::stringstream message;
       message << methodName;
-      message << " calling shake1\n";
+      message << " calling constrain1\n";
       SimTKOpenMMLog::printMessage( message );
 */
 
-      referenceShakeAlgorithm->applyShake( numberOfAtoms, atomCoordinates, xPrime,
+      referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime,
                                            inverseMasses );
 
    }
@@ -556,16 +556,16 @@ int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoo
       SimTKOpenMMLog::printMessage( message );
    }
 
-   if( referenceShakeAlgorithm ){
+   if( referenceConstraintAlgorithm ){
 
 /*
       std::stringstream message;
       message << methodName;
-      message << " calling shake2\n";
+      message << " calling constrain2\n";
       SimTKOpenMMLog::printMessage( message );
 */
 
-      referenceShakeAlgorithm->applyShake( numberOfAtoms, atomCoordinates, xPrime,
+      referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime,
                                            inverseMasses );
    }
 

platforms/reference/src/SimTKReference/ReferenceVerletDynamics.cpp

@@ -170,9 +170,9 @@ int ReferenceVerletDynamics::update( int numberOfAtoms, RealOpenMM** atomCoordin
            xPrime[i][j] = atomCoordinates[i][j] + velocities[i][j]*getDeltaT();
        }
    }
-   ReferenceShakeAlgorithm* referenceShakeAlgorithm = getReferenceShakeAlgorithm();
-   if( referenceShakeAlgorithm )
-      referenceShakeAlgorithm->applyShake( numberOfAtoms, atomCoordinates, xPrime, inverseMasses );
+   ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
+   if( referenceConstraintAlgorithm )
+      referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime, inverseMasses );
    
    // Update the positions and velocities.