Began implementing new mixed precision model that does integration in double...

Began implementing new mixed precision model that does integration in double precision and force evaluation in single precision

platforms/cuda/src/CudaContext.cpp

@@ -67,22 +67,22 @@ bool CudaContext::hasInitializedCuda = false;
 CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlockingSync, const string& precision, const string& compiler,
         const string& tempDir, CudaPlatform::PlatformData& platformData) : system(system), compiler(compiler),
         time(0.0), platformData(platformData), stepCount(0), computeForceCount(0), contextIsValid(false), atomsWereReordered(false), pinnedBuffer(NULL), posq(NULL),
-        velm(NULL), force(NULL), energyBuffer(NULL), integration(NULL), expression(NULL), bonded(NULL), nonbonded(NULL), thread(NULL) {
+        posqCorrection(NULL), velm(NULL), force(NULL), energyBuffer(NULL), integration(NULL), expression(NULL), bonded(NULL), nonbonded(NULL), thread(NULL) {
     if (!hasInitializedCuda) {
         CHECK_RESULT2(cuInit(0), "Error initializing CUDA");
         hasInitializedCuda = true;
     }
     if (precision == "single") {
         useDoublePrecision = false;
-        accumulateInDouble = false;
+        useMixedPrecision = false;
     }
     else if (precision == "mixed") {
         useDoublePrecision = false;
-        accumulateInDouble = true;
+        useMixedPrecision = true;
     }
     else if (precision == "double") {
         useDoublePrecision = true;
-        accumulateInDouble = true;
+        useMixedPrecision = false;
     }
     else
         throw OpenMMException("Illegal value for CudaPrecision: "+precision);
@@ -150,16 +150,37 @@ CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlocking
         compilationDefines["make_real2"] = "make_double2";
         compilationDefines["make_real3"] = "make_double3";
         compilationDefines["make_real4"] = "make_double4";
+        compilationDefines["make_mixed2"] = "make_double2";
+        compilationDefines["make_mixed3"] = "make_double3";
+        compilationDefines["make_mixed4"] = "make_double4";
         energyBuffer = CudaArray::create<double>(*this, numEnergyBuffers, "energyBuffer");
         int pinnedBufferSize = max(paddedNumAtoms*4, numEnergyBuffers);
         CHECK_RESULT(cuMemHostAlloc(&pinnedBuffer, pinnedBufferSize*sizeof(double), 0));
     }
+    else if (useMixedPrecision) {
+        posq = CudaArray::create<float4>(*this, paddedNumAtoms, "posq");
+        posqCorrection = CudaArray::create<float4>(*this, paddedNumAtoms, "posqCorrection");
+        velm = CudaArray::create<double4>(*this, paddedNumAtoms, "velm");
+        compilationDefines["USE_MIXED_PRECISION"] = "1";
+        compilationDefines["make_real2"] = "make_float2";
+        compilationDefines["make_real3"] = "make_float3";
+        compilationDefines["make_real4"] = "make_float4";
+        compilationDefines["make_mixed2"] = "make_double2";
+        compilationDefines["make_mixed3"] = "make_double3";
+        compilationDefines["make_mixed4"] = "make_double4";
+        energyBuffer = CudaArray::create<float>(*this, numEnergyBuffers, "energyBuffer");
+        int pinnedBufferSize = max(paddedNumAtoms*4, numEnergyBuffers);
+        CHECK_RESULT(cuMemHostAlloc(&pinnedBuffer, pinnedBufferSize*sizeof(double), 0));
+    }
     else {
         posq = CudaArray::create<float4>(*this, paddedNumAtoms, "posq");
         velm = CudaArray::create<float4>(*this, paddedNumAtoms, "velm");
         compilationDefines["make_real2"] = "make_float2";
         compilationDefines["make_real3"] = "make_float3";
         compilationDefines["make_real4"] = "make_float4";
+        compilationDefines["make_mixed2"] = "make_float2";
+        compilationDefines["make_mixed3"] = "make_float3";
+        compilationDefines["make_mixed4"] = "make_float4";
         energyBuffer = CudaArray::create<float>(*this, numEnergyBuffers, "energyBuffer");
         int pinnedBufferSize = max(paddedNumAtoms*6, numEnergyBuffers);
         CHECK_RESULT(cuMemHostAlloc(&pinnedBuffer, pinnedBufferSize*sizeof(float), 0));
@@ -211,6 +232,8 @@ CudaContext::~CudaContext() {
         cuMemFreeHost(pinnedBuffer);
     if (posq != NULL)
         delete posq;
+    if (posqCorrection != NULL)
+        delete posqCorrection;
     if (velm != NULL)
         delete velm;
     if (force != NULL)
@@ -237,7 +260,7 @@ void CudaContext::initialize() {
     string errorMessage = "Error initializing Context";
     for (int i = 0; i < numAtoms; i++) {
         double mass = system.getParticleMass(i);
-        if (useDoublePrecision)
+        if (useDoublePrecision || useMixedPrecision)
             ((double4*) pinnedBuffer)[i] = make_double4(0.0, 0.0, 0.0, mass == 0.0 ? 0.0 : 1.0/mass);
         else
             ((float4*) pinnedBuffer)[i] = make_float4(0.0f, 0.0f, 0.0f, mass == 0.0 ? 0.0f : (float) (1.0/mass));
@@ -308,6 +331,18 @@ CUmodule CudaContext::createModule(const string source, const map<string, string
         src << "typedef float3 real3;\n";
         src << "typedef float4 real4;\n";
     }
+    if (useDoublePrecision || useMixedPrecision) {
+        src << "typedef double mixed;\n";
+        src << "typedef double2 mixed2;\n";
+        src << "typedef double3 mixed3;\n";
+        src << "typedef double4 mixed4;\n";
+    }
+    else {
+        src << "typedef float mixed;\n";
+        src << "typedef float2 mixed2;\n";
+        src << "typedef float3 mixed3;\n";
+        src << "typedef float4 mixed4;\n";
+    }
     for (map<string, string>::const_iterator iter = defines.begin(); iter != defines.end(); ++iter) {
         src << "#define " << iter->first;
         if (!iter->second.empty())
@@ -789,6 +824,22 @@ void CudaContext::validateMolecules() {
         posq->upload(newPosq);
         velm->upload(newVelm);
     }
+    else if (useMixedPrecision) {
+        vector<float4> oldPosq(paddedNumAtoms);
+        vector<float4> newPosq(paddedNumAtoms);
+        vector<double4> oldVelm(paddedNumAtoms);
+        vector<double4> newVelm(paddedNumAtoms);
+        posq->download(oldPosq);
+        velm->download(oldVelm);
+        for (int i = 0; i < numAtoms; i++) {
+            int index = atomIndex[i];
+            newPosq[index] = oldPosq[i];
+            newVelm[index] = oldVelm[i];
+            newCellOffsets[index] = posCellOffsets[i];
+        }
+        posq->upload(newPosq);
+        velm->upload(newVelm);
+    }
     else {
         vector<float4> oldPosq(paddedNumAtoms);
         vector<float4> newPosq(paddedNumAtoms);
@@ -822,19 +873,24 @@ void CudaContext::reorderAtoms(bool enforcePeriodic) {
         validateMolecules();
     atomsWereReordered = true;
     if (useDoublePrecision)
-        reorderAtomsImpl<double, double4>(enforcePeriodic);
+        reorderAtomsImpl<double, double4, double, double4>(enforcePeriodic);
+    else if (useMixedPrecision)
+        reorderAtomsImpl<float, float4, double, double4>(enforcePeriodic);
     else
-        reorderAtomsImpl<float, float4>(enforcePeriodic);
+        reorderAtomsImpl<float, float4, float, float4>(enforcePeriodic);
 }
 
-template <class Real, class Real4>
+template <class Real, class Real4, class Mixed, class Mixed4>
 void CudaContext::reorderAtomsImpl(bool enforcePeriodic) {
     // Find the range of positions and the number of bins along each axis.
 
     vector<Real4> oldPosq(paddedNumAtoms);
-    vector<Real4> oldVelm(paddedNumAtoms);
+    vector<Real4> oldPosqCorrection(paddedNumAtoms);
+    vector<Mixed4> oldVelm(paddedNumAtoms);
     posq->download(oldPosq);
     velm->download(oldVelm);
+    if (useMixedPrecision)
+        posqCorrection->download(oldPosqCorrection);
     Real minx = oldPosq[0].x, maxx = oldPosq[0].x;
     Real miny = oldPosq[0].y, maxy = oldPosq[0].y;
     Real minz = oldPosq[0].z, maxz = oldPosq[0].z;
@@ -860,7 +916,8 @@ void CudaContext::reorderAtomsImpl(bool enforcePeriodic) {
 
     vector<int> originalIndex(numAtoms);
     vector<Real4> newPosq(paddedNumAtoms);
-    vector<Real4> newVelm(paddedNumAtoms);
+    vector<Real4> newPosqCorrection(paddedNumAtoms);
+    vector<Mixed4> newVelm(paddedNumAtoms);
     vector<int4> newCellOffsets(numAtoms);
     for (int group = 0; group < (int) moleculeGroups.size(); group++) {
         // Find the center of each molecule.
@@ -959,6 +1016,8 @@ void CudaContext::reorderAtomsImpl(bool enforcePeriodic) {
                 int newIndex = mol.offsets[i]+atoms[j];
                 originalIndex[newIndex] = atomIndex[oldIndex];
                 newPosq[newIndex] = oldPosq[oldIndex];
+                if (useMixedPrecision)
+                    newPosqCorrection[newIndex] = oldPosqCorrection[oldIndex];
                 newVelm[newIndex] = oldVelm[oldIndex];
                 newCellOffsets[newIndex] = posCellOffsets[oldIndex];
             }
@@ -972,6 +1031,8 @@ void CudaContext::reorderAtomsImpl(bool enforcePeriodic) {
         posCellOffsets[i] = newCellOffsets[i];
     }
     posq->upload(newPosq);
+    if (useMixedPrecision)
+        posqCorrection->upload(newPosqCorrection);
     velm->upload(newVelm);
     atomIndexDevice->upload(atomIndex);
     for (int i = 0; i < (int) reorderListeners.size(); i++)

platforms/cuda/src/CudaContext.h

@@ -135,6 +135,12 @@ public:
     CudaArray& getPosq() {
         return *posq;
     }
+    /**
+     * Get the array which contains a correction to the position of each atom.  This only exists if getUseMixedPrecision() returns true.
+     */
+    CudaArray& getPosqCorrection() {
+        return *posqCorrection;
+    }
     /**
      * Get the array which contains the velocity (the xyz components) and inverse mass (the w component) of each atom.
      */
@@ -314,10 +320,10 @@ public:
         return useDoublePrecision;
     }
     /**
-     * Get whether accumulation is being done in double precision.
+     * Get whether mixed precision is being used.
      */
-    bool getAccumulateInDouble() {
-        return accumulateInDouble;
+    bool getUseMixedPrecision() {
+        return useMixedPrecision;
     }
     /**
      * Convert a number to a string in a format suitable for including in a kernel.
@@ -455,7 +461,7 @@ private:
     /**
      * This is the internal implementation of reorderAtoms(), templatized by the numerical precision in use.
      */
-    template <class Real, class Real4>
+    template <class Real, class Real4, class Mixed, class Mixed4>
     void reorderAtomsImpl(bool enforcePeriodic);
     static bool hasInitializedCuda;
     const System& system;
@@ -469,7 +475,7 @@ private:
     int paddedNumAtoms;
     int numAtomBlocks;
     int numThreadBlocks;
-    bool useBlockingSync, useDoublePrecision, accumulateInDouble, contextIsValid, atomsWereReordered, moleculesInvalid;
+    bool useBlockingSync, useDoublePrecision, useMixedPrecision, contextIsValid, atomsWereReordered, moleculesInvalid;
     std::string compiler, tempDir, gpuArchitecture;
     float4 periodicBoxSizeFloat, invPeriodicBoxSizeFloat;
     double4 periodicBoxSize, invPeriodicBoxSize;
@@ -489,6 +495,7 @@ private:
     std::vector<int4> posCellOffsets;
     void* pinnedBuffer;
     CudaArray* posq;
+    CudaArray* posqCorrection;
     CudaArray* velm;
     CudaArray* force;
     CudaArray* energyBuffer;

platforms/cuda/src/CudaIntegrationUtilities.cpp

@@ -104,7 +104,7 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
         vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL) {
     // Create workspace arrays.
 
-    if (context.getUseDoublePrecision()) {
+    if (context.getUseDoublePrecision() || context.getUseMixedPrecision()) {
         posDelta = CudaArray::create<double4>(context, context.getPaddedNumAtoms(), "posDelta");
         vector<double4> deltas(posDelta->getSize(), make_double4(0.0, 0.0, 0.0, 0.0));
         posDelta->upload(deltas);
@@ -473,7 +473,7 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
         vector<int> atomConstraintsVec(ccmaAtomConstraints->getSize());
         vector<int> numAtomConstraintsVec(ccmaNumAtomConstraints->getSize());
         vector<int> constraintMatrixColumnVec(ccmaConstraintMatrixColumn->getSize());
-        if (context.getUseDoublePrecision()) {
+        if (context.getUseDoublePrecision() || context.getUseMixedPrecision()) {
             ccmaDistance = CudaArray::create<double4>(context, numCCMA, "CcmaDistance");
             ccmaDelta1 = CudaArray::create<double>(context, numCCMA, "CcmaDelta1");
             ccmaDelta2 = CudaArray::create<double>(context, numCCMA, "CcmaDelta2");
@@ -717,23 +717,24 @@ void CudaIntegrationUtilities::applyConstraints(bool constrainVelocities, double
         ccmaForceKernel = ccmaPosForceKernel;
     }
     float floatTol = (float) tol;
-    void* tolPointer = (context.getUseDoublePrecision() ? (void*) &tol : (void*) &floatTol);
+    void* tolPointer = (context.getUseDoublePrecision() || context.getUseMixedPrecision() ? (void*) &tol : (void*) &floatTol);
+    CUdeviceptr posCorrection = (context.getUseMixedPrecision() ? context.getPosqCorrection().getDevicePointer() : 0);
     if (settleAtoms != NULL) {
         int numClusters = settleAtoms->getSize();
-        void* args[] = {&numClusters, tolPointer, &context.getPosq().getDevicePointer(),
+        void* args[] = {&numClusters, tolPointer, &context.getPosq().getDevicePointer(), &posCorrection,
                 &posDelta->getDevicePointer(), &context.getVelm().getDevicePointer(),
                 &settleAtoms->getDevicePointer(), &settleParams->getDevicePointer()};
         context.executeKernel(settleKernel, args, settleAtoms->getSize());
     }
     if (shakeAtoms != NULL) {
         int numClusters = shakeAtoms->getSize();
-        void* args[] = {&numClusters, tolPointer, &context.getPosq().getDevicePointer(),
+        void* args[] = {&numClusters, tolPointer, &context.getPosq().getDevicePointer(), &posCorrection,
                 constrainVelocities ? &context.getVelm().getDevicePointer() : &posDelta->getDevicePointer(),
                 &shakeAtoms->getDevicePointer(), &shakeParams->getDevicePointer()};
         context.executeKernel(shakeKernel, args, shakeAtoms->getSize());
     }
     if (ccmaAtoms != NULL) {
-        void* directionsArgs[] = {&ccmaAtoms->getDevicePointer(), &ccmaDistance->getDevicePointer(), &context.getPosq().getDevicePointer()};
+        void* directionsArgs[] = {&ccmaAtoms->getDevicePointer(), &ccmaDistance->getDevicePointer(), &context.getPosq().getDevicePointer(), &posCorrection};
         context.executeKernel(ccmaDirectionsKernel, directionsArgs, ccmaAtoms->getSize());
         int i;
         void* forceArgs[] = {&ccmaAtoms->getDevicePointer(), &ccmaDistance->getDevicePointer(),
@@ -768,7 +769,8 @@ void CudaIntegrationUtilities::applyConstraints(bool constrainVelocities, double
 
 void CudaIntegrationUtilities::computeVirtualSites() {
     if (numVsites > 0) {
-        void* args[] = {&context.getPosq().getDevicePointer(), &vsite2AvgAtoms->getDevicePointer(), &vsite2AvgWeights->getDevicePointer(),
+        CUdeviceptr posCorrection = (context.getUseMixedPrecision() ? context.getPosqCorrection().getDevicePointer() : 0);
+        void* args[] = {&context.getPosq().getDevicePointer(), &posCorrection, &vsite2AvgAtoms->getDevicePointer(), &vsite2AvgWeights->getDevicePointer(),
                 &vsite3AvgAtoms->getDevicePointer(), &vsite3AvgWeights->getDevicePointer(),
                 &vsiteOutOfPlaneAtoms->getDevicePointer(), &vsiteOutOfPlaneWeights->getDevicePointer()};
         context.executeKernel(vsitePositionKernel, args, numVsites);
@@ -777,7 +779,8 @@ void CudaIntegrationUtilities::computeVirtualSites() {
 
 void CudaIntegrationUtilities::distributeForcesFromVirtualSites() {
     if (numVsites > 0) {
-        void* args[] = {&context.getPosq().getDevicePointer(), &context.getForce().getDevicePointer(),
+        CUdeviceptr posCorrection = (context.getUseMixedPrecision() ? context.getPosqCorrection().getDevicePointer() : 0);
+        void* args[] = {&context.getPosq().getDevicePointer(), &posCorrection, &context.getForce().getDevicePointer(),
                 &vsite2AvgAtoms->getDevicePointer(), &vsite2AvgWeights->getDevicePointer(),
                 &vsite3AvgAtoms->getDevicePointer(), &vsite3AvgWeights->getDevicePointer(),
                 &vsiteOutOfPlaneAtoms->getDevicePointer(), &vsiteOutOfPlaneWeights->getDevicePointer()};

platforms/cuda/src/CudaKernels.cpp

@@ -148,6 +148,18 @@ void CudaUpdateStateDataKernel::getPositions(ContextImpl& context, vector<Vec3>&
             positions[order[i]] = Vec3(pos.x-offset.x*periodicBoxSize.x, pos.y-offset.y*periodicBoxSize.y, pos.z-offset.z*periodicBoxSize.z);
         }
     }
+    else if (cu.getUseMixedPrecision()) {
+        float4* posq = (float4*) cu.getPinnedBuffer();
+        vector<float4> posCorrection;
+        cu.getPosq().download(posq);
+        cu.getPosqCorrection().download(posCorrection);
+        for (int i = 0; i < numParticles; ++i) {
+            float4 pos1 = posq[i];
+            float4 pos2 = posCorrection[i];
+            int4 offset = cu.getPosCellOffsets()[i];
+            positions[order[i]] = Vec3((double)pos1.x+(double)pos2.x-offset.x*periodicBoxSize.x, (double)pos1.y+(double)pos2.y-offset.y*periodicBoxSize.y, (double)pos1.z+(double)pos2.z-offset.z*periodicBoxSize.z);
+        }
+    }
     else {
         float4* posq = (float4*) cu.getPinnedBuffer();
         cu.getPosq().download(posq);
@@ -183,14 +195,28 @@ void CudaUpdateStateDataKernel::setPositions(ContextImpl& context, const vector<
         for (int i = 0; i < numParticles; ++i) {
             float4& pos = posq[i];
             const Vec3& p = positions[order[i]];
-            pos.x = p[0];
-            pos.y = p[1];
-            pos.z = p[2];
+            pos.x = (float) p[0];
+            pos.y = (float) p[1];
+            pos.z = (float) p[2];
         }
         for (int i = numParticles; i < cu.getPaddedNumAtoms(); i++)
             posq[i] = make_float4(0.0, 0.0, 0.0, 0.0);
         cu.getPosq().upload(posq);
     }
+    if (cu.getUseMixedPrecision()) {
+        float4* posCorrection = (float4*) cu.getPinnedBuffer();
+        for (int i = 0; i < numParticles; ++i) {
+            float4& c = posCorrection[i];
+            const Vec3& p = positions[order[i]];
+            c.x = (float) (p[0]-(float)p[0]);
+            c.y = (float) (p[1]-(float)p[1]);
+            c.z = (float) (p[2]-(float)p[2]);
+            c.w = 0;
+        }
+        for (int i = numParticles; i < cu.getPaddedNumAtoms(); i++)
+            posCorrection[i] = make_float4(0.0, 0.0, 0.0, 0.0);
+        cu.getPosqCorrection().upload(posCorrection);
+    }
     for (int i = 0; i < (int) cu.getPosCellOffsets().size(); i++)
         cu.getPosCellOffsets()[i] = make_int4(0, 0, 0, 0);
 }
@@ -200,7 +226,7 @@ void CudaUpdateStateDataKernel::getVelocities(ContextImpl& context, vector<Vec3>
     const vector<int>& order = cu.getAtomIndex();
     int numParticles = context.getSystem().getNumParticles();
     velocities.resize(numParticles);
-    if (cu.getUseDoublePrecision()) {
+    if (cu.getUseDoublePrecision() || cu.getUseMixedPrecision()) {
         double4* velm = (double4*) cu.getPinnedBuffer();
         cu.getVelm().download(velm);
         for (int i = 0; i < numParticles; ++i) {
@@ -224,7 +250,7 @@ void CudaUpdateStateDataKernel::setVelocities(ContextImpl& context, const vector
     cu.setAsCurrent();
     const vector<int>& order = cu.getAtomIndex();
     int numParticles = context.getSystem().getNumParticles();
-    if (cu.getUseDoublePrecision()) {
+    if (cu.getUseDoublePrecision() || cu.getUseMixedPrecision()) {
         double4* velm = (double4*) cu.getPinnedBuffer();
         cu.getVelm().download(velm);
         for (int i = 0; i < numParticles; ++i) {
@@ -290,12 +316,11 @@ void CudaUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostream&
     stream.write((char*) &stepCount, sizeof(int));
     int computeForceCount = cu.getComputeForceCount();
     stream.write((char*) &computeForceCount, sizeof(int));
-    int bufferSize = cu.getPaddedNumAtoms()*(cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4));
     char* buffer = (char*) cu.getPinnedBuffer();
     cu.getPosq().download(buffer);
-    stream.write(buffer, bufferSize);
+    stream.write(buffer, cu.getPosq().getSize()*cu.getPosq().getElementSize());
     cu.getVelm().download(buffer);
-    stream.write(buffer, bufferSize);
+    stream.write(buffer, cu.getVelm().getSize()*cu.getVelm().getElementSize());
     stream.write((char*) &cu.getAtomIndex()[0], sizeof(int)*cu.getAtomIndex().size());
     stream.write((char*) &cu.getPosCellOffsets()[0], sizeof(int4)*cu.getPosCellOffsets().size());
     double4 box = cu.getPeriodicBoxSize();
@@ -321,11 +346,10 @@ void CudaUpdateStateDataKernel::loadCheckpoint(ContextImpl& context, istream& st
         contexts[i]->setStepCount(stepCount);
         contexts[i]->setComputeForceCount(computeForceCount);
     }
-    int bufferSize = cu.getPaddedNumAtoms()*(cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4));
     char* buffer = (char*) cu.getPinnedBuffer();
-    stream.read(buffer, bufferSize);
+    stream.read(buffer, cu.getPosq().getSize()*cu.getPosq().getElementSize());
     cu.getPosq().upload(buffer);
-    stream.read(buffer, bufferSize);
+    stream.read(buffer, cu.getVelm().getSize()*cu.getVelm().getElementSize());
     cu.getVelm().upload(buffer);
     stream.read((char*) &cu.getAtomIndex()[0], sizeof(int)*cu.getAtomIndex().size());
     cu.getAtomIndexArray().upload(cu.getAtomIndex());
@@ -2016,7 +2040,7 @@ double CudaCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeFor
         defines["FORCE_WORK_GROUP_SIZE"] = cu.intToString(nb.getForceThreadBlockSize());
         map<string, string> replacements;
         stringstream defineAccum;
-        if (cu.getAccumulateInDouble()) {
+        if (cu.getUseMixedPrecision()) {
             defineAccum << "typedef double accum;\n";
             defineAccum << "typedef double4 accum4;\n";
             defines["make_accum4"] = "make_double4";
@@ -2531,7 +2555,7 @@ void CudaCalcCustomGBForceKernel::initialize(const System& system, const CustomG
         replacements["STORE_DERIVATIVES_2"] = storeDerivs2.str();
         map<string, string> defines;
         stringstream defineAccum;
-        if (cu.getAccumulateInDouble()) {
+        if (cu.getUseMixedPrecision()) {
             defineAccum << "typedef double accum;\n";
             defineAccum << "typedef double3 accum3;\n";
             defines["make_accum3"] = "make_double3";
@@ -3971,7 +3995,7 @@ void CudaIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIn
     int numAtoms = cu.getNumAtoms();
     double dt = integrator.getStepSize();
     if (dt != prevStepSize) {
-        if (cu.getUseDoublePrecision()) {
+        if (cu.getUseDoublePrecision() || cu.getUseMixedPrecision()) {
             vector<double2> stepSizeVec(1);
             stepSizeVec[0] = make_double2(dt, dt);
             cu.getIntegrationUtilities().getStepSize().upload(stepSizeVec);
@@ -3986,7 +4010,8 @@ void CudaIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIn
 
     // Call the first integration kernel.
 
-    void* args1[] = {&cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(),
+    CUdeviceptr posCorrection = (cu.getUseMixedPrecision() ? cu.getPosqCorrection().getDevicePointer() : 0);
+    void* args1[] = {&cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(), &posCorrection,
             &cu.getVelm().getDevicePointer(), &cu.getForce().getDevicePointer(), &integration.getPosDelta().getDevicePointer()};
     cu.executeKernel(kernel1, args1, numAtoms);
 
@@ -3996,7 +4021,7 @@ void CudaIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIn
 
     // Call the second integration kernel.
 
-    void* args2[] = {&cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(),
+    void* args2[] = {&cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(), &posCorrection,
             &cu.getVelm().getDevicePointer(), &integration.getPosDelta().getDevicePointer()};
     cu.executeKernel(kernel2, args2, numAtoms);
     integration.computeVirtualSites();
@@ -4023,7 +4048,7 @@ void CudaIntegrateLangevinStepKernel::initialize(const System& system, const Lan
     CUmodule module = cu.createModule(CudaKernelSources::langevin, defines, "");
     kernel1 = cu.getKernel(module, "integrateLangevinPart1");
     kernel2 = cu.getKernel(module, "integrateLangevinPart2");
-    params = new CudaArray(cu, 3, cu.getUseDoublePrecision() ? sizeof(double) : sizeof(float), "langevinParams");
+    params = new CudaArray(cu, 3, cu.getUseDoublePrecision() || cu.getUseMixedPrecision() ? sizeof(double) : sizeof(float), "langevinParams");
     prevStepSize = -1.0;
 }
 
@@ -4042,7 +4067,7 @@ void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const Langev
         double vscale = exp(-stepSize/tau);
         double fscale = (1-vscale)*tau;
         double noisescale = sqrt(2*kT/tau)*sqrt(0.5*(1-vscale*vscale)*tau);
-        if (cu.getUseDoublePrecision()) {
+        if (cu.getUseDoublePrecision() || cu.getUseMixedPrecision()) {
             vector<double> p(params->getSize());
             p[0] = vscale;
             p[1] = fscale;
@@ -4078,7 +4103,8 @@ void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const Langev
 
     // Call the second integration kernel.
 
-    void* args2[] = {&cu.getPosq().getDevicePointer(), &integration.getPosDelta().getDevicePointer(),
+    CUdeviceptr posCorrection = (cu.getUseMixedPrecision() ? cu.getPosqCorrection().getDevicePointer() : 0);
+    void* args2[] = {&cu.getPosq().getDevicePointer(), &posCorrection, &integration.getPosDelta().getDevicePointer(),
             &cu.getVelm().getDevicePointer(), &integration.getStepSize().getDevicePointer()};
     cu.executeKernel(kernel2, args2, numAtoms);
     integration.computeVirtualSites();
@@ -4172,16 +4198,18 @@ double CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, cons
     float maxStepSizeFloat = (float) maxStepSize;
     double tol = integrator.getErrorTolerance();
     float tolFloat = (float) tol;
-    void* argsSelect[] = {cu.getUseDoublePrecision() ? (void*) &maxStepSize : (void*) &maxStepSizeFloat,
-            cu.getUseDoublePrecision() ? (void*) &tol : (void*) &tolFloat,
+    bool useDouble = cu.getUseDoublePrecision() || cu.getUseMixedPrecision();
+    void* argsSelect[] = {useDouble ? (void*) &maxStepSize : (void*) &maxStepSizeFloat,
+            useDouble ? (void*) &tol : (void*) &tolFloat,
             &cu.getIntegrationUtilities().getStepSize().getDevicePointer(),
             &cu.getVelm().getDevicePointer(), &cu.getForce().getDevicePointer()};
-    int sharedSize = blockSize*(cu.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
+    int sharedSize = blockSize*(useDouble ? sizeof(double) : sizeof(float));
     cu.executeKernel(selectSizeKernel, argsSelect, blockSize, blockSize, sharedSize);
 
     // Call the first integration kernel.
 
-    void* args1[] = {&cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(),
+    CUdeviceptr posCorrection = (cu.getUseMixedPrecision() ? cu.getPosqCorrection().getDevicePointer() : 0);
+    void* args1[] = {&cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(), &posCorrection,
             &cu.getVelm().getDevicePointer(), &cu.getForce().getDevicePointer(), &integration.getPosDelta().getDevicePointer()};
     cu.executeKernel(kernel1, args1, numAtoms);
 
@@ -4191,7 +4219,7 @@ double CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, cons
 
     // Call the second integration kernel.
 
-    void* args2[] = {&cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(),
+    void* args2[] = {&cu.getIntegrationUtilities().getStepSize().getDevicePointer(), &cu.getPosq().getDevicePointer(), &posCorrection,
             &cu.getVelm().getDevicePointer(), &integration.getPosDelta().getDevicePointer()};
     cu.executeKernel(kernel2, args2, numAtoms);
     integration.computeVirtualSites();
@@ -4199,7 +4227,7 @@ double CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, cons
     // Update the time and step count.
 
     double dt, time;
-    if (cu.getUseDoublePrecision()) {
+    if (useDouble) {
         double2 stepSize;
         cu.getIntegrationUtilities().getStepSize().download(&stepSize);
         dt = stepSize.y;
@@ -4237,7 +4265,7 @@ void CudaIntegrateVariableLangevinStepKernel::initialize(const System& system, c
     kernel1 = cu.getKernel(module, "integrateLangevinPart1");
     kernel2 = cu.getKernel(module, "integrateLangevinPart2");
     selectSizeKernel = cu.getKernel(module, "selectLangevinStepSize");
-    params = CudaArray::create<float>(cu, 3, "langevinParams");
+    params = new CudaArray(cu, 3, cu.getUseDoublePrecision() || cu.getUseMixedPrecision() ? sizeof(double) : sizeof(float), "langevinParams");
     blockSize = min(256, system.getNumParticles());
     blockSize = max(blockSize, params->getSize());
 }
@@ -4257,13 +4285,14 @@ double CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, co
     float tauFloat = (float) tau;
     double kT = BOLTZ*integrator.getTemperature();
     float kTFloat = (float) kT;
-    void* argsSelect[] = {cu.getUseDoublePrecision() ? (void*) &maxStepSize : (void*) &maxStepSizeFloat,
-            cu.getUseDoublePrecision() ? (void*) &tol : (void*) &tolFloat,
-            cu.getUseDoublePrecision() ? (void*) &tau : (void*) &tauFloat,
-            cu.getUseDoublePrecision() ? (void*) &kT : (void*) &kTFloat,
+    bool useDouble = cu.getUseDoublePrecision() || cu.getUseMixedPrecision();
+    void* argsSelect[] = {useDouble ? (void*) &maxStepSize : (void*) &maxStepSizeFloat,
+            useDouble ? (void*) &tol : (void*) &tolFloat,
+            useDouble ? (void*) &tau : (void*) &tauFloat,
+            useDouble ? (void*) &kT : (void*) &kTFloat,
             &cu.getIntegrationUtilities().getStepSize().getDevicePointer(),
             &cu.getVelm().getDevicePointer(), &cu.getForce().getDevicePointer(), &params->getDevicePointer()};
-    int sharedSize = blockSize*(cu.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
+    int sharedSize = blockSize*(useDouble ? sizeof(double) : sizeof(float));
     cu.executeKernel(selectSizeKernel, argsSelect, blockSize, blockSize, sharedSize);
 
     // Call the first integration kernel.
@@ -4279,7 +4308,8 @@ double CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, co
 
     // Call the second integration kernel.
 
-    void* args2[] = {&cu.getPosq().getDevicePointer(), &integration.getPosDelta().getDevicePointer(),
+    CUdeviceptr posCorrection = (cu.getUseMixedPrecision() ? cu.getPosqCorrection().getDevicePointer() : 0);
+    void* args2[] = {&cu.getPosq().getDevicePointer(), &posCorrection, &integration.getPosDelta().getDevicePointer(),
             &cu.getVelm().getDevicePointer(), &integration.getStepSize().getDevicePointer()};
     cu.executeKernel(kernel2, args2, numAtoms);
     integration.computeVirtualSites();
@@ -4287,7 +4317,7 @@ double CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, co
     // Update the time and step count.
 
     double dt, time;
-    if (cu.getUseDoublePrecision()) {
+    if (useDouble) {
         double2 stepSize;
         cu.getIntegrationUtilities().getStepSize().download(&stepSize);
         dt = stepSize.y;
@@ -5099,7 +5129,7 @@ double CudaCalcKineticEnergyKernel::execute(ContextImpl& context) {
 
     const vector<int>& order = cu.getAtomIndex();
     double energy = 0.0;
-    if (cu.getUseDoublePrecision()) {
+    if (cu.getUseDoublePrecision() || cu.getUseMixedPrecision()) {
         double4* velm = (double4*) cu.getPinnedBuffer();
         cu.getVelm().download(velm);
         for (size_t i = 0; i < masses.size(); ++i) {

platforms/cuda/src/CudaNonbondedUtilities.cpp

@@ -446,12 +446,12 @@ CUfunction CudaNonbondedUtilities::createInteractionKernel(const string& source,
     defines["NUM_ATOMS"] = context.intToString(context.getNumAtoms());
     defines["PADDED_NUM_ATOMS"] = context.intToString(context.getPaddedNumAtoms());
     defines["NUM_BLOCKS"] = context.intToString(context.getNumAtomBlocks());
-    if ((localDataSize/4)%2 == 0 && !context.getUseDoublePrecision() && !context.getAccumulateInDouble())
+    if ((localDataSize/4)%2 == 0 && !context.getUseDoublePrecision() && !context.getUseMixedPrecision())
         defines["PARAMETER_SIZE_IS_EVEN"] = "1";
     if (context.getComputeCapability() >= 3.0 && !context.getUseDoublePrecision())
         defines["ENABLE_SHUFFLE"] = "1";
     stringstream defineAccum;
-    if (context.getAccumulateInDouble()) {
+    if (context.getUseMixedPrecision()) {
         defineAccum << "typedef double accum;\n";
         defineAccum << "typedef double3 accum3;\n";
         defines["make_accum3"] = "make_double3";

platforms/cuda/src/kernels/andersenThermostat.cu

@@ -2,12 +2,12 @@
  * Apply the Andersen thermostat to adjust particle velocities.
  */
 
-extern "C" __global__ void applyAndersenThermostat(float collisionFrequency, float kT, real4* velm, const real2* __restrict__ stepSize, const float4* __restrict__ random,
+extern "C" __global__ void applyAndersenThermostat(float collisionFrequency, float kT, mixed4* velm, const mixed4* __restrict__ stepSize, const float4* __restrict__ random,
         unsigned int randomIndex, const int* __restrict__ atomGroups) {
-    float collisionProbability = 1.0f-expf(-collisionFrequency*stepSize[0].y);
+    float collisionProbability = 1.0f-expf(-(float) (collisionFrequency*stepSize[0].y));
     float randomRange = erff(collisionProbability/sqrtf(2.0f));
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_ATOMS; index += blockDim.x*gridDim.x) {
-        real4 velocity = velm[index];
+        mixed4 velocity = velm[index];
         float4 selectRand = random[randomIndex+atomGroups[index]];
         float4 velRand = random[randomIndex+index];
         real scale = (selectRand.w > -randomRange && selectRand.w < randomRange ? 0 : 1);

platforms/cuda/src/kernels/integrationUtilities.cu

@@ -116,48 +116,73 @@ extern "C" __global__ void generateRandomNumbers(int numValues, float4* __restri
     seed[blockIdx.x*blockDim.x+threadIdx.x] = state;
 }
 
+/**
+ * Load the position of a particle.
+ */
+inline __device__ mixed4 loadPos(const real4* __restrict__ posq, const real4* __restrict__ posqCorrection, int index) {
+#ifdef USE_MIXED_PRECISION
+    real4 pos1 = posq[index];
+    real4 pos2 = posqCorrection[index];
+    return make_mixed4(pos1.x+(mixed)pos2.x, pos1.y+(mixed)pos2.y, pos1.z+(mixed)pos2.z, pos1.w);
+#else
+    return posq[index];
+#endif
+}
+
+/**
+ * Store the position of a particle.
+ */
+inline __device__ void storePos(real4* __restrict__ posq, real4* __restrict__ posqCorrection, int index, mixed4 pos) {
+#ifdef USE_MIXED_PRECISION
+    posq[index] = make_real4((real) pos.x, (real) pos.y, (real) pos.z, (real) pos.w);
+    posqCorrection[index] = make_real4(pos.x-(real) pos.x, pos.y-(real) pos.y, pos.z-(real) pos.z, 0);
+#else
+    posq[index] = pos;
+#endif
+}
+
 /**
  * Enforce constraints on SHAKE clusters
  */
-extern "C" __global__ void applyShakeToPositions(int numClusters, real tol, const real4* __restrict__ oldPos, real4* __restrict__ posDelta, const int4* __restrict__ clusterAtoms, const float4* __restrict__ clusterParams) {
+extern "C" __global__ void applyShakeToPositions(int numClusters, mixed tol, const real4* __restrict__ oldPos, real4* __restrict__ posCorrection, mixed4* __restrict__ posDelta, const int4* __restrict__ clusterAtoms, const float4* __restrict__ clusterParams) {
     int index = blockIdx.x*blockDim.x+threadIdx.x;
     while (index < numClusters) {
         // Load the data for this cluster.
 
         int4 atoms = clusterAtoms[index];
         float4 params = clusterParams[index];
-        real4 pos = oldPos[atoms.x];
-        real4 xpi = posDelta[atoms.x];
-        real4 pos1 = oldPos[atoms.y];
-        real4 xpj1 = posDelta[atoms.y];
-        real4 pos2 = make_real4(0);
-        real4 xpj2 = make_real4(0);
+        mixed4 pos = loadPos(oldPos, posCorrection, atoms.x);
+        mixed4 xpi = posDelta[atoms.x];
+        mixed4 pos1 = loadPos(oldPos, posCorrection, atoms.y);
+        mixed4 xpj1 = posDelta[atoms.y];
+        mixed4 pos2 = make_mixed4(0);
+        mixed4 xpj2 = make_mixed4(0);
         real invMassCentral = params.x;
         real avgMass = params.y;
         float d2 = params.z;
         float invMassPeripheral = params.w;
         if (atoms.z != -1) {
-            pos2 = oldPos[atoms.z];
+            pos2 = loadPos(oldPos, posCorrection, atoms.z);
             xpj2 = posDelta[atoms.z];
         }
-        real4 pos3 = make_real4(0);
-        real4 xpj3 = make_real4(0);
+        mixed4 pos3 = make_mixed4(0);
+        mixed4 xpj3 = make_mixed4(0);
         if (atoms.w != -1) {
-            pos3 = oldPos[atoms.w];
+            pos3 = loadPos(oldPos, posCorrection, atoms.w);
             xpj3 = posDelta[atoms.w];
         }
 
         // Precompute quantities.
 
-        real3 rij1 = make_real3(pos.x-pos1.x, pos.y-pos1.y, pos.z-pos1.z);
-        real3 rij2 = make_real3(pos.x-pos2.x, pos.y-pos2.y, pos.z-pos2.z);
-        real3 rij3 = make_real3(pos.x-pos3.x, pos.y-pos3.y, pos.z-pos3.z);
-        real rij1sq = rij1.x*rij1.x + rij1.y*rij1.y + rij1.z*rij1.z;
-        real rij2sq = rij2.x*rij2.x + rij2.y*rij2.y + rij2.z*rij2.z;
-        real rij3sq = rij3.x*rij3.x + rij3.y*rij3.y + rij3.z*rij3.z;
-        real ld1 = d2-rij1sq;
-        real ld2 = d2-rij2sq;
-        real ld3 = d2-rij3sq;
+        mixed3 rij1 = make_mixed3(pos.x-pos1.x, pos.y-pos1.y, pos.z-pos1.z);
+        mixed3 rij2 = make_mixed3(pos.x-pos2.x, pos.y-pos2.y, pos.z-pos2.z);
+        mixed3 rij3 = make_mixed3(pos.x-pos3.x, pos.y-pos3.y, pos.z-pos3.z);
+        mixed rij1sq = rij1.x*rij1.x + rij1.y*rij1.y + rij1.z*rij1.z;
+        mixed rij2sq = rij2.x*rij2.x + rij2.y*rij2.y + rij2.z*rij2.z;
+        mixed rij3sq = rij3.x*rij3.x + rij3.y*rij3.y + rij3.z*rij3.z;
+        mixed ld1 = d2-rij1sq;
+        mixed ld2 = d2-rij2sq;
+        mixed ld3 = d2-rij3sq;
 
         // Iterate until convergence.
 
@@ -165,13 +190,13 @@ extern "C" __global__ void applyShakeToPositions(int numClusters, real tol, cons
         int iteration = 0;
         while (iteration < 15 && !converged) {
             converged = true;
-            real3 rpij = make_real3(xpi.x-xpj1.x, xpi.y-xpj1.y, xpi.z-xpj1.z);
-            real rpsqij = rpij.x*rpij.x + rpij.y*rpij.y + rpij.z*rpij.z;
-            real rrpr = rij1.x*rpij.x + rij1.y*rpij.y + rij1.z*rpij.z;
-            real diff = fabs(ld1-2.0f*rrpr-rpsqij) / (d2*tol);
+            mixed3 rpij = make_mixed3(xpi.x-xpj1.x, xpi.y-xpj1.y, xpi.z-xpj1.z);
+            mixed rpsqij = rpij.x*rpij.x + rpij.y*rpij.y + rpij.z*rpij.z;
+            mixed rrpr = rij1.x*rpij.x + rij1.y*rpij.y + rij1.z*rpij.z;
+            mixed diff = fabs(ld1-2.0f*rrpr-rpsqij) / (d2*tol);
             if (diff >= 1.0f) {
-                real acor  = (ld1-2.0f*rrpr-rpsqij)*avgMass / (rrpr+rij1sq);
-                real3 dr = rij1*acor;
+                mixed acor  = (ld1-2.0f*rrpr-rpsqij)*avgMass / (rrpr+rij1sq);
+                mixed3 dr = rij1*acor;
                 xpi.x += dr.x*invMassCentral;
                 xpi.y += dr.y*invMassCentral;
                 xpi.z += dr.z*invMassCentral;
@@ -181,13 +206,13 @@ extern "C" __global__ void applyShakeToPositions(int numClusters, real tol, cons
                 converged = false;
             }
             if (atoms.z != -1) {
-                rpij = make_real3(xpi.x-xpj2.x, xpi.y-xpj2.y, xpi.z-xpj2.z);
+                rpij = make_mixed3(xpi.x-xpj2.x, xpi.y-xpj2.y, xpi.z-xpj2.z);
                 rpsqij = rpij.x*rpij.x + rpij.y*rpij.y + rpij.z*rpij.z;
                 rrpr = rij2.x*rpij.x + rij2.y*rpij.y + rij2.z*rpij.z;
                 diff = fabs(ld2-2.0f*rrpr-rpsqij) / (d2*tol);
                 if (diff >= 1.0f) {
-                    real acor  = (ld2 - 2.0f*rrpr - rpsqij)*avgMass / (rrpr + rij2sq);
-                    real3 dr = rij2*acor;
+                    mixed acor  = (ld2 - 2.0f*rrpr - rpsqij)*avgMass / (rrpr + rij2sq);
+                    mixed3 dr = rij2*acor;
                     xpi.x += dr.x*invMassCentral;
                     xpi.y += dr.y*invMassCentral;
                     xpi.z += dr.z*invMassCentral;
@@ -198,13 +223,13 @@ extern "C" __global__ void applyShakeToPositions(int numClusters, real tol, cons
                 }
             }
             if (atoms.w != -1) {
-                rpij = make_real3(xpi.x-xpj3.x, xpi.y-xpj3.y, xpi.z-xpj3.z);
+                rpij = make_mixed3(xpi.x-xpj3.x, xpi.y-xpj3.y, xpi.z-xpj3.z);
                 rpsqij = rpij.x*rpij.x + rpij.y*rpij.y + rpij.z*rpij.z;
                 rrpr = rij3.x*rpij.x + rij3.y*rpij.y + rij3.z*rpij.z;
                 diff = fabs(ld3 - 2.0f*rrpr - rpsqij) / (d2*tol);
                 if (diff >= 1.0f) {
-                    real acor  = (ld3-2.0f*rrpr-rpsqij)*avgMass / (rrpr+rij3sq);
-                    real3 dr = rij3*acor;
+                    mixed acor  = (ld3-2.0f*rrpr-rpsqij)*avgMass / (rrpr+rij3sq);
+                    mixed3 dr = rij3*acor;
                     xpi.x += dr.x*invMassCentral;
                     xpi.y += dr.y*invMassCentral;
                     xpi.z += dr.z*invMassCentral;
@@ -232,45 +257,45 @@ extern "C" __global__ void applyShakeToPositions(int numClusters, real tol, cons
 /**
  * Enforce velocity constraints on SHAKE clusters
  */
-extern "C" __global__ void applyShakeToVelocities(int numClusters, real tol, const real4* __restrict__ oldPos, real4* __restrict__ posDelta, const int4* __restrict__ clusterAtoms, const float4* __restrict__ clusterParams) {
+extern "C" __global__ void applyShakeToVelocities(int numClusters, mixed tol, const real4* __restrict__ oldPos, real4* __restrict__ posCorrection, mixed4* __restrict__ posDelta, const int4* __restrict__ clusterAtoms, const float4* __restrict__ clusterParams) {
     int index = blockIdx.x*blockDim.x+threadIdx.x;
     while (index < numClusters) {
         // Load the data for this cluster.
 
         int4 atoms = clusterAtoms[index];
         float4 params = clusterParams[index];
-        real4 pos = oldPos[atoms.x];
-        real4 xpi = posDelta[atoms.x];
-        real4 pos1 = oldPos[atoms.y];
-        real4 xpj1 = posDelta[atoms.y];
-        real4 pos2 = make_real4(0);
-        real4 xpj2 = make_real4(0);
+        mixed4 pos = loadPos(oldPos, posCorrection, atoms.x);
+        mixed4 xpi = posDelta[atoms.x];
+        mixed4 pos1 = loadPos(oldPos, posCorrection, atoms.y);
+        mixed4 xpj1 = posDelta[atoms.y];
+        mixed4 pos2 = make_mixed4(0);
+        mixed4 xpj2 = make_mixed4(0);
         real invMassCentral = params.x;
         real avgMass = params.y;
         float d2 = params.z;
         float invMassPeripheral = params.w;
         if (atoms.z != -1) {
-            pos2 = oldPos[atoms.z];
+            pos2 = loadPos(oldPos, posCorrection, atoms.z);
             xpj2 = posDelta[atoms.z];
         }
-        real4 pos3 = make_real4(0);
-        real4 xpj3 = make_real4(0);
+        mixed4 pos3 = make_mixed4(0);
+        mixed4 xpj3 = make_mixed4(0);
         if (atoms.w != -1) {
-            pos3 = oldPos[atoms.w];
+            pos3 = loadPos(oldPos, posCorrection, atoms.w);
             xpj3 = posDelta[atoms.w];
         }
 
         // Precompute quantities.
 
-        real3 rij1 = make_real3(pos.x-pos1.x, pos.y-pos1.y, pos.z-pos1.z);
-        real3 rij2 = make_real3(pos.x-pos2.x, pos.y-pos2.y, pos.z-pos2.z);
-        real3 rij3 = make_real3(pos.x-pos3.x, pos.y-pos3.y, pos.z-pos3.z);
-        real rij1sq = rij1.x*rij1.x + rij1.y*rij1.y + rij1.z*rij1.z;
-        real rij2sq = rij2.x*rij2.x + rij2.y*rij2.y + rij2.z*rij2.z;
-        real rij3sq = rij3.x*rij3.x + rij3.y*rij3.y + rij3.z*rij3.z;
-        real ld1 = d2-rij1sq;
-        real ld2 = d2-rij2sq;
-        real ld3 = d2-rij3sq;
+        mixed3 rij1 = make_mixed3(pos.x-pos1.x, pos.y-pos1.y, pos.z-pos1.z);
+        mixed3 rij2 = make_mixed3(pos.x-pos2.x, pos.y-pos2.y, pos.z-pos2.z);
+        mixed3 rij3 = make_mixed3(pos.x-pos3.x, pos.y-pos3.y, pos.z-pos3.z);
+        mixed rij1sq = rij1.x*rij1.x + rij1.y*rij1.y + rij1.z*rij1.z;
+        mixed rij2sq = rij2.x*rij2.x + rij2.y*rij2.y + rij2.z*rij2.z;
+        mixed rij3sq = rij3.x*rij3.x + rij3.y*rij3.y + rij3.z*rij3.z;
+        mixed ld1 = d2-rij1sq;
+        mixed ld2 = d2-rij2sq;
+        mixed ld3 = d2-rij3sq;
 
         // Iterate until convergence.
 
@@ -278,10 +303,10 @@ extern "C" __global__ void applyShakeToVelocities(int numClusters, real tol, con
         int iteration = 0;
         while (iteration < 15 && !converged) {
             converged = true;
-            real3 rpij = make_real3(xpi.x-xpj1.x, xpi.y-xpj1.y, xpi.z-xpj1.z);
-            real rrpr = rpij.x*rij1.x + rpij.y*rij1.y + rpij.z*rij1.z;
-            real delta = -2.0f*avgMass*rrpr/rij1sq;
-            real3 dr = rij1*delta;
+            mixed3 rpij = make_mixed3(xpi.x-xpj1.x, xpi.y-xpj1.y, xpi.z-xpj1.z);
+            mixed rrpr = rpij.x*rij1.x + rpij.y*rij1.y + rpij.z*rij1.z;
+            mixed delta = -2.0f*avgMass*rrpr/rij1sq;
+            mixed3 dr = rij1*delta;
             xpi.x += dr.x*invMassCentral;
             xpi.y += dr.y*invMassCentral;
             xpi.z += dr.z*invMassCentral;
@@ -291,7 +316,7 @@ extern "C" __global__ void applyShakeToVelocities(int numClusters, real tol, con
             if (fabs(delta) > tol)
                 converged = false;
             if (atoms.z != -1) {
-                rpij = make_real3(xpi.x-xpj2.x, xpi.y-xpj2.y, xpi.z-xpj2.z);
+                rpij = make_mixed3(xpi.x-xpj2.x, xpi.y-xpj2.y, xpi.z-xpj2.z);
                 rrpr = rpij.x*rij2.x + rpij.y*rij2.y + rpij.z*rij2.z;
                 delta = -2.0f*avgMass*rrpr/rij2sq;
                 dr = rij2*delta;
@@ -305,7 +330,7 @@ extern "C" __global__ void applyShakeToVelocities(int numClusters, real tol, con
                     converged = false;
             }
             if (atoms.w != -1) {
-                rpij = make_real3(xpi.x-xpj3.x, xpi.y-xpj3.y, xpi.z-xpj3.z);
+                rpij = make_mixed3(xpi.x-xpj3.x, xpi.y-xpj3.y, xpi.z-xpj3.z);
                 rrpr = rpij.x*rij3.x + rpij.y*rij3.y + rpij.z*rij3.z;
                 delta = -2.0f*avgMass*rrpr/rij3sq;
                 dr = rij3*delta;
@@ -336,135 +361,135 @@ extern "C" __global__ void applyShakeToVelocities(int numClusters, real tol, con
 /**
  * Enforce constraints on SETTLE clusters
  */
-extern "C" __global__ void applySettleToPositions(int numClusters, float tol, const real4* __restrict__ oldPos, real4* __restrict__ posDelta, const real4* __restrict__ velm, const int4* __restrict__ clusterAtoms, const float2* __restrict__ clusterParams) {
+extern "C" __global__ void applySettleToPositions(int numClusters, float tol, const real4* __restrict__ oldPos, real4* __restrict__ posCorrection, mixed4* __restrict__ posDelta, const mixed4* __restrict__ velm, const int4* __restrict__ clusterAtoms, const float2* __restrict__ clusterParams) {
     int index = blockIdx.x*blockDim.x+threadIdx.x;
     while (index < numClusters) {
         // Load the data for this cluster.
 
         int4 atoms = clusterAtoms[index];
         float2 params = clusterParams[index];
-        real4 apos0 = oldPos[atoms.x];
-        real4 xp0 = posDelta[atoms.x];
-        real4 apos1 = oldPos[atoms.y];
-        real4 xp1 = posDelta[atoms.y];
-        real4 apos2 = oldPos[atoms.z];
-        real4 xp2 = posDelta[atoms.z];
-        real m0 = RECIP(velm[atoms.x].w);
-        real m1 = RECIP(velm[atoms.y].w);
-        real m2 = RECIP(velm[atoms.z].w);
+        mixed4 apos0 = loadPos(oldPos, posCorrection, atoms.x);
+        mixed4 xp0 = posDelta[atoms.x];
+        mixed4 apos1 = loadPos(oldPos, posCorrection, atoms.y);
+        mixed4 xp1 = posDelta[atoms.y];
+        mixed4 apos2 = loadPos(oldPos, posCorrection, atoms.z);
+        mixed4 xp2 = posDelta[atoms.z];
+        mixed m0 = 1/velm[atoms.x].w;
+        mixed m1 = 1/velm[atoms.y].w;
+        mixed m2 = 1/velm[atoms.z].w;
 
         // Apply the SETTLE algorithm.
 
-        real xb0 = apos1.x-apos0.x;
-        real yb0 = apos1.y-apos0.y;
-        real zb0 = apos1.z-apos0.z;
-        real xc0 = apos2.x-apos0.x;
-        real yc0 = apos2.y-apos0.y;
-        real zc0 = apos2.z-apos0.z;
-
-        real invTotalMass = RECIP(m0+m1+m2);
-        real xcom = (xp0.x*m0 + (xb0+xp1.x)*m1 + (xc0+xp2.x)*m2) * invTotalMass;
-        real ycom = (xp0.y*m0 + (yb0+xp1.y)*m1 + (yc0+xp2.y)*m2) * invTotalMass;
-        real zcom = (xp0.z*m0 + (zb0+xp1.z)*m1 + (zc0+xp2.z)*m2) * invTotalMass;
-
-        real xa1 = xp0.x - xcom;
-        real ya1 = xp0.y - ycom;
-        real za1 = xp0.z - zcom;
-        real xb1 = xb0 + xp1.x - xcom;
-        real yb1 = yb0 + xp1.y - ycom;
-        real zb1 = zb0 + xp1.z - zcom;
-        real xc1 = xc0 + xp2.x - xcom;
-        real yc1 = yc0 + xp2.y - ycom;
-        real zc1 = zc0 + xp2.z - zcom;
-
-        real xaksZd = yb0*zc0 - zb0*yc0;
-        real yaksZd = zb0*xc0 - xb0*zc0;
-        real zaksZd = xb0*yc0 - yb0*xc0;
-        real xaksXd = ya1*zaksZd - za1*yaksZd;
-        real yaksXd = za1*xaksZd - xa1*zaksZd;
-        real zaksXd = xa1*yaksZd - ya1*xaksZd;
-        real xaksYd = yaksZd*zaksXd - zaksZd*yaksXd;
-        real yaksYd = zaksZd*xaksXd - xaksZd*zaksXd;
-        real zaksYd = xaksZd*yaksXd - yaksZd*xaksXd;
-
-        real axlng = SQRT(xaksXd*xaksXd + yaksXd*yaksXd + zaksXd*zaksXd);
-        real aylng = SQRT(xaksYd*xaksYd + yaksYd*yaksYd + zaksYd*zaksYd);
-        real azlng = SQRT(xaksZd*xaksZd + yaksZd*yaksZd + zaksZd*zaksZd);
-        real trns11 = xaksXd / axlng;
-        real trns21 = yaksXd / axlng;
-        real trns31 = zaksXd / axlng;
-        real trns12 = xaksYd / aylng;
-        real trns22 = yaksYd / aylng;
-        real trns32 = zaksYd / aylng;
-        real trns13 = xaksZd / azlng;
-        real trns23 = yaksZd / azlng;
-        real trns33 = zaksZd / azlng;
-
-        real xb0d = trns11*xb0 + trns21*yb0 + trns31*zb0;
-        real yb0d = trns12*xb0 + trns22*yb0 + trns32*zb0;
-        real xc0d = trns11*xc0 + trns21*yc0 + trns31*zc0;
-        real yc0d = trns12*xc0 + trns22*yc0 + trns32*zc0;
-        real za1d = trns13*xa1 + trns23*ya1 + trns33*za1;
-        real xb1d = trns11*xb1 + trns21*yb1 + trns31*zb1;
-        real yb1d = trns12*xb1 + trns22*yb1 + trns32*zb1;
-        real zb1d = trns13*xb1 + trns23*yb1 + trns33*zb1;
-        real xc1d = trns11*xc1 + trns21*yc1 + trns31*zc1;
-        real yc1d = trns12*xc1 + trns22*yc1 + trns32*zc1;
-        real zc1d = trns13*xc1 + trns23*yc1 + trns33*zc1;
+        mixed xb0 = apos1.x-apos0.x;
+        mixed yb0 = apos1.y-apos0.y;
+        mixed zb0 = apos1.z-apos0.z;
+        mixed xc0 = apos2.x-apos0.x;
+        mixed yc0 = apos2.y-apos0.y;
+        mixed zc0 = apos2.z-apos0.z;
+
+        mixed invTotalMass = 1/(m0+m1+m2);
+        mixed xcom = (xp0.x*m0 + (xb0+xp1.x)*m1 + (xc0+xp2.x)*m2) * invTotalMass;
+        mixed ycom = (xp0.y*m0 + (yb0+xp1.y)*m1 + (yc0+xp2.y)*m2) * invTotalMass;
+        mixed zcom = (xp0.z*m0 + (zb0+xp1.z)*m1 + (zc0+xp2.z)*m2) * invTotalMass;
+
+        mixed xa1 = xp0.x - xcom;
+        mixed ya1 = xp0.y - ycom;
+        mixed za1 = xp0.z - zcom;
+        mixed xb1 = xb0 + xp1.x - xcom;
+        mixed yb1 = yb0 + xp1.y - ycom;
+        mixed zb1 = zb0 + xp1.z - zcom;
+        mixed xc1 = xc0 + xp2.x - xcom;
+        mixed yc1 = yc0 + xp2.y - ycom;
+        mixed zc1 = zc0 + xp2.z - zcom;
+
+        mixed xaksZd = yb0*zc0 - zb0*yc0;
+        mixed yaksZd = zb0*xc0 - xb0*zc0;
+        mixed zaksZd = xb0*yc0 - yb0*xc0;
+        mixed xaksXd = ya1*zaksZd - za1*yaksZd;
+        mixed yaksXd = za1*xaksZd - xa1*zaksZd;
+        mixed zaksXd = xa1*yaksZd - ya1*xaksZd;
+        mixed xaksYd = yaksZd*zaksXd - zaksZd*yaksXd;
+        mixed yaksYd = zaksZd*xaksXd - xaksZd*zaksXd;
+        mixed zaksYd = xaksZd*yaksXd - yaksZd*xaksXd;
+
+        mixed axlng = sqrt(xaksXd*xaksXd + yaksXd*yaksXd + zaksXd*zaksXd);
+        mixed aylng = sqrt(xaksYd*xaksYd + yaksYd*yaksYd + zaksYd*zaksYd);
+        mixed azlng = sqrt(xaksZd*xaksZd + yaksZd*yaksZd + zaksZd*zaksZd);
+        mixed trns11 = xaksXd / axlng;
+        mixed trns21 = yaksXd / axlng;
+        mixed trns31 = zaksXd / axlng;
+        mixed trns12 = xaksYd / aylng;
+        mixed trns22 = yaksYd / aylng;
+        mixed trns32 = zaksYd / aylng;
+        mixed trns13 = xaksZd / azlng;
+        mixed trns23 = yaksZd / azlng;
+        mixed trns33 = zaksZd / azlng;
+
+        mixed xb0d = trns11*xb0 + trns21*yb0 + trns31*zb0;
+        mixed yb0d = trns12*xb0 + trns22*yb0 + trns32*zb0;
+        mixed xc0d = trns11*xc0 + trns21*yc0 + trns31*zc0;
+        mixed yc0d = trns12*xc0 + trns22*yc0 + trns32*zc0;
+        mixed za1d = trns13*xa1 + trns23*ya1 + trns33*za1;
+        mixed xb1d = trns11*xb1 + trns21*yb1 + trns31*zb1;
+        mixed yb1d = trns12*xb1 + trns22*yb1 + trns32*zb1;
+        mixed zb1d = trns13*xb1 + trns23*yb1 + trns33*zb1;
+        mixed xc1d = trns11*xc1 + trns21*yc1 + trns31*zc1;
+        mixed yc1d = trns12*xc1 + trns22*yc1 + trns32*zc1;
+        mixed zc1d = trns13*xc1 + trns23*yc1 + trns33*zc1;
 
         //                                        --- Step2  A2' ---
 
         float rc = 0.5f*params.y;
-        float rb = SQRT(params.x*params.x-rc*rc);
-        real ra = rb*(m1+m2)*invTotalMass;
+        float rb = sqrt(params.x*params.x-rc*rc);
+        mixed ra = rb*(m1+m2)*invTotalMass;
         rb -= ra;
-        real sinphi = za1d/ra;
-        real cosphi = SQRT(1-sinphi*sinphi);
-        real sinpsi = (zb1d-zc1d) / (2*rc*cosphi);
-        real cospsi = SQRT(1-sinpsi*sinpsi);
-
-        real ya2d =   ra*cosphi;
-        real xb2d = - rc*cospsi;
-        real yb2d = - rb*cosphi - rc*sinpsi*sinphi;
-        real yc2d = - rb*cosphi + rc*sinpsi*sinphi;
-        real xb2d2 = xb2d*xb2d;
-        real hh2 = 4.0f*xb2d2 + (yb2d-yc2d)*(yb2d-yc2d) + (zb1d-zc1d)*(zb1d-zc1d);
-        real deltx = 2.0f*xb2d + SQRT(4.0f*xb2d2 - hh2 + params.y*params.y);
+        mixed sinphi = za1d/ra;
+        mixed cosphi = sqrt(1-sinphi*sinphi);
+        mixed sinpsi = (zb1d-zc1d) / (2*rc*cosphi);
+        mixed cospsi = sqrt(1-sinpsi*sinpsi);
+
+        mixed ya2d =   ra*cosphi;
+        mixed xb2d = - rc*cospsi;
+        mixed yb2d = - rb*cosphi - rc*sinpsi*sinphi;
+        mixed yc2d = - rb*cosphi + rc*sinpsi*sinphi;
+        mixed xb2d2 = xb2d*xb2d;
+        mixed hh2 = 4.0f*xb2d2 + (yb2d-yc2d)*(yb2d-yc2d) + (zb1d-zc1d)*(zb1d-zc1d);
+        mixed deltx = 2.0f*xb2d + sqrt(4.0f*xb2d2 - hh2 + params.y*params.y);
         xb2d -= deltx*0.5f;
 
         //                                        --- Step3  al,be,ga ---
 
-        real alpha = (xb2d*(xb0d-xc0d) + yb0d*yb2d + yc0d*yc2d);
-        real beta = (xb2d*(yc0d-yb0d) + xb0d*yb2d + xc0d*yc2d);
-        real gamma = xb0d*yb1d - xb1d*yb0d + xc0d*yc1d - xc1d*yc0d;
+        mixed alpha = (xb2d*(xb0d-xc0d) + yb0d*yb2d + yc0d*yc2d);
+        mixed beta = (xb2d*(yc0d-yb0d) + xb0d*yb2d + xc0d*yc2d);
+        mixed gamma = xb0d*yb1d - xb1d*yb0d + xc0d*yc1d - xc1d*yc0d;
 
-        real al2be2 = alpha*alpha + beta*beta;
-        real sintheta = (alpha*gamma - beta*SQRT(al2be2 - gamma*gamma)) / al2be2;
+        mixed al2be2 = alpha*alpha + beta*beta;
+        mixed sintheta = (alpha*gamma - beta*sqrt(al2be2 - gamma*gamma)) / al2be2;
 
         //                                        --- Step4  A3' ---
 
-        real costheta = SQRT(1-sintheta*sintheta);
-        real xa3d = - ya2d*sintheta;
-        real ya3d =   ya2d*costheta;
-        real za3d = za1d;
-        real xb3d =   xb2d*costheta - yb2d*sintheta;
-        real yb3d =   xb2d*sintheta + yb2d*costheta;
-        real zb3d = zb1d;
-        real xc3d = - xb2d*costheta - yc2d*sintheta;
-        real yc3d = - xb2d*sintheta + yc2d*costheta;
-        real zc3d = zc1d;
+        mixed costheta = sqrt(1-sintheta*sintheta);
+        mixed xa3d = - ya2d*sintheta;
+        mixed ya3d =   ya2d*costheta;
+        mixed za3d = za1d;
+        mixed xb3d =   xb2d*costheta - yb2d*sintheta;
+        mixed yb3d =   xb2d*sintheta + yb2d*costheta;
+        mixed zb3d = zb1d;
+        mixed xc3d = - xb2d*costheta - yc2d*sintheta;
+        mixed yc3d = - xb2d*sintheta + yc2d*costheta;
+        mixed zc3d = zc1d;
 
         //                                        --- Step5  A3 ---
 
-        real xa3 = trns11*xa3d + trns12*ya3d + trns13*za3d;
-        real ya3 = trns21*xa3d + trns22*ya3d + trns23*za3d;
-        real za3 = trns31*xa3d + trns32*ya3d + trns33*za3d;
-        real xb3 = trns11*xb3d + trns12*yb3d + trns13*zb3d;
-        real yb3 = trns21*xb3d + trns22*yb3d + trns23*zb3d;
-        real zb3 = trns31*xb3d + trns32*yb3d + trns33*zb3d;
-        real xc3 = trns11*xc3d + trns12*yc3d + trns13*zc3d;
-        real yc3 = trns21*xc3d + trns22*yc3d + trns23*zc3d;
-        real zc3 = trns31*xc3d + trns32*yc3d + trns33*zc3d;
+        mixed xa3 = trns11*xa3d + trns12*ya3d + trns13*za3d;
+        mixed ya3 = trns21*xa3d + trns22*ya3d + trns23*za3d;
+        mixed za3 = trns31*xa3d + trns32*ya3d + trns33*za3d;
+        mixed xb3 = trns11*xb3d + trns12*yb3d + trns13*zb3d;
+        mixed yb3 = trns21*xb3d + trns22*yb3d + trns23*zb3d;
+        mixed zb3 = trns31*xb3d + trns32*yb3d + trns33*zb3d;
+        mixed xc3 = trns11*xc3d + trns12*yc3d + trns13*zc3d;
+        mixed yc3 = trns21*xc3d + trns22*yc3d + trns23*zc3d;
+        mixed zc3 = trns31*xc3d + trns32*yc3d + trns33*zc3d;
 
         xp0.x = xcom + xa3;
         xp0.y = ycom + ya3;
@@ -488,49 +513,49 @@ extern "C" __global__ void applySettleToPositions(int numClusters, float tol, co
 /**
  * Enforce velocity constraints on SETTLE clusters
  */
-extern "C" __global__ void applySettleToVelocities(int numClusters, float tol, const real4* __restrict__ oldPos, real4* __restrict__ posDelta, real4* __restrict__ velm, const int4* __restrict__ clusterAtoms, const float2* __restrict__ clusterParams) {
+extern "C" __global__ void applySettleToVelocities(int numClusters, float tol, const real4* __restrict__ oldPos, real4* __restrict__ posCorrection, mixed4* __restrict__ posDelta, mixed4* __restrict__ velm, const int4* __restrict__ clusterAtoms, const float2* __restrict__ clusterParams) {
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numClusters; index += blockDim.x*gridDim.x) {
         // Load the data for this cluster.
 
         int4 atoms = clusterAtoms[index];
-        real4 apos0 = oldPos[atoms.x];
-        real4 apos1 = oldPos[atoms.y];
-        real4 apos2 = oldPos[atoms.z];
-        real4 v0 = velm[atoms.x];
-        real4 v1 = velm[atoms.y];
-        real4 v2 = velm[atoms.z];
+        mixed4 apos0 = loadPos(oldPos, posCorrection, atoms.x);
+        mixed4 apos1 = loadPos(oldPos, posCorrection, atoms.y);
+        mixed4 apos2 = loadPos(oldPos, posCorrection, atoms.z);
+        mixed4 v0 = velm[atoms.x];
+        mixed4 v1 = velm[atoms.y];
+        mixed4 v2 = velm[atoms.z];
         
         // Compute intermediate quantities: the atom masses, the bond directions, the relative velocities,
         // and the angle cosines and sines.
         
-        real mA = RECIP(v0.w);
-        real mB = RECIP(v1.w);
-        real mC = RECIP(v2.w);
-        real3 eAB = make_real3(apos1.x-apos0.x, apos1.y-apos0.y, apos1.z-apos0.z);
-        real3 eBC = make_real3(apos2.x-apos1.x, apos2.y-apos1.y, apos2.z-apos1.z);
-        real3 eCA = make_real3(apos0.x-apos2.x, apos0.y-apos2.y, apos0.z-apos2.z);
-        eAB *= RECIP(SQRT(eAB.x*eAB.x + eAB.y*eAB.y + eAB.z*eAB.z));
-        eBC *= RECIP(SQRT(eBC.x*eBC.x + eBC.y*eBC.y + eBC.z*eBC.z));
-        eCA *= RECIP(SQRT(eCA.x*eCA.x + eCA.y*eCA.y + eCA.z*eCA.z));
-        real vAB = (v1.x-v0.x)*eAB.x + (v1.y-v0.y)*eAB.y + (v1.z-v0.z)*eAB.z;
-        real vBC = (v2.x-v1.x)*eBC.x + (v2.y-v1.y)*eBC.y + (v2.z-v1.z)*eBC.z;
-        real vCA = (v0.x-v2.x)*eCA.x + (v0.y-v2.y)*eCA.y + (v0.z-v2.z)*eCA.z;
-        real cA = -(eAB.x*eCA.x + eAB.y*eCA.y + eAB.z*eCA.z);
-        real cB = -(eAB.x*eBC.x + eAB.y*eBC.y + eAB.z*eBC.z);
-        real cC = -(eBC.x*eCA.x + eBC.y*eCA.y + eBC.z*eCA.z);
-        real s2A = 1-cA*cA;
-        real s2B = 1-cB*cB;
-        real s2C = 1-cC*cC;
+        mixed mA = 1/v0.w;
+        mixed mB = 1/v1.w;
+        mixed mC = 1/v2.w;
+        mixed3 eAB = make_mixed3(apos1.x-apos0.x, apos1.y-apos0.y, apos1.z-apos0.z);
+        mixed3 eBC = make_mixed3(apos2.x-apos1.x, apos2.y-apos1.y, apos2.z-apos1.z);
+        mixed3 eCA = make_mixed3(apos0.x-apos2.x, apos0.y-apos2.y, apos0.z-apos2.z);
+        eAB *= rsqrt(eAB.x*eAB.x + eAB.y*eAB.y + eAB.z*eAB.z);
+        eBC *= rsqrt(eBC.x*eBC.x + eBC.y*eBC.y + eBC.z*eBC.z);
+        eCA *= rsqrt(eCA.x*eCA.x + eCA.y*eCA.y + eCA.z*eCA.z);
+        mixed vAB = (v1.x-v0.x)*eAB.x + (v1.y-v0.y)*eAB.y + (v1.z-v0.z)*eAB.z;
+        mixed vBC = (v2.x-v1.x)*eBC.x + (v2.y-v1.y)*eBC.y + (v2.z-v1.z)*eBC.z;
+        mixed vCA = (v0.x-v2.x)*eCA.x + (v0.y-v2.y)*eCA.y + (v0.z-v2.z)*eCA.z;
+        mixed cA = -(eAB.x*eCA.x + eAB.y*eCA.y + eAB.z*eCA.z);
+        mixed cB = -(eAB.x*eBC.x + eAB.y*eBC.y + eAB.z*eBC.z);
+        mixed cC = -(eBC.x*eCA.x + eBC.y*eCA.y + eBC.z*eCA.z);
+        mixed s2A = 1-cA*cA;
+        mixed s2B = 1-cB*cB;
+        mixed s2C = 1-cC*cC;
         
         // Solve the equations.  These are different from those in the SETTLE paper (JCC 13(8), pp. 952-962, 1992), because
         // in going from equations B1 to B2, they make the assumption that mB=mC (but don't bother to mention they're
         // making that assumption).  We allow all three atoms to have different masses.
         
-        real mABCinv = RECIP(mA*mB*mC);
-        real denom = (((s2A*mB+s2B*mA)*mC+(s2A*mB*mB+2*(cA*cB*cC+1)*mA*mB+s2B*mA*mA))*mC+s2C*mA*mB*(mA+mB))*mABCinv;
-        real tab = ((cB*cC*mA-cA*mB-cA*mC)*vCA + (cA*cC*mB-cB*mC-cB*mA)*vBC + (s2C*mA*mA*mB*mB*mABCinv+(mA+mB+mC))*vAB)/denom;
-        real tbc = ((cA*cB*mC-cC*mB-cC*mA)*vCA + (s2A*mB*mB*mC*mC*mABCinv+(mA+mB+mC))*vBC + (cA*cC*mB-cB*mA-cB*mC)*vAB)/denom;
-        real tca = ((s2B*mA*mA*mC*mC*mABCinv+(mA+mB+mC))*vCA + (cA*cB*mC-cC*mB-cC*mA)*vBC + (cB*cC*mA-cA*mB-cA*mC)*vAB)/denom;
+        mixed mABCinv = 1/(mA*mB*mC);
+        mixed denom = (((s2A*mB+s2B*mA)*mC+(s2A*mB*mB+2*(cA*cB*cC+1)*mA*mB+s2B*mA*mA))*mC+s2C*mA*mB*(mA+mB))*mABCinv;
+        mixed tab = ((cB*cC*mA-cA*mB-cA*mC)*vCA + (cA*cC*mB-cB*mC-cB*mA)*vBC + (s2C*mA*mA*mB*mB*mABCinv+(mA+mB+mC))*vAB)/denom;
+        mixed tbc = ((cA*cB*mC-cC*mB-cC*mA)*vCA + (s2A*mB*mB*mC*mC*mABCinv+(mA+mB+mC))*vBC + (cA*cC*mB-cB*mA-cB*mC)*vAB)/denom;
+        mixed tca = ((s2B*mA*mA*mC*mC*mABCinv+(mA+mB+mC))*vCA + (cA*cB*mC-cC*mB-cC*mA)*vBC + (cB*cC*mA-cA*mB-cA*mC)*vAB)/denom;
         v0.x += (tab*eAB.x - tca*eCA.x)*v0.w;
         v0.y += (tab*eAB.y - tca*eCA.y)*v0.w;
         v0.z += (tab*eAB.z - tca*eCA.z)*v0.w;
@@ -549,14 +574,14 @@ extern "C" __global__ void applySettleToVelocities(int numClusters, float tol, c
 /**
  * Compute the direction each CCMA constraint is pointing in.  This is called once at the beginning of constraint evaluation.
  */
-extern "C" __global__ void computeCCMAConstraintDirections(const int2* __restrict__ constraintAtoms, real4* __restrict__ constraintDistance, const real4* __restrict__ atomPositions) {
+extern "C" __global__ void computeCCMAConstraintDirections(const int2* __restrict__ constraintAtoms, mixed4* __restrict__ constraintDistance, const real4* __restrict__ atomPositions, const real4* __restrict__ posqCorrection) {
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_CCMA_CONSTRAINTS; index += blockDim.x*gridDim.x) {
         // Compute the direction for this constraint.
 
         int2 atoms = constraintAtoms[index];
-        real4 dir = constraintDistance[index];
-        real4 oldPos1 = atomPositions[atoms.x];
-        real4 oldPos2 = atomPositions[atoms.y];
+        mixed4 dir = constraintDistance[index];
+        mixed4 oldPos1 = loadPos(atomPositions, posqCorrection, atoms.x);
+        mixed4 oldPos2 = loadPos(atomPositions, posqCorrection, atoms.y);
         dir.x = oldPos1.x-oldPos2.x;
         dir.y = oldPos1.y-oldPos2.y;
         dir.z = oldPos1.z-oldPos2.z;
@@ -567,8 +592,8 @@ extern "C" __global__ void computeCCMAConstraintDirections(const int2* __restric
 /**
  * Compute the force applied by each CCMA position constraint.
  */
-extern "C" __global__ void computeCCMAPositionConstraintForce(const int2* __restrict__ constraintAtoms, const real4* __restrict__ constraintDistance, const real4* __restrict__ atomPositions,
-        const real* __restrict__ reducedMass, real* __restrict__ delta1, int* __restrict__ converged, float tol, int iteration) {
+extern "C" __global__ void computeCCMAPositionConstraintForce(const int2* __restrict__ constraintAtoms, const mixed4* __restrict__ constraintDistance, const mixed4* __restrict__ atomPositions,
+        const mixed* __restrict__ reducedMass, mixed* __restrict__ delta1, int* __restrict__ converged, mixed tol, int iteration) {
     __shared__ int groupConverged;
     if (converged[1-iteration%2]) {
         if (blockIdx.x == 0 && threadIdx.x == 0)
@@ -578,22 +603,22 @@ extern "C" __global__ void computeCCMAPositionConstraintForce(const int2* __rest
     if (threadIdx.x == 0)
         groupConverged = 1;
     __syncthreads();
-    real lowerTol = 1.0f-2.0f*tol+tol*tol;
-    real upperTol = 1.0f+2.0f*tol+tol*tol;
+    mixed lowerTol = 1.0f-2.0f*tol+tol*tol;
+    mixed upperTol = 1.0f+2.0f*tol+tol*tol;
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_CCMA_CONSTRAINTS; index += blockDim.x*gridDim.x) {
         // Compute the force due to this constraint.
 
         int2 atoms = constraintAtoms[index];
-        real4 dir = constraintDistance[index];
-        real4 rp_ij = atomPositions[atoms.x]-atomPositions[atoms.y];
+        mixed4 dir = constraintDistance[index];
+        mixed4 rp_ij = atomPositions[atoms.x]-atomPositions[atoms.y];
         rp_ij.x += dir.x;
         rp_ij.y += dir.y;
         rp_ij.z += dir.z;
-        real rrpr = rp_ij.x*dir.x + rp_ij.y*dir.y + rp_ij.z*dir.z;
-        real d_ij2 = dir.x*dir.x + dir.y*dir.y + dir.z*dir.z;
-        real rp2 = rp_ij.x*rp_ij.x + rp_ij.y*rp_ij.y + rp_ij.z*rp_ij.z;
-        real dist2 = dir.w*dir.w;
-        real diff = dist2 - rp2;
+        mixed rrpr = rp_ij.x*dir.x + rp_ij.y*dir.y + rp_ij.z*dir.z;
+        mixed d_ij2 = dir.x*dir.x + dir.y*dir.y + dir.z*dir.z;
+        mixed rp2 = rp_ij.x*rp_ij.x + rp_ij.y*rp_ij.y + rp_ij.z*rp_ij.z;
+        mixed dist2 = dir.w*dir.w;
+        mixed diff = dist2 - rp2;
         delta1[index] = (rrpr > d_ij2*1e-6f ? reducedMass[index]*diff/rrpr : 0.0f);
 
         // See whether it has converged.
@@ -608,8 +633,8 @@ extern "C" __global__ void computeCCMAPositionConstraintForce(const int2* __rest
 /**
  * Compute the force applied by each CCMA velocity constraint.
  */
-extern "C" __global__ void computeCCMAVelocityConstraintForce(const int2* __restrict__ constraintAtoms, const real4* __restrict__ constraintDistance, const real4* __restrict__ atomPositions,
-        const real* __restrict__ reducedMass, real* __restrict__ delta1, int* __restrict__ converged, float tol, int iteration) {
+extern "C" __global__ void computeCCMAVelocityConstraintForce(const int2* __restrict__ constraintAtoms, const mixed4* __restrict__ constraintDistance, const mixed4* __restrict__ atomPositions,
+        const mixed* __restrict__ reducedMass, mixed* __restrict__ delta1, int* __restrict__ converged, float tol, int iteration) {
     __shared__ int groupConverged;
     if (converged[1-iteration%2]) {
         if (blockIdx.x == 0 && threadIdx.x == 0)
@@ -619,16 +644,16 @@ extern "C" __global__ void computeCCMAVelocityConstraintForce(const int2* __rest
     if (threadIdx.x == 0)
         groupConverged = 1;
     __syncthreads();
-    real lowerTol = 1.0f-2.0f*tol+tol*tol;
-    real upperTol = 1.0f+2.0f*tol+tol*tol;
+    mixed lowerTol = 1.0f-2.0f*tol+tol*tol;
+    mixed upperTol = 1.0f+2.0f*tol+tol*tol;
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_CCMA_CONSTRAINTS; index += blockDim.x*gridDim.x) {
         // Compute the force due to this constraint.
 
         int2 atoms = constraintAtoms[index];
-        real4 dir = constraintDistance[index];
-        real4 rp_ij = atomPositions[atoms.x]-atomPositions[atoms.y];
-        real rrpr = rp_ij.x*dir.x + rp_ij.y*dir.y + rp_ij.z*dir.z;
-        real d_ij2 = dir.x*dir.x + dir.y*dir.y + dir.z*dir.z;
+        mixed4 dir = constraintDistance[index];
+        mixed4 rp_ij = atomPositions[atoms.x]-atomPositions[atoms.y];
+        mixed rrpr = rp_ij.x*dir.x + rp_ij.y*dir.y + rp_ij.z*dir.z;
+        mixed d_ij2 = dir.x*dir.x + dir.y*dir.y + dir.z*dir.z;
         delta1[index] = -2.0f*reducedMass[index]*rrpr/d_ij2;
 
         // See whether it has converged.
@@ -643,15 +668,15 @@ extern "C" __global__ void computeCCMAVelocityConstraintForce(const int2* __rest
 /**
  * Multiply the vector of CCMA constraint forces by the constraint matrix.
  */
-extern "C" __global__ void multiplyByCCMAConstraintMatrix(const real* __restrict__ delta1, real* __restrict__ delta2, const int* __restrict__ constraintMatrixColumn,
-        const real* __restrict__ constraintMatrixValue, const int* __restrict__ converged, int iteration) {
+extern "C" __global__ void multiplyByCCMAConstraintMatrix(const mixed* __restrict__ delta1, mixed* __restrict__ delta2, const int* __restrict__ constraintMatrixColumn,
+        const mixed* __restrict__ constraintMatrixValue, const int* __restrict__ converged, int iteration) {
     if (converged[iteration%2])
         return; // The constraint iteration has already converged.
 
     // Multiply by the inverse constraint matrix.
 
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_CCMA_CONSTRAINTS; index += blockDim.x*gridDim.x) {
-        real sum = 0.0f;
+        mixed sum = 0;
         for (int i = 0; ; i++) {
             int element = index+i*NUM_CCMA_CONSTRAINTS;
             int column = constraintMatrixColumn[element];
@@ -666,26 +691,26 @@ extern "C" __global__ void multiplyByCCMAConstraintMatrix(const real* __restrict
 /**
  * Update the atom positions based on CCMA constraint forces.
  */
-extern "C" __global__ void updateCCMAAtomPositions(const int* __restrict__ numAtomConstraints, const int* __restrict__ atomConstraints, const real4* __restrict__ constraintDistance,
-        real4* __restrict__ atomPositions, const real4* __restrict__ velm, const real* __restrict__ delta1, const real* __restrict__ delta2, int* __restrict__ converged, int iteration) {
+extern "C" __global__ void updateCCMAAtomPositions(const int* __restrict__ numAtomConstraints, const int* __restrict__ atomConstraints, const mixed4* __restrict__ constraintDistance,
+        mixed4* __restrict__ atomPositions, const mixed4* __restrict__ velm, const mixed* __restrict__ delta1, const mixed* __restrict__ delta2, int* __restrict__ converged, int iteration) {
     if (blockIdx.x == 0 && threadIdx.x == 0)
         converged[1-iteration%2] = 1;
     if (converged[iteration%2])
         return; // The constraint iteration has already converged.
-    real damping = (iteration < 2 ? 0.5f : 1.0f);
+    mixed damping = (iteration < 2 ? 0.5f : 1.0f);
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_ATOMS; index += blockDim.x*gridDim.x) {
         // Compute the new position of this atom.
 
-        real4 atomPos = atomPositions[index];
-        real invMass = velm[index].w;
+        mixed4 atomPos = atomPositions[index];
+        mixed invMass = velm[index].w;
         int num = numAtomConstraints[index];
         for (int i = 0; i < num; i++) {
             int constraint = atomConstraints[index+i*NUM_ATOMS];
             bool forward = (constraint > 0);
             constraint = (forward ? constraint-1 : -constraint-1);
-            real constraintForce = damping*invMass*delta2[constraint];
+            mixed constraintForce = damping*invMass*delta2[constraint];
             constraintForce = (forward ? constraintForce : -constraintForce);
-            real4 dir = constraintDistance[constraint];
+            mixed4 dir = constraintDistance[constraint];
             atomPos.x += constraintForce*dir.x;
             atomPos.y += constraintForce*dir.y;
             atomPos.z += constraintForce*dir.z;
@@ -697,7 +722,7 @@ extern "C" __global__ void updateCCMAAtomPositions(const int* __restrict__ numAt
 /**
  * Compute the positions of virtual sites
  */
-extern "C" __global__ void computeVirtualSites(real4* __restrict__ posq, const int4* __restrict__ avg2Atoms, const real2* __restrict__ avg2Weights,
+extern "C" __global__ void computeVirtualSites(real4* __restrict__ posq, real4* __restrict__ posqCorrection, const int4* __restrict__ avg2Atoms, const real2* __restrict__ avg2Weights,
         const int4* __restrict__ avg3Atoms, const real4* __restrict__ avg3Weights,
         const int4* __restrict__ outOfPlaneAtoms, const real4* __restrict__ outOfPlaneWeights) {
     
@@ -706,13 +731,13 @@ extern "C" __global__ void computeVirtualSites(real4* __restrict__ posq, const i
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_2_AVERAGE; index += blockDim.x*gridDim.x) {
         int4 atoms = avg2Atoms[index];
         real2 weights = avg2Weights[index];
-        real4 pos = posq[atoms.x];
-        real4 pos1 = posq[atoms.y];
-        real4 pos2 = posq[atoms.z];
+        mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
+        mixed4 pos1 = loadPos(posq, posqCorrection, atoms.y);
+        mixed4 pos2 = loadPos(posq, posqCorrection, atoms.z);
         pos.x = pos1.x*weights.x + pos2.x*weights.y;
         pos.y = pos1.y*weights.x + pos2.y*weights.y;
         pos.z = pos1.z*weights.x + pos2.z*weights.y;
-        posq[atoms.x] = pos;
+        storePos(posq, posqCorrection, atoms.x, pos);
     }
     
     // Three particle average sites.
@@ -720,14 +745,14 @@ extern "C" __global__ void computeVirtualSites(real4* __restrict__ posq, const i
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_3_AVERAGE; index += blockDim.x*gridDim.x) {
         int4 atoms = avg3Atoms[index];
         real4 weights = avg3Weights[index];
-        real4 pos = posq[atoms.x];
-        real4 pos1 = posq[atoms.y];
-        real4 pos2 = posq[atoms.z];
-        real4 pos3 = posq[atoms.w];
+        mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
+        mixed4 pos1 = loadPos(posq, posqCorrection, atoms.y);
+        mixed4 pos2 = loadPos(posq, posqCorrection, atoms.z);
+        mixed4 pos3 = loadPos(posq, posqCorrection, atoms.w);
         pos.x = pos1.x*weights.x + pos2.x*weights.y + pos3.x*weights.z;
         pos.y = pos1.y*weights.x + pos2.y*weights.y + pos3.y*weights.z;
         pos.z = pos1.z*weights.x + pos2.z*weights.y + pos3.z*weights.z;
-        posq[atoms.x] = pos;
+        storePos(posq, posqCorrection, atoms.x, pos);
     }
     
     // Out of plane sites.
@@ -735,22 +760,22 @@ extern "C" __global__ void computeVirtualSites(real4* __restrict__ posq, const i
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_OUT_OF_PLANE; index += blockDim.x*gridDim.x) {
         int4 atoms = outOfPlaneAtoms[index];
         real4 weights = outOfPlaneWeights[index];
-        real4 pos = posq[atoms.x];
-        real4 pos1 = posq[atoms.y];
-        real4 pos2 = posq[atoms.z];
-        real4 pos3 = posq[atoms.w];
-        real4 v12 = pos2-pos1;
-        real4 v13 = pos3-pos1;
-        real3 cr = cross(v12, v13);
+        mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
+        mixed4 pos1 = loadPos(posq, posqCorrection, atoms.y);
+        mixed4 pos2 = loadPos(posq, posqCorrection, atoms.z);
+        mixed4 pos3 = loadPos(posq, posqCorrection, atoms.w);
+        mixed4 v12 = pos2-pos1;
+        mixed4 v13 = pos3-pos1;
+        mixed3 cr = cross(v12, v13);
         pos.x = pos1.x + v12.x*weights.x + v13.x*weights.y + cr.x*weights.z;
         pos.y = pos1.y + v12.y*weights.x + v13.y*weights.y + cr.y*weights.z;
         pos.z = pos1.z + v12.z*weights.x + v13.z*weights.y + cr.z*weights.z;
-        posq[atoms.x] = pos;
+        storePos(posq, posqCorrection, atoms.x, pos);
     }
 }
 
 inline __device__ real3 loadForce(int index, long long* __restrict__ force) {
-    real scale = RECIP((real) 0xFFFFFFFF);
+    real scale = 1/((real) 0xFFFFFFFF);
     return make_real3(scale*force[index], scale*force[index+PADDED_NUM_ATOMS], scale*force[index+PADDED_NUM_ATOMS*2]);
 }
 
@@ -763,7 +788,7 @@ inline __device__ void storeForce(int index, long long* __restrict__ force, real
 /**
  * Distribute forces from virtual sites to the atoms they are based on.
  */
-extern "C" __global__ void distributeVirtualSiteForces(const real4* __restrict__ posq, long long* __restrict__ force,
+extern "C" __global__ void distributeVirtualSiteForces(const real4* __restrict__ posq, const real4* __restrict__ posqCorrection, long long* __restrict__ force,
         const int4* __restrict__ avg2Atoms, const real2* __restrict__ avg2Weights,
         const int4* __restrict__ avg3Atoms, const real4* __restrict__ avg3Weights,
         const int4* __restrict__ outOfPlaneAtoms, const real4* __restrict__ outOfPlaneWeights) {
@@ -804,21 +829,21 @@ extern "C" __global__ void distributeVirtualSiteForces(const real4* __restrict__
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_OUT_OF_PLANE; index += blockDim.x*gridDim.x) {
         int4 atoms = outOfPlaneAtoms[index];
         real4 weights = outOfPlaneWeights[index];
-        real4 pos1 = posq[atoms.y];
-        real4 pos2 = posq[atoms.z];
-        real4 pos3 = posq[atoms.w];
-        real4 v12 = pos2-pos1;
-        real4 v13 = pos3-pos1;
+        mixed4 pos1 = loadPos(posq, posqCorrection, atoms.y);
+        mixed4 pos2 = loadPos(posq, posqCorrection, atoms.z);
+        mixed4 pos3 = loadPos(posq, posqCorrection, atoms.w);
+        mixed4 v12 = pos2-pos1;
+        mixed4 v13 = pos3-pos1;
         real3 f = loadForce(atoms.x, force);
         real3 f1 = loadForce(atoms.y, force);
         real3 f2 = loadForce(atoms.z, force);
         real3 f3 = loadForce(atoms.w, force);
-        real3 fp2 = make_real3(weights.x*f.x - weights.z*v13.z*f.y + weights.z*v13.y*f.z,
-                   weights.z*v13.z*f.x + weights.x*f.y - weights.z*v13.x*f.z,
-                  -weights.z*v13.y*f.x + weights.z*v13.x*f.y + weights.x*f.z);
-        real3 fp3 = make_real3(weights.y*f.x + weights.z*v12.z*f.y - weights.z*v12.y*f.z,
-                  -weights.z*v12.z*f.x + weights.y*f.y + weights.z*v12.x*f.z,
-                   weights.z*v12.y*f.x - weights.z*v12.x*f.y + weights.y*f.z);
+        real3 fp2 = make_real3((real) (weights.x*f.x - weights.z*v13.z*f.y + weights.z*v13.y*f.z),
+                   (real) (weights.z*v13.z*f.x + weights.x*f.y - weights.z*v13.x*f.z),
+                   (real) (-weights.z*v13.y*f.x + weights.z*v13.x*f.y + weights.x*f.z));
+        real3 fp3 = make_real3((real) (weights.y*f.x + weights.z*v12.z*f.y - weights.z*v12.y*f.z),
+                   (real) (-weights.z*v12.z*f.x + weights.y*f.y + weights.z*v12.x*f.z),
+                   (real) (weights.z*v12.y*f.x - weights.z*v12.x*f.y + weights.y*f.z));
         f1 += f-fp2-fp3;
         f2 += fp2;
         f3 += fp3;

platforms/cuda/src/kernels/langevin.cu

@@ -4,23 +4,23 @@ enum {VelScale, ForceScale, NoiseScale, MaxParams};
  * Perform the first step of Langevin integration.
  */
 
-extern "C" __global__ void integrateLangevinPart1(real4* __restrict__ velm, const long long* __restrict__ force, real4* __restrict__ posDelta,
-        const real* __restrict__ paramBuffer, const real2* __restrict__ dt, const float4* __restrict__ random, unsigned int randomIndex) {
-    real vscale = paramBuffer[VelScale];
-    real fscale = paramBuffer[ForceScale]/(real) 0xFFFFFFFF;
-    real noisescale = paramBuffer[NoiseScale];
-    real stepSize = dt[0].y;
+extern "C" __global__ void integrateLangevinPart1(mixed4* __restrict__ velm, const long long* __restrict__ force, mixed4* __restrict__ posDelta,
+        const mixed* __restrict__ paramBuffer, const mixed2* __restrict__ dt, const float4* __restrict__ random, unsigned int randomIndex) {
+    mixed vscale = paramBuffer[VelScale];
+    mixed fscale = paramBuffer[ForceScale]/(mixed) 0xFFFFFFFF;
+    mixed noisescale = paramBuffer[NoiseScale];
+    mixed stepSize = dt[0].y;
     int index = blockIdx.x*blockDim.x+threadIdx.x;
     randomIndex += index;
     while (index < NUM_ATOMS) {
-        real4 velocity = velm[index];
+        mixed4 velocity = velm[index];
         if (velocity.w != 0) {
-            real sqrtInvMass = SQRT(velocity.w);
+            mixed sqrtInvMass = SQRT(velocity.w);
             velocity.x = vscale*velocity.x + fscale*velocity.w*force[index] + noisescale*sqrtInvMass*random[randomIndex].x;
             velocity.y = vscale*velocity.y + fscale*velocity.w*force[index+PADDED_NUM_ATOMS] + noisescale*sqrtInvMass*random[randomIndex].y;
             velocity.z = vscale*velocity.z + fscale*velocity.w*force[index+PADDED_NUM_ATOMS*2] + noisescale*sqrtInvMass*random[randomIndex].z;
             velm[index] = velocity;
-            posDelta[index] = make_real4(stepSize*velocity.x, stepSize*velocity.y, stepSize*velocity.z, 0);
+            posDelta[index] = make_mixed4(stepSize*velocity.x, stepSize*velocity.y, stepSize*velocity.z, 0);
         }
         randomIndex += blockDim.x*gridDim.x;
         index += blockDim.x*gridDim.x;
@@ -31,21 +31,32 @@ extern "C" __global__ void integrateLangevinPart1(real4* __restrict__ velm, cons
  * Perform the second step of Langevin integration.
  */
 
-extern "C" __global__ void integrateLangevinPart2(real4* __restrict__ posq, const real4* __restrict__ posDelta, real4* __restrict__ velm, const real2* __restrict__ dt) {
+extern "C" __global__ void integrateLangevinPart2(real4* __restrict__ posq, real4* __restrict__ posqCorrection, const mixed4* __restrict__ posDelta, mixed4* __restrict__ velm, const mixed2* __restrict__ dt) {
     double invStepSize = 1.0/dt[0].y;
     int index = blockIdx.x*blockDim.x+threadIdx.x;
     while (index < NUM_ATOMS) {
-        real4 vel = velm[index];
+        mixed4 vel = velm[index];
         if (vel.w != 0) {
+#ifdef USE_MIXED_PRECISION
+            real4 pos1 = posq[index];
+            real4 pos2 = posqCorrection[index];
+            mixed4 pos = make_mixed4(pos1.x+(mixed)pos2.x, pos1.y+(mixed)pos2.y, pos1.z+(mixed)pos2.z, pos1.w);
+#else
             real4 pos = posq[index];
-            real4 delta = posDelta[index];
+#endif
+            mixed4 delta = posDelta[index];
             pos.x += delta.x;
             pos.y += delta.y;
             pos.z += delta.z;
-            vel.x = (real) invStepSize*delta.x;
-            vel.y = (real) invStepSize*delta.y;
-            vel.z = (real) invStepSize*delta.z;
+            vel.x = (mixed) (invStepSize*delta.x);
+            vel.y = (mixed) (invStepSize*delta.y);
+            vel.z = (mixed) (invStepSize*delta.z);
+#ifdef USE_MIXED_PRECISION
+            posq[index] = make_real4((real) pos.x, (real) pos.y, (real) pos.z, (real) pos.w);
+            posqCorrection[index] = make_real4(pos.x-(real) pos.x, pos.y-(real) pos.y, pos.z-(real) pos.z, 0);
+#else
             posq[index] = pos;
+#endif
             velm[index] = vel;
         }
         index += blockDim.x*gridDim.x;
@@ -56,18 +67,18 @@ extern "C" __global__ void integrateLangevinPart2(real4* __restrict__ posq, cons
  * Select the step size to use for the next step.
  */
 
-extern "C" __global__ void selectLangevinStepSize(real maxStepSize, real errorTol, real tau, real kT, real2* __restrict__ dt,
-        const real4* __restrict__ velm, const long long* __restrict__ force, real* __restrict__ paramBuffer) {
+extern "C" __global__ void selectLangevinStepSize(mixed maxStepSize, mixed errorTol, mixed tau, mixed kT, mixed2* __restrict__ dt,
+        const mixed4* __restrict__ velm, const long long* __restrict__ force, mixed* __restrict__ paramBuffer) {
     // Calculate the error.
 
-    extern __shared__ real params[];
-    real* error = &params[MaxParams];
-    real err = 0;
+    extern __shared__ mixed params[];
+    mixed* error = &params[MaxParams];
+    mixed err = 0;
     unsigned int index = threadIdx.x;
-    const real scale = RECIP((real) 0xFFFFFFFF);
+    const mixed scale = RECIP((mixed) 0xFFFFFFFF);
     while (index < NUM_ATOMS) {
-        real3 f = make_real3(scale*force[index], scale*force[index+PADDED_NUM_ATOMS], scale*force[index+PADDED_NUM_ATOMS*2]);
-        real invMass = velm[index].w;
+        mixed3 f = make_mixed3(scale*force[index], scale*force[index+PADDED_NUM_ATOMS], scale*force[index+PADDED_NUM_ATOMS*2]);
+        mixed invMass = velm[index].w;
         err += (f.x*f.x + f.y*f.y + f.z*f.z)*invMass;
         index += blockDim.x*gridDim.x;
     }
@@ -84,9 +95,9 @@ extern "C" __global__ void selectLangevinStepSize(real maxStepSize, real errorTo
     if (blockIdx.x*blockDim.x+threadIdx.x == 0) {
         // Select the new step size.
 
-        real totalError = sqrt(error[0]/(NUM_ATOMS*3));
-        real newStepSize = sqrt(errorTol/totalError);
-        real oldStepSize = dt[0].y;
+        mixed totalError = sqrt(error[0]/(NUM_ATOMS*3));
+        mixed newStepSize = sqrt(errorTol/totalError);
+        mixed oldStepSize = dt[0].y;
         if (oldStepSize > 0.0f)
             newStepSize = min(newStepSize, oldStepSize*2.0f); // For safety, limit how quickly dt can increase.
         if (newStepSize > oldStepSize && newStepSize < 1.1f*oldStepSize)
@@ -97,9 +108,9 @@ extern "C" __global__ void selectLangevinStepSize(real maxStepSize, real errorTo
 
         // Recalculate the integration parameters.
 
-        real vscale = exp(-newStepSize/tau);
-        real fscale = (1-vscale)*tau;
-        real noisescale = sqrt(2*kT/tau)*sqrt(0.5f*(1-vscale*vscale)*tau);
+        mixed vscale = exp(-newStepSize/tau);
+        mixed fscale = (1-vscale)*tau;
+        mixed noisescale = sqrt(2*kT/tau)*sqrt(0.5f*(1-vscale*vscale)*tau);
         params[VelScale] = vscale;
         params[ForceScale] = fscale;
         params[NoiseScale] = noisescale;

platforms/cuda/src/kernels/removeCM.cu

@@ -2,13 +2,13 @@
  * Calculate the center of mass momentum.
  */
 
-extern "C" __global__ void calcCenterOfMassMomentum(int numAtoms, const real4* __restrict__ velm, float4* __restrict__ cmMomentum) {
+extern "C" __global__ void calcCenterOfMassMomentum(int numAtoms, const mixed4* __restrict__ velm, float4* __restrict__ cmMomentum) {
     extern __shared__ volatile float3 temp[];
     float3 cm = make_float3(0, 0, 0);
     for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < numAtoms; index += blockDim.x*gridDim.x) {
-        real4 velocity = velm[index];
-        if (velocity.w != 0.0) {
-            real mass = RECIP(velocity.w);
+        mixed4 velocity = velm[index];
+        if (velocity.w != 0) {
+            mixed mass = RECIP(velocity.w);
             cm.x += (float) velocity.x*mass;
             cm.y += (float) velocity.y*mass;
             cm.z += (float) velocity.z*mass;
@@ -57,7 +57,7 @@ extern "C" __global__ void calcCenterOfMassMomentum(int numAtoms, const real4* _
  * Remove center of mass motion.
  */
 
-extern "C" __global__ void removeCenterOfMassMomentum(unsigned int numAtoms, real4* __restrict__ velm, const float4* __restrict__ cmMomentum) {
+extern "C" __global__ void removeCenterOfMassMomentum(unsigned int numAtoms, mixed4* __restrict__ velm, const float4* __restrict__ cmMomentum) {
     // First sum all of the momenta that were calculated by individual groups.
 
     extern volatile float3 temp[];
@@ -104,7 +104,7 @@ extern "C" __global__ void removeCenterOfMassMomentum(unsigned int numAtoms, rea
     // Now remove the center of mass velocity from each atom.
 
     for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < numAtoms; index += blockDim.x*gridDim.x) {
-        real4 velocity = velm[index];
+        mixed4 velocity = velm[index];
         velocity.x -= cm.x;
         velocity.y -= cm.y;
         velocity.z -= cm.z;

platforms/cuda/src/kernels/verlet.cu

@@ -2,15 +2,22 @@
  * Perform the first step of Verlet integration.
  */
 
-extern "C" __global__ void integrateVerletPart1(const real2* __restrict__ dt, const real4* __restrict__ posq, real4* __restrict__ velm, const long long* __restrict__ force, real4* __restrict__ posDelta) {
-    const real2 stepSize = dt[0];
-    const real dtPos = stepSize.y;
-    const real dtVel = 0.5f*(stepSize.x+stepSize.y);
-    const real scale = dtVel/(real) 0xFFFFFFFF;
+extern "C" __global__ void integrateVerletPart1(const mixed2* __restrict__ dt, const real4* __restrict__ posq,
+        const real4* __restrict__ posqCorrection, mixed4* __restrict__ velm, const long long* __restrict__ force, mixed4* __restrict__ posDelta) {
+    const mixed2 stepSize = dt[0];
+    const mixed dtPos = stepSize.y;
+    const mixed dtVel = 0.5f*(stepSize.x+stepSize.y);
+    const mixed scale = dtVel/(mixed) 0xFFFFFFFF;
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_ATOMS; index += blockDim.x*gridDim.x) {
-        real4 velocity = velm[index];
+        mixed4 velocity = velm[index];
         if (velocity.w != 0.0) {
+#ifdef USE_MIXED_PRECISION
+            real4 pos1 = posq[index];
+            real4 pos2 = posqCorrection[index];
+            mixed4 pos = make_mixed4(pos1.x+(mixed)pos2.x, pos1.y+(mixed)pos2.y, pos1.z+(mixed)pos2.z, pos1.w);
+#else
             real4 pos = posq[index];
+#endif
             velocity.x += scale*force[index]*velocity.w;
             velocity.y += scale*force[index+PADDED_NUM_ATOMS]*velocity.w;
             velocity.z += scale*force[index+PADDED_NUM_ATOMS*2]*velocity.w;
@@ -27,22 +34,34 @@ extern "C" __global__ void integrateVerletPart1(const real2* __restrict__ dt, co
  * Perform the second step of Verlet integration.
  */
 
-extern "C" __global__ void integrateVerletPart2(real2* __restrict__ dt, real4* __restrict__ posq, real4* __restrict__ velm, const real4* __restrict__ posDelta) {
-    real2 stepSize = dt[0];
+extern "C" __global__ void integrateVerletPart2(mixed2* __restrict__ dt, real4* __restrict__ posq,
+        real4* __restrict__ posqCorrection, mixed4* __restrict__ velm, const mixed4* __restrict__ posDelta) {
+    mixed2 stepSize = dt[0];
     double oneOverDt = 1.0/stepSize.y;
     int index = blockIdx.x*blockDim.x+threadIdx.x;
     if (index == 0)
         dt[0].x = stepSize.y;
     for (; index < NUM_ATOMS; index += blockDim.x*gridDim.x) {
-        real4 velocity = velm[index];
+        mixed4 velocity = velm[index];
         if (velocity.w != 0.0) {
+#ifdef USE_MIXED_PRECISION
+            real4 pos1 = posq[index];
+            real4 pos2 = posqCorrection[index];
+            mixed4 pos = make_mixed4(pos1.x+(mixed)pos2.x, pos1.y+(mixed)pos2.y, pos1.z+(mixed)pos2.z, pos1.w);
+#else
             real4 pos = posq[index];
-            real4 delta = posDelta[index];
+#endif
+            mixed4 delta = posDelta[index];
             pos.x += delta.x;
             pos.y += delta.y;
             pos.z += delta.z;
-            velocity = make_real4((real) (delta.x*oneOverDt), (real) (delta.y*oneOverDt), (real) (delta.z*oneOverDt), velocity.w);
+            velocity = make_mixed4((mixed) (delta.x*oneOverDt), (mixed) (delta.y*oneOverDt), (mixed) (delta.z*oneOverDt), velocity.w);
+#ifdef USE_MIXED_PRECISION
+            posq[index] = make_real4((real) pos.x, (real) pos.y, (real) pos.z, (real) pos.w);
+            posqCorrection[index] = make_real4(pos.x-(real) pos.x, pos.y-(real) pos.y, pos.z-(real) pos.z, 0);
+#else
             posq[index] = pos;
+#endif
             velm[index] = velocity;
         }
     }
@@ -52,15 +71,15 @@ extern "C" __global__ void integrateVerletPart2(real2* __restrict__ dt, real4* _
  * Select the step size to use for the next step.
  */
 
-extern "C" __global__ void selectVerletStepSize(real maxStepSize, real errorTol, real2* __restrict__ dt, const real4* __restrict__ velm, const long long* __restrict__ force) {
+extern "C" __global__ void selectVerletStepSize(mixed maxStepSize, mixed errorTol, mixed2* __restrict__ dt, const mixed4* __restrict__ velm, const long long* __restrict__ force) {
     // Calculate the error.
 
-    extern __shared__ real error[];
-    real err = 0.0f;
-    const real scale = RECIP((real) 0xFFFFFFFF);
+    extern __shared__ mixed error[];
+    mixed err = 0.0f;
+    const mixed scale = RECIP((mixed) 0xFFFFFFFF);
     for (int index = threadIdx.x; index < NUM_ATOMS; index += blockDim.x*gridDim.x) {
-        real3 f = make_real3(scale*force[index], scale*force[index+PADDED_NUM_ATOMS], scale*force[index+PADDED_NUM_ATOMS*2]);
-        real invMass = velm[index].w;
+        mixed3 f = make_mixed3(scale*force[index], scale*force[index+PADDED_NUM_ATOMS], scale*force[index+PADDED_NUM_ATOMS*2]);
+        mixed invMass = velm[index].w;
         err += (f.x*f.x + f.y*f.y + f.z*f.z)*invMass;
     }
     error[threadIdx.x] = err;
@@ -74,9 +93,9 @@ extern "C" __global__ void selectVerletStepSize(real maxStepSize, real errorTol,
         __syncthreads();
     }
     if (threadIdx.x == 0) {
-        real totalError = SQRT(error[0]/(NUM_ATOMS*3));
-        real newStepSize = SQRT(errorTol/totalError);
-        real oldStepSize = dt[0].y;
+        mixed totalError = sqrt(error[0]/(NUM_ATOMS*3));
+        mixed newStepSize = sqrt(errorTol/totalError);
+        mixed oldStepSize = dt[0].y;
         if (oldStepSize > 0.0f)
             newStepSize = min(newStepSize, oldStepSize*2.0f); // For safety, limit how quickly dt can increase.
         if (newStepSize > oldStepSize && newStepSize < 1.1f*oldStepSize)

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp

@@ -1410,11 +1410,12 @@ void CudaCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextImpl&
     }
 }
 
-template <class T, class T4>
+template <class T, class T4, class M4>
 void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextImpl& context, vector<double>& outputMultipoleMoments) {
     // Compute the local coordinates relative to the center of mass.
     int numAtoms = cu.getNumAtoms();
-    vector<T4> posq, velm;
+    vector<T4> posq;
+    vector<M4> velm;
     cu.getPosq().download(posq);
     cu.getVelm().download(velm);
     double totalMass = 0.0;
@@ -1524,9 +1525,11 @@ void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextIm
 void CudaCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, vector<double>& outputMultipoleMoments) {
     context.calcForcesAndEnergy(false, false, -1);
     if (cu.getUseDoublePrecision())
-        computeSystemMultipoleMoments<double, double4>(context, outputMultipoleMoments);
+        computeSystemMultipoleMoments<double, double4, double4>(context, outputMultipoleMoments);
+    else if (cu.getUseMixedPrecision())
+        computeSystemMultipoleMoments<float, float4, double4>(context, outputMultipoleMoments);
     else
-        computeSystemMultipoleMoments<float, float4>(context, outputMultipoleMoments);
+        computeSystemMultipoleMoments<float, float4, float4>(context, outputMultipoleMoments);
 }
 
 /* -------------------------------------------------------------------------- *

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h

@@ -321,7 +321,7 @@ private:
         const char* getSortKey() const {return "value.y";}
     };
     void initializeScaleFactors();
-    template <class T, class T4> void computeSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
+    template <class T, class T4, class M4> void computeSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
     int numMultipoles, maxInducedIterations;
     double inducedEpsilon;
     bool hasInitializedScaleFactors, hasInitializedFFT;