Add Drude Nose-Hoover capability, with tests

olla/include/openmm/kernels.h

@@ -1365,8 +1365,9 @@ public:
      *
      * @param context the context in which to execute this kernel
      * @param noseHooverChain the chain whose energy is to be determined.
+     * @param downloadValue whether the computed value should be downloaded and returned.
      */
-    virtual double computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain) = 0;
+    virtual double computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue) = 0;
     /**
      * Execute the kernel that scales the velocities of particles associated with a nose hoover chain
      *

openmmapi/include/openmm/VelocityVerletIntegrator.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "Integrator.h"
+#include "openmm/State.h"
 #include "openmm/Kernel.h"
 #include "openmm/NoseHooverChain.h"
 #include "internal/windowsExport.h"
@@ -174,10 +175,10 @@ protected:
     /**
      * Compute the kinetic energy of the system at the current time.
      */
-    double computeKineticEnergy();
+    virtual double computeKineticEnergy();
+
     std::vector<NoseHooverChain> noseHooverChains;
     bool forcesAreValid;
-private:
     Kernel vvKernel, nhcKernel;
 };
 

openmmapi/src/VelocityVerletIntegrator.cpp

@@ -38,9 +38,9 @@
 #include "openmm/CMMotionRemover.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/kernels.h"
+#include <iostream>
 #include <string>
 #include <algorithm>
-#include <iostream>
 
 using namespace OpenMM;
 using std::string;
@@ -95,8 +95,8 @@ int VelocityVerletIntegrator::addMaskedNoseHooverChainThermostat(System& system,
         int particle1, particle2;
         double distance;
         system.getConstraintParameters(constraintNum, particle1, particle2, distance);
-        bool particle1_in_mask = (std::find(mask.begin(), mask.end(), particle1) == mask.end());
-        bool particle2_in_mask = (std::find(mask.begin(), mask.end(), particle2) == mask.end());
+        bool particle1_in_mask = (std::find(mask.begin(), mask.end(), particle1) != mask.end());
+        bool particle2_in_mask = (std::find(mask.begin(), mask.end(), particle2) != mask.end());
         if ((system.getParticleMass(particle1) > 0) && (system.getParticleMass(particle2) > 0)){
             if ((particle1_in_mask && !particle2_in_mask) || (!particle1_in_mask && particle2_in_mask)){
                 throw OpenMMException("Cannot add only one of particles " + std::to_string(particle1) + " and " + std::to_string(particle2)
@@ -173,8 +173,17 @@ void VelocityVerletIntegrator::setCollisionFrequency(double frequency, int chain
 }
 
 double VelocityVerletIntegrator::computeKineticEnergy() {
-    return vvKernel.getAs<IntegrateVelocityVerletStepKernel>().computeKineticEnergy(*context, *this);
-
+    double kE = 0.0;
+    if(noseHooverChains.size()) {
+        for (const auto &nhc: noseHooverChains){
+            if (nhc.getParentAtoms().size() == 0) {
+                kE += nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, true);
+            }
+        }
+    } else {
+        kE = vvKernel.getAs<IntegrateVelocityVerletStepKernel>().computeKineticEnergy(*context, *this);
+    }
+    return kE;
 }
 
 double VelocityVerletIntegrator::computeHeatBathEnergy() {
@@ -216,13 +225,13 @@ void VelocityVerletIntegrator::step(int steps) {
     for (int i = 0; i < steps; ++i) {
         context->updateContextState();
         for(auto &nhc : noseHooverChains) {
-            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc);
+            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, false);
             scale = nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, nhc, kineticEnergy, getStepSize());
             nhcKernel.getAs<NoseHooverChainKernel>().scaleVelocities(*context, nhc, scale);
         }
         vvKernel.getAs<IntegrateVelocityVerletStepKernel>().execute(*context, *this, forcesAreValid);
         for(auto &nhc : noseHooverChains) {
-            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc);
+            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, false);
             scale = nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, nhc, kineticEnergy, getStepSize());
             nhcKernel.getAs<NoseHooverChainKernel>().scaleVelocities(*context, nhc, scale);
         }

platforms/cuda/include/CudaArray.h

@@ -71,10 +71,6 @@ public:
      * @param name              the name of the array
      */
     CudaArray(CudaContext& context, int size, int elementSize, const std::string& name);
-    CudaArray(const CudaArray &other) noexcept;
-    CudaArray& operator=(CudaArray &&other) noexcept;
-    CudaArray& operator=(const CudaArray& other) noexcept;
-    CudaArray(CudaArray &&other) noexcept ; 
     ~CudaArray();
     /**
      * Initialize this object.

platforms/cuda/include/CudaIntegrationUtilities.h

@@ -30,6 +30,7 @@
 #include "openmm/System.h"
 #include "CudaContext.h"
 #include "windowsExportCuda.h"
+#include <map>
 #include <iosfwd>
 
 namespace OpenMM {
@@ -125,7 +126,7 @@ public:
      * 
      * @return  vector of chain states
      */
-    std::vector<CudaArray>& getNoseHooverChainState();
+    std::map<int, CudaArray>& getNoseHooverChainState();
 private:
     void applyConstraints(bool constrainVelocities, double tol);
     CudaContext& context;
@@ -174,7 +175,7 @@ private:
     double2 lastStepSize;
     struct ShakeCluster;
     struct ConstraintOrderer;
-    std::vector<CudaArray> noseHooverChainState;
+    std::map<int, CudaArray> noseHooverChainState;
 };
 
 } // namespace OpenMM

platforms/cuda/include/CudaKernels.h

@@ -1750,8 +1750,10 @@ public:
      *
      * @param context the context in which to execute this kernel
      * @param noseHooverChain the chain whose energy is to be determined.
+     * @param downloadValue whether the computed value should be downloaded and returned.
+     *
      */
-     virtual double computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain);
+     virtual double computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue);
 
     /**
      * Execute the kernel that scales the velocities of particles associated with a nose hoover chain
@@ -1765,13 +1767,12 @@ public:
 private:
     CudaContext& cu;
     CudaArray scaleFactorBuffer, kineticEnergyBuffer, chainMasses, chainForces, heatBathEnergy;
-    std::vector<CudaArray> masks;
+    std::map<int, CudaArray> masks;
     std::map<int, CUfunction> propagateKernels;
     CUfunction reduceEnergyKernel;
     CUfunction computeHeatBathEnergyKernel;
     CUfunction computeMaskedKineticEnergyKernel;
     CUfunction scaleVelocitiesKernel;
-    CUfunction zeroEnergyBuffersKernel;
 };
 
 /**

platforms/cuda/src/CudaArray.cpp

@@ -35,47 +35,6 @@ using namespace OpenMM;
 CudaArray::CudaArray() : pointer(0), ownsMemory(false) {
 }
 
-CudaArray::CudaArray(const CudaArray &other) noexcept {
-    context = other.context;
-    pointer = other.pointer;
-    size = other.size;
-    elementSize = other.elementSize;
-    ownsMemory = false;
-    name = other.name;
-}
-
-CudaArray::CudaArray(CudaArray &&other) noexcept :
-    context(std::move(other.context)),
-    pointer(std::move(other.pointer)),
-    size(std::move(other.size)),
-    elementSize(std::move(other.elementSize)),
-    ownsMemory(false),
-    name(std::move(other.name)) { }
-
-CudaArray& CudaArray::operator =(const CudaArray& other) noexcept{
-    if(this != &other) {
-        context = other.context;
-        pointer = other.pointer;
-        size = other.size;
-        elementSize = other.elementSize;
-        ownsMemory = false;
-        name = other.name;
-    }
-    return *this;
-}
-
-CudaArray& CudaArray::operator =(CudaArray &&other) noexcept{
-    if(this != &other) {
-        context = other.context;
-        pointer = other.pointer;
-        size = other.size;
-        elementSize = other.elementSize;
-        ownsMemory = false;
-        name = other.name;
-    }
-    return *this;
-}
-
 CudaArray::CudaArray(CudaContext& context, int size, int elementSize, const std::string& name) : pointer(0) {
     initialize(context, size, elementSize, name);
 }

platforms/cuda/src/CudaIntegrationUtilities.cpp

@@ -707,6 +707,24 @@ int CudaIntegrationUtilities::prepareRandomNumbers(int numValues) {
 }
 
 void CudaIntegrationUtilities::createCheckpoint(ostream& stream) {
+    size_t numChains = noseHooverChainState.size();
+    bool useDouble = context.getUseDoublePrecision() || context.getUseMixedPrecision();
+    stream.write((char*) &numChains, sizeof(size_t));
+    for (auto &chainState: noseHooverChainState){
+        int chainID = chainState.first;
+        size_t chainLength = chainState.second.getSize();
+        stream.write((char*) &chainID, sizeof(int));
+        stream.write((char*) &chainLength, sizeof(size_t));
+        if (useDouble) {
+            vector<double2> stateVec;
+            chainState.second.download(stateVec);
+            stream.write((char*) stateVec.data(), sizeof(double2)*chainLength);
+        } else {
+            vector<float2> stateVec;
+            chainState.second.download(stateVec);
+            stream.write((char*) stateVec.data(), sizeof(float2)*chainLength);
+        }
+    }
     if (!random.isInitialized()) 
         return;
     stream.write((char*) &randomPos, sizeof(int));
@@ -716,18 +734,31 @@ void CudaIntegrationUtilities::createCheckpoint(ostream& stream) {
     vector<int4> randomSeedVec;
     randomSeed.download(randomSeedVec);
     stream.write((char*) &randomSeedVec[0], sizeof(int4)*randomSeed.getSize());
-    size_t numChains = noseHooverChainState.size();
-    stream.write((char*) &numChains, sizeof(size_t));
-    for (auto &chainState: noseHooverChainState){
-        vector<float2> stateVec;
-        chainState.download(stateVec);
-        size_t vecLength = stateVec.size();
-        stream.write((char*) &vecLength, sizeof(size_t));
-        stream.write((char*) stateVec.data(), sizeof(float2)*stateVec.size());
-    }
 }
 
 void CudaIntegrationUtilities::loadCheckpoint(istream& stream) {
+    size_t numChains, chainLength;
+    bool useDouble = context.getUseDoublePrecision() || context.getUseMixedPrecision();
+    stream.read((char*) &numChains, sizeof(size_t));
+    noseHooverChainState.clear();
+    for (size_t i=0; i<numChains; i++){
+        int chainID;
+        stream.read((char*) &chainID, sizeof(int));
+        stream.read((char*) &chainLength, sizeof(size_t));
+        if (useDouble) {
+            noseHooverChainState[chainID] = CudaArray();
+            noseHooverChainState[chainID].initialize<double2>(context, chainLength, "chainState" + std::to_string(chainID));
+            std::vector<double2> stateVec(chainLength);
+            stream.read((char*) &stateVec[0], sizeof(double2)*chainLength);
+            noseHooverChainState[chainID].upload(stateVec);
+        } else {
+            noseHooverChainState[chainID] = CudaArray();
+            noseHooverChainState[chainID].initialize<float2>(context, chainLength, "chainState" + std::to_string(chainID));
+            std::vector<float2> stateVec(chainLength);
+            stream.read((char*) &stateVec[0], sizeof(float2)*chainLength);
+            noseHooverChainState[chainID].upload(stateVec);
+        }
+    }
     if (!random.isInitialized()) 
         return;
     stream.read((char*) &randomPos, sizeof(int));
@@ -737,17 +768,6 @@ void CudaIntegrationUtilities::loadCheckpoint(istream& stream) {
     vector<int4> randomSeedVec(randomSeed.getSize());
     stream.read((char*) &randomSeedVec[0], sizeof(int4)*randomSeed.getSize());
     randomSeed.upload(randomSeedVec);
-    size_t numChains, chainLength;
-    stream.read((char*) &numChains, sizeof(size_t));
-    noseHooverChainState.clear();
-    for (size_t i=0; i<numChains; i++){
-        stream.read((char*) &chainLength, sizeof(size_t));
-        std::vector<float2> stateVec(chainLength);
-        stream.read((char*) &stateVec[0], sizeof(float2)*chainLength);
-        CudaArray state;
-        state.upload(stateVec);
-        noseHooverChainState.push_back(state);
-    }
 }
 
 double CudaIntegrationUtilities::computeKineticEnergy(double timeShift) {
@@ -796,6 +816,6 @@ double CudaIntegrationUtilities::computeKineticEnergy(double timeShift) {
     return 0.5*energy;
 }
 
-std::vector<CudaArray>& CudaIntegrationUtilities::getNoseHooverChainState(){
+std::map<int, CudaArray>& CudaIntegrationUtilities::getNoseHooverChainState(){
         return noseHooverChainState;
 };

platforms/cuda/src/CudaKernels.cpp

@@ -8391,8 +8391,10 @@ double CudaNoseHooverChainKernel::propagateChain(ContextImpl& context, const Nos
     double frequency = nhc.getDefaultCollisionFrequency();
 
     auto & chainState = cu.getIntegrationUtilities().getNoseHooverChainState();
-    if (chainID >= chainState.size()) chainState.resize(chainID+1);
-    if (!chainState.at(chainID).isInitialized() || chainState.at(chainID).getSize() != chainLength) {
+    if (!chainState.count(chainID)) {
+        chainState[chainID] = CudaArray();
+    }
+    if (chainState.at(chainID).getSize() != chainLength) {
         // We need to upload the CUDA array
         if(useDouble){
             chainState.at(chainID).initialize<double2>(cu, chainLength, "chainState" + std::to_string(chainID));
@@ -8504,7 +8506,7 @@ double CudaNoseHooverChainKernel::computeHeatBathEnergy(ContextImpl& context, co
     }
 }
 
-double CudaNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &nhc) {
+double CudaNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &nhc, bool downloadValue) {
     cu.setAsCurrent();
 
     bool useDouble = cu.getUseDoublePrecision() || cu.getUseMixedPrecision();
@@ -8512,9 +8514,8 @@ double CudaNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& contex
     int chainID = nhc.getDefaultChainID();
     const auto & parents = nhc.getParentAtoms();
     const auto & thermostatedAtoms = nhc.getThermostatedAtoms();
-    if (chainID >= masks.size()) masks.resize(chainID+1);
     auto nAtoms = cu.getPaddedNumAtoms();
-    if (!masks.at(chainID).isInitialized()) { 
+    if (!masks.count(chainID)) { 
         // We need to upload the CUDA array
         std::vector<int> maskVec(nAtoms);
         for (int i = 0; i < nAtoms; ++i) maskVec[i] = i;
@@ -8530,12 +8531,9 @@ double CudaNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& contex
                 maskVec[atom] = -1L;
             }
         }
-        // Account for padding in the paddedNumAtoms loops
-        for(int i = thermostatedAtoms.size(); i < nAtoms; ++i){
-            maskVec[i] = i;
-        }
-        masks.at(chainID).initialize<int>(cu, nAtoms, "mask" + std::to_string(chainID));
-        masks.at(chainID).upload(maskVec);
+        masks[chainID] = CudaArray();
+        masks[chainID].initialize<int>(cu, nAtoms, "mask" + std::to_string(chainID));
+        masks[chainID].upload(maskVec);
     }
     if (masks[chainID].getSize() != nAtoms) {
         throw OpenMMException("Number of atoms changed. Cannot be handled by the same Nose-Hoover thermostat.");
@@ -8561,7 +8559,13 @@ double CudaNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& contex
     void* args2[] = {&cu.getEnergyBuffer().getDevicePointer(), &kineticEnergyBuffer.getDevicePointer(), &bufferSize, &workGroupSize};
     cu.executeKernel(reduceEnergyKernel, args2, workGroupSize, workGroupSize, workGroupSize*cu.getEnergyBuffer().getElementSize());
 
-    return 0;
+    double value = 0;
+    if (downloadValue) {
+        void * pinnedBuffer = cu.getPinnedBuffer();
+        kineticEnergyBuffer.download(pinnedBuffer);
+        value = useDouble ? *((double*) pinnedBuffer) : *((float*) pinnedBuffer);
+    }
+    return value;
 }
 
 void CudaNoseHooverChainKernel::scaleVelocities(ContextImpl& context, const NoseHooverChain &nhc, double scaleFactor) {

platforms/reference/include/ReferenceKernels.h

@@ -1465,8 +1465,10 @@ public:
      *
      * @param context the context in which to execute this kernel
      * @param noseHooverChain the chain whose energy is to be determined.
+     * @param downloadValue whether the computed value should be downloaded and returned.
+     *
      */
-     virtual double computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain);
+     virtual double computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue);
 
     /**
      * Execute the kernel that scales the velocities of particles associated with a nose hoover chain

platforms/reference/src/ReferenceKernels.cpp

@@ -2527,7 +2527,7 @@ double ReferenceNoseHooverChainKernel::computeHeatBathEnergy(ContextImpl& contex
     return kineticEnergy + potentialEnergy;
 }
 
-double ReferenceNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain) {
+double ReferenceNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue) {
     const std::vector<int>& mask = noseHooverChain.getThermostatedAtoms();
     const std::vector<int>& parents = noseHooverChain.getParentAtoms();
     std::vector<Vec3>& velocities = extractVelocities(context);
@@ -2554,6 +2554,7 @@ double ReferenceNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& c
             ke += 0.5 * mass * velocity.dot(velocity);
         }
     }
+    // We ignore the downloadValue argument here and always return the correct value
     return ke;
 }
 

plugins/drude/openmmapi/include/openmm/DrudeVelocityVerletIntegrator.h

@@ -76,13 +76,6 @@ public:
     int addDrudeNoseHooverChainThermostat(System& system, double temperature, double collisionFrequency,
                                           double drudeTemperature, double drudeCollisionFrequency,
                                           int chainLength, int numMTS, int numYoshidaSuzuki);
-    /**
-     * Advance a simulation through time by taking a series of time steps.
-     *
-     * @param steps   the number of time steps to take
-     */
-    //void step(int steps);
-//protected:
     /**
      * This will be called by the Context when it is created.  It informs the Integrator
      * of what context it will be integrating, and gives it a chance to do any necessary initialization.
@@ -90,9 +83,13 @@ public:
      */
     void initialize(ContextImpl& context);
     /**
-     * Compute the kinetic energy of the system at the current time.
+     * Compute the kinetic energy of the drude particles at the current time.
+     */
+    double computeDrudeKineticEnergy();
+    /**
+     * Compute the kinetic energy of all (real and drude) particles at the current time.
      */
-    //double computeKineticEnergy();
+    double computeTotalKineticEnergy();
 };
 
 } // namespace OpenMM

plugins/drude/openmmapi/src/DrudeVelocityVerletIntegrator.cpp

@@ -35,7 +35,10 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/DrudeKernels.h"
 #include "openmm/DrudeForce.h"
+#include "openmm/CMMotionRemover.h"
+#include "openmm/kernels.h"
 #include <ctime>
+#include <iostream>
 #include <string>
 #include <set>
 
@@ -63,14 +66,13 @@ int DrudeVelocityVerletIntegrator::addDrudeNoseHooverChainThermostat(System& sys
     std::set<int> realParticlesSet; 
     vector<int> realParticles, drudeParticles, drudeParents;
     for (int i = 0; i < system.getNumParticles(); i++) {
-        realParticlesSet.insert(i); 
+        if (system.getParticleMass(i) > 0.0) realParticlesSet.insert(i); 
     }
     for (int i = 0; i < drudeForce->getNumParticles(); i++) {
         int p, p1, p2, p3, p4;
         double charge, polarizability, aniso12, aniso34;
         drudeForce->getParticleParameters(i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
         realParticlesSet.erase(p);
-        realParticlesSet.erase(p1);
         drudeParticles.push_back(p);
         drudeParents.push_back(p1);
     }
@@ -84,19 +86,44 @@ int DrudeVelocityVerletIntegrator::addDrudeNoseHooverChainThermostat(System& sys
 }
 
 void DrudeVelocityVerletIntegrator::initialize(ContextImpl& contextRef) {
+    VelocityVerletIntegrator::initialize(contextRef);
     if (owner != NULL && &contextRef.getOwner() != owner)
         throw OpenMMException("This Integrator is already bound to a context");
     const DrudeForce* drudeForce = NULL;
     const System& system = contextRef.getSystem();
-    for (int i = 0; i < system.getNumForces(); i++)
+    bool hasCMMotionRemover = false;
+    for (int i = 0; i < system.getNumForces(); i++){
         if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
             if (drudeForce == NULL)
                 drudeForce = dynamic_cast<const DrudeForce*>(&system.getForce(i));
             else
                 throw OpenMMException("The System contains multiple DrudeForces");
         }
+        if (dynamic_cast<const CMMotionRemover*>(&system.getForce(i))) {
+            hasCMMotionRemover = true;
+        }
+    }
     if (drudeForce == NULL)
         throw OpenMMException("The System does not contain a DrudeForce");
+    if (!hasCMMotionRemover) {
+        std::cout << "Warning: Did not find a center-of-mass motion remover in the system. "
+                     "This is problematic when using Drude." << std::endl;
+    }
     context = &contextRef;
     owner = &contextRef.getOwner();
 }
+
+double DrudeVelocityVerletIntegrator::computeDrudeKineticEnergy() {
+    double kE = 0.0;
+    for (const auto &nhc: noseHooverChains){
+        if (nhc.getParentAtoms().size() != 0) {
+            kE += nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, true);
+        }
+    }
+    return kE;
+}
+
+double DrudeVelocityVerletIntegrator::computeTotalKineticEnergy() {
+    return vvKernel.getAs<IntegrateVelocityVerletStepKernel>().computeKineticEnergy(*context, *this);
+}
+

plugins/drude/platforms/cuda/tests/CMakeLists.txt

@@ -5,6 +5,7 @@
 ENABLE_TESTING()
 
 INCLUDE_DIRECTORIES(${CUDA_INCLUDE_DIR})
+INCLUDE_DIRECTORIES(${OPENMM_DIR}/plugins/drude/tests)
 
 # Automatically create tests using files named "Test*.cpp"
 FILE(GLOB TEST_PROGS "*Test*.cpp")

plugins/drude/platforms/cuda/tests/TestCudaDrudeNoseHoover.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2013 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+//#include "ReferenceTests.h"
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/Platform.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "openmm/DrudeForce.h"
+#include "CudaPlatform.h"
+#include "SimTKOpenMMUtilities.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+extern "C" OPENMM_EXPORT void registerDrudeCudaKernelFactories();
+
+//OpenMM::CudaPlatform platform; 
+
+
+Platform& initializePlatform(int argc, char* argv[]) {
+    registerDrudeCudaKernelFactories();
+    if (argc > 1) Platform::getPlatformByName("CUDA").setPropertyDefaultValue("Precision", std::string(argv[1]));   
+    return Platform::getPlatformByName("CUDA");
+}
+
+#include "TestDrudeNoseHoover.h"
+
+void runPlatformTests() { }
+

plugins/drude/platforms/reference/tests/CMakeLists.txt

@@ -5,6 +5,7 @@ ENABLE_TESTING()
 INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/include)
 INCLUDE_DIRECTORIES(${OPENMM_DIR}/openmmapi/include/openmm)
 INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/reference/src)
+INCLUDE_DIRECTORIES(${OPENMM_DIR}/plugins/drude/tests)
 
 SET(SHARED_OPENMM_DRUDE_TARGET OpenMMDrude)
 

plugins/drude/platforms/reference/tests/TestReferenceDrudeNoseHoover.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2013 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+//#include "ReferenceTests.h"
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/Platform.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "openmm/DrudeForce.h"
+#include "SimTKOpenMMUtilities.h"
+#include "ReferencePlatform.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+extern "C" OPENMM_EXPORT void registerDrudeReferenceKernelFactories();
+
+
+Platform& initializePlatform(int argc, char* argv[]) {
+    registerDrudeReferenceKernelFactories();
+    return Platform::getPlatformByName("Reference");
+}
+
+#include "TestDrudeNoseHoover.h"
+
+void runPlatformTests() { }
+

tests/TestNoseHooverThermostat.h

@@ -265,7 +265,8 @@ void testPropagateChainConsistentWithPythonReference() {
 
 
 void testCheckpoints() {
-    VelocityVerletIntegrator integrator(0.001);
+    double timeStep = 0.001;
+    VelocityVerletIntegrator integrator(timeStep), newIntegrator(timeStep);
     System system;
     double mass = 1;
     system.addParticle(mass);
@@ -278,11 +279,17 @@ void testCheckpoints() {
     double temperature=300, collisionFrequency=1, chainLength=3, numMTS=3, numYS=3;
     integrator.addMaskedNoseHooverChainThermostat(system, std::vector<int>(1,0), std::vector<int>(), temperature, collisionFrequency,
                                                                   chainLength, numMTS, numYS);
+    newIntegrator.addMaskedNoseHooverChainThermostat(system, std::vector<int>(1,0), std::vector<int>(), temperature, collisionFrequency,
+                                                                  chainLength, numMTS, numYS);
     chainLength = 10;
     integrator.addMaskedNoseHooverChainThermostat(system, std::vector<int>(1,1),  std::vector<int>(1,0), 
                                                                   temperature, collisionFrequency,
                                                                   chainLength, numMTS, numYS);
+    newIntegrator.addMaskedNoseHooverChainThermostat(system, std::vector<int>(1,1),  std::vector<int>(1,0), 
+                                                                  temperature, collisionFrequency,
+                                                                  chainLength, numMTS, numYS);
     Context context(system, integrator, platform);
+    Context newContext(system, newIntegrator, platform);
     std::vector<Vec3> positions(2);
     positions[1] = {0.1,0.1,0.1};
     context.setPositions(positions);
@@ -304,16 +311,16 @@ void testCheckpoints() {
 #if DEBUG
     std::cout << std::endl << std::endl << "loading checkpoint" << std::endl;
 #endif
-    context.loadCheckpoint(checkpoint);
+    newContext.loadCheckpoint(checkpoint);
 
     
     State state2 = context.getState(State::Positions | State::Velocities);
     for (size_t i=0; i<100; i++){
-        state2 = context.getState(State::Positions | State::Velocities);
+        state2 = newContext.getState(State::Positions | State::Velocities);
 #if DEBUG
         std::cout << "posvel" << state2.getPositions()[0] << " " << state2.getVelocities()[0] << std::endl;
 #endif
-        integrator.step(1);
+        newIntegrator.step(1);
     }
 
     for (int i=0; i<3;i++){

wrappers/generateWrappers.py

@@ -78,6 +78,7 @@ class WrapperGenerator:
                             'Vec3 OpenMM::LocalCoordinatesSite::getXWeights',
                             'Vec3 OpenMM::LocalCoordinatesSite::getYWeights',
                             'std::vector<double> OpenMM::NoseHooverChain::getDefaultYoshidaSuzukiWeights',
+                            'virtual void OpenMM::VelocityVerletIntegrator::stateChanged',
                            ] 
         self.hideClasses = ['Kernel', 'KernelImpl', 'KernelFactory', 'ContextImpl', 'SerializationNode', 'SerializationProxy']
         self.nodeByID={}