Remove system from NHC constructor

openmmapi/include/openmm/NoseHooverIntegrator.h

@@ -41,7 +41,6 @@
 namespace OpenMM {
 
 class System;
-
 /**
  * This is an Integrator which simulates a System using one or more Nose Hoover chain
  * thermostats, using the velocity Verlet propagation algorithm.
@@ -60,8 +59,6 @@ public:
      * Create a NoseHooverIntegrator.
      *
      * @param stepSize the step size with which to integrate the system (in picoseconds)
-     * @param system the system to be thermostated.  Note: this must be setup, i.e. all
-     *        particles should have been added, before calling this function.
      * @param temperature the target temperature for the system.
      * @param collisionFrequency the frequency of the interaction with the heat bath (in 1/ps).
      * @param chainLength the number of beads in the Nose-Hoover chain.
@@ -69,7 +66,7 @@ public:
      * @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
      *        used in the chain propagation algorithm (must be 1, 3, or 5).
      */
-    explicit NoseHooverIntegrator(double stepSize, System &system, double temperature, int collisionFrequnency,
+    explicit NoseHooverIntegrator(double stepSize, double temperature, int collisionFrequnency,
                                   int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 3);
 
     virtual ~NoseHooverIntegrator();
@@ -82,8 +79,6 @@ public:
     /**
      * Add a simple Nose-Hoover Chain thermostat to control the temperature of the full system
      *
-     * @param system the system to be thermostated.  Note: this must be setup, i.e. all
-     *        particles should have been added, before calling this function.
      * @param temperature the target temperature for the system.
      * @param collisionFrequency the frequency of the interaction with the heat bath (in 1/ps).
      * @param chainLength the number of beads in the Nose-Hoover chain.
@@ -91,7 +86,7 @@ public:
      * @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
      *        used in the chain propagation algorithm (must be 1, 3, or 5).
      */
-     int addThermostat(System& system, double temperature, double collisionFrequency,
+     int addThermostat(double temperature, double collisionFrequency,
                        int chainLength, int numMTS, int numYoshidaSuzuki);
     /**
      * Add a Nose-Hoover Chain thermostat to control the temperature of a collection of atoms and/or pairs of
@@ -100,8 +95,6 @@ public:
      * connected atoms may be provided; in this case both the center of mass absolute motion of each pair is
      * controlled as well as their motion relative to each other, which is independently thermostated.
      *
-     * @param system the system to be thermostated.  Note: this must be setup, i.e. all
-     *        particles should have been added, before calling this function.
      * @param thermostatedParticles list of particle ids to be thermostated.
      * @param thermostatedPairs a list of pairs of connected atoms whose absolute center of mass motion
      *        and motion relative to one another will be independently thermostated.
@@ -116,7 +109,7 @@ public:
      * @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
      *        used in the chain propagation algorithm (must be 1, 3, or 5).
      */
-     int addSubsystemThermostat(System& system, const std::vector<int>& thermostatedParticles,
+     int addSubsystemThermostat(const std::vector<int>& thermostatedParticles,
                                 const std::vector< std::pair< int, int> >& thermostatedPairs,
                                 double temperature, double relativeTemperature,
                                 double collisionFrequency, double relativeCollisionFrequency,
@@ -221,6 +214,10 @@ protected:
      * It will also get called again if the application calls reinitialize() on the Context.
      */
     void initialize(ContextImpl& context);
+    /**
+     * Goes through the list of requested thermostat data and creates the thermostat chains.
+     */
+    void createThermostats(const System& system);
     /**
      * This will be called by the Context when it is destroyed to let the Integrator do any necessary
      * cleanup.  It will also get called again if the application calls reinitialize() on the Context.
@@ -235,6 +232,25 @@ protected:
      */
     virtual double computeKineticEnergy();
 
+    struct ThermostatData {
+        std::vector<int> thermostatedParticles;
+        std::vector< std::pair< int, int> > thermostatedPairs;
+        double temperature;
+        double relativeTemperature;
+        double collisionFrequency;
+        double relativeCollisionFrequency;
+        int chainLength;
+        int numMTS;
+        int numYoshidaSuzuki;
+        ThermostatData(const std::vector<int>&particles, const std::vector<std::pair<int, int>> &pairs, double temp,
+                       double relTemp, double freq, double relFreq, int length, int MTS, int YS) :
+                               thermostatedParticles(particles), thermostatedPairs(pairs), temperature(temp),
+                               relativeTemperature(relTemp), collisionFrequency(freq), relativeCollisionFrequency(relFreq),
+                               chainLength(length), numMTS(MTS), numYoshidaSuzuki(YS) {}
+    };
+
+    std::vector<ThermostatData> thermostatData;
+
     std::vector<NoseHooverChain> noseHooverChains;
     bool forcesAreValid;
     Kernel vvKernel, nhcKernel;

openmmapi/src/NoseHooverIntegrator.cpp

@@ -52,11 +52,11 @@ NoseHooverIntegrator::NoseHooverIntegrator(double stepSize):
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
 }
-NoseHooverIntegrator::NoseHooverIntegrator(double stepSize, System &system, double temperature, int collisionFrequnency,
+NoseHooverIntegrator::NoseHooverIntegrator(double stepSize, double temperature, int collisionFrequnency,
                                            int chainLength, int numMTS, int numYoshidaSuzuki) : forcesAreValid(false) {
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
-    addThermostat(system, temperature, collisionFrequnency, chainLength, numMTS, numYoshidaSuzuki);
+    addThermostat(temperature, collisionFrequnency, chainLength, numMTS, numYoshidaSuzuki);
 }
 
 NoseHooverIntegrator::~NoseHooverIntegrator() {}
@@ -65,89 +65,133 @@ std::pair<double, double> NoseHooverIntegrator::propagateChain(std::pair<double,
     return nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, noseHooverChains.at(chainID), kineticEnergy, getStepSize());
 }
 
-int NoseHooverIntegrator::addThermostat(System& system, double temperature, double collisionFrequency,
+
+int NoseHooverIntegrator::addThermostat(double temperature, double collisionFrequency,
                                                            int chainLength, int numMTS, int numYoshidaSuzuki) {
-    std::vector<int> thermostatedParticles;
-    for(int particle = 0; particle < system.getNumParticles(); ++particle) {
-        double mass = system.getParticleMass(particle);
-        if ( (mass > 0) && (mass < 0.8) ){
-            std::cout << "Warning: Found particles with mass between 0.0 and 0.8 dalton. Did you mean to make a DrudeVelocityVerletIntegrator instead? "
-                         "The thermostat you are about to use will not treat these particles as Drude particles!" << std::endl;
-        }
-        if(system.getParticleMass(particle) > 0) {
-            thermostatedParticles.push_back(particle);
-        }
-    }
 
-    return addSubsystemThermostat(system, thermostatedParticles, std::vector<std::pair<int, int>>(), temperature, temperature,
+    return addSubsystemThermostat(std::vector<int>(), std::vector<std::pair<int, int>>(), temperature, temperature,
                                   collisionFrequency, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
 }
 
-int NoseHooverIntegrator::addSubsystemThermostat(System& system, const std::vector<int>& thermostatedParticles,
+int NoseHooverIntegrator::addSubsystemThermostat(const std::vector<int>& thermostatedParticles,
                                                  const std::vector< std::pair< int, int> > &thermostatedPairs,
                                                  double temperature, double relativeTemperature,
                                                  double collisionFrequency, double relativeCollisionFrequency,
                                                  int chainLength, int numMTS, int numYoshidaSuzuki) {
-    if (context) {
-        throw OpenMMException("Nose-Hoover chains cannot be added after binding this integrator to a context.");
-    }
-    // figure out the number of DOFs
-    int nDOF = 3*(thermostatedParticles.size() + thermostatedPairs.size());
-    for (int constraintNum = 0; constraintNum < system.getNumConstraints(); constraintNum++) {
-        int particle1, particle2;
-        double distance;
-        system.getConstraintParameters(constraintNum, particle1, particle2, distance);
-        bool particle1_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(), thermostatedParticles.end(), particle1)
-                                                    != thermostatedParticles.end())) ||
-                                                  (std::find_if(thermostatedPairs.begin(), thermostatedPairs.end(),
-                                                   [&particle1](const std::pair<int, int>& pair){ return pair.first == particle1 || pair.second == particle1;})
-                                                      != thermostatedPairs.end());
-        bool particle2_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(), thermostatedParticles.end(), particle2)
-                                                    != thermostatedParticles.end())) ||
-                                                  (std::find_if(thermostatedPairs.begin(), thermostatedPairs.end(),
-                                                   [&particle2](const std::pair<int, int>& pair){ return pair.first == particle2 || pair.second == particle2;})
-                                                      != thermostatedPairs.end());
-        if ((system.getParticleMass(particle1) > 0) && (system.getParticleMass(particle2) > 0)){
-            if ((particle1_in_thermostatedParticles && !particle2_in_thermostatedParticles) ||
-                 (!particle1_in_thermostatedParticles && particle2_in_thermostatedParticles)){
-                throw OpenMMException("Cannot add only one of particles " + std::to_string(particle1) + " and " + std::to_string(particle2)
-                                        + " to NoseHooverChain, because they are connected by a constraint.");
+    auto data = ThermostatData(thermostatedParticles, thermostatedPairs, temperature,
+                               relativeTemperature, collisionFrequency,
+                               relativeCollisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
+    thermostatData.push_back(data);
+    return thermostatData.size() - 1;
+}
+
+
+void NoseHooverIntegrator::createThermostats(const System &system) {
+
+    for (const auto &thermostat : thermostatData) {
+
+        // figure out the number of DOFs
+        int nDOF = 3*(thermostat.thermostatedParticles.size() + thermostat.thermostatedPairs.size());
+        for (int constraintNum = 0; constraintNum < system.getNumConstraints(); constraintNum++) {
+            int particle1, particle2;
+            double distance;
+            system.getConstraintParameters(constraintNum, particle1, particle2, distance);
+            bool particle1_in_thermostatedParticles = ((std::find(thermostat.thermostatedParticles.begin(),
+                                                                  thermostat.thermostatedParticles.end(), particle1)
+                                                                    != thermostat.thermostatedParticles.end())) ||
+                                                      (std::find_if(thermostat.thermostatedPairs.begin(),
+                                                                    thermostat.thermostatedPairs.end(),
+                                                                    [&particle1](const std::pair<int, int>& pair){
+                                                                           return pair.first == particle1 || pair.second == particle1;})
+                                                                      != thermostat.thermostatedPairs.end());
+            bool particle2_in_thermostatedParticles = ((std::find(thermostat.thermostatedParticles.begin(),
+                                                                  thermostat.thermostatedParticles.end(), particle2)
+                                                                    != thermostat.thermostatedParticles.end())) ||
+                                                      (std::find_if(thermostat.thermostatedPairs.begin(),
+                                                                    thermostat.thermostatedPairs.end(),
+                                                                    [&particle2](const std::pair<int, int>& pair){
+                                                                           return pair.first == particle2 || pair.second == particle2;})
+                                                                      != thermostat.thermostatedPairs.end());
+            if ((system.getParticleMass(particle1) > 0) && (system.getParticleMass(particle2) > 0)){
+                if ((particle1_in_thermostatedParticles && !particle2_in_thermostatedParticles) ||
+                     (!particle1_in_thermostatedParticles && particle2_in_thermostatedParticles)){
+                    throw OpenMMException("Cannot add only one of particles " + std::to_string(particle1) + " and " + std::to_string(particle2)
+                                            + " to NoseHooverChain, because they are connected by a constraint.");
+                }
+                if (particle1_in_thermostatedParticles && particle2_in_thermostatedParticles){
+                    nDOF -= 1;
+                }
             }
-            if (particle1_in_thermostatedParticles && particle2_in_thermostatedParticles){
-                nDOF -= 1;
+        }
+
+        // remove 3 degrees of freedom from thermostats that act on absolute motions
+        int numForces = system.getNumForces();
+        if (thermostat.thermostatedPairs.size() == 0){
+            for (int forceNum = 0; forceNum < numForces; ++forceNum) {
+                if (dynamic_cast<const CMMotionRemover*>(&system.getForce(forceNum))) nDOF -= 3;
             }
         }
+
+        // create and add new chain
+        int chainID = noseHooverChains.size();
+        auto chain = NoseHooverChain(thermostat.temperature, thermostat.relativeTemperature,
+                                     thermostat.collisionFrequency, thermostat.relativeCollisionFrequency,
+                                     nDOF, thermostat.chainLength, thermostat.numMTS,
+                                     thermostat.numYoshidaSuzuki, chainID,
+                                     thermostat.thermostatedParticles, thermostat.thermostatedPairs);
+        noseHooverChains.push_back(chain);
+
     }
-    int numForces = system.getNumForces();
-    if (thermostatedPairs.size() == 0){ // remove 3 degrees of freedom from thermostats that act on absolute motions
-        for (int forceNum = 0; forceNum < numForces; ++forceNum) {
-            if (dynamic_cast<CMMotionRemover*>(&system.getForce(forceNum))) nDOF -= 3;
+
+
+    for (int chain1 = 0; chain1 < noseHooverChains.size(); ++chain1){
+        const auto& nhc = noseHooverChains[chain1];
+
+       // make sure that thermostats do not overlap
+        for (int chain2 = 0; chain2 < chain1; ++chain2){
+            const auto& other_nhc = noseHooverChains[chain2];
+            for (auto &particle: nhc.getThermostatedAtoms()){
+                bool isParticleInOtherChain = (std::find(other_nhc.getThermostatedAtoms().begin(),
+                                                         other_nhc.getThermostatedAtoms().end(),
+                                                         particle) != other_nhc.getThermostatedAtoms().end()) ||
+                                              (std::find_if(other_nhc.getThermostatedPairs().begin(),
+                                                            other_nhc.getThermostatedPairs().end(),
+                                               [&particle](const std::pair<int, int>& pair){ return pair.first == particle || pair.second == particle;})
+                                                  != other_nhc.getThermostatedPairs().end());
+                if (isParticleInOtherChain){
+                    throw OpenMMException("Found particle " + std::to_string(particle) + "in a different NoseHooverChain, "
+                                          "but particles can only be thermostated by one thermostat.");
+                }
+            }
+            for (auto &pair: nhc.getThermostatedPairs()){
+                bool isParticleInOtherChain = (std::find(other_nhc.getThermostatedAtoms().begin(),
+                                                         other_nhc.getThermostatedAtoms().end(),
+                                                         pair.first) != other_nhc.getThermostatedAtoms().end()) ||
+                                              (std::find(other_nhc.getThermostatedAtoms().begin(),
+                                                         other_nhc.getThermostatedAtoms().end(),
+                                                         pair.second) != other_nhc.getThermostatedAtoms().end()) ||
+                                              (std::find_if(other_nhc.getThermostatedPairs().begin(),
+                                                            other_nhc.getThermostatedPairs().end(),
+                                               [&pair](const std::pair<int, int>& other_pair){
+                                                        return pair.first == other_pair.first || pair.first == other_pair.second ||
+                                                               pair.second == other_pair.first || pair.second == other_pair.second;})
+                                                  != other_nhc.getThermostatedPairs().end());
+                if (isParticleInOtherChain){
+                    throw OpenMMException("Found pair " + std::to_string(pair.first) + "," +
+                                          std::to_string(pair.second) + " in a different NoseHooverChain, "
+                                          "but particles can only be thermostated by one thermostat.");
+                }
+            }
         }
-    }
 
-    // make sure that thermostats do not overlap
-    for (auto &other_nhc: noseHooverChains) {
-        for (auto &particle: thermostatedParticles){
-            bool isParticleInOtherChain = (std::find(other_nhc.getThermostatedAtoms().begin(),
-                                                     other_nhc.getThermostatedAtoms().end(),
-                                                     particle) != other_nhc.getThermostatedAtoms().end()) ||
-                                          (std::find_if(other_nhc.getThermostatedPairs().begin(),
-                                                        other_nhc.getThermostatedPairs().end(),
-                                           [&particle](const std::pair<int, int>& pair){ return pair.first == particle || pair.second == particle;})
-                                              != other_nhc.getThermostatedPairs().end());
-            if (isParticleInOtherChain){
-                throw OpenMMException("Found particle " + std::to_string(particle) + "in a different NoseHooverChain, "
-                                      "but particles can only be thermostated by one thermostat.");
+        // make sure that massless particles are not thermostated
+        for(auto particle: nhc.getThermostatedAtoms()){
+            double mass = system.getParticleMass(particle);
+            if (mass < 1e-8) {
+                throw OpenMMException("Found a particle with no mass (" + std::to_string(particle) + ") in a thermostat. Massless particles cannot be thermostated.");
             }
         }
     }
-
-    // create and add new chain
-    int chainID = noseHooverChains.size();
-    auto nhc = NoseHooverChain(temperature, relativeTemperature, collisionFrequency, relativeCollisionFrequency,
-                               nDOF, chainLength, numMTS, numYoshidaSuzuki, chainID, thermostatedParticles, thermostatedPairs);
-    noseHooverChains.push_back(nhc);
-    return chainID;
 }
 
 double NoseHooverIntegrator::getTemperature(int chainID) const {
@@ -247,6 +291,7 @@ double NoseHooverIntegrator::computeHeatBathEnergy() {
 void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
     if (owner != NULL && &contextRef.getOwner() != owner)
         throw OpenMMException("This Integrator is already bound to a context");
+
     context = &contextRef;
     owner = &contextRef.getOwner();
     vvKernel = context->getPlatform().createKernel(IntegrateVelocityVerletStepKernel::Name(), contextRef);
@@ -254,6 +299,26 @@ void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
     nhcKernel = context->getPlatform().createKernel(NoseHooverChainKernel::Name(), contextRef);
     nhcKernel.getAs<NoseHooverChainKernel>().initialize();
     forcesAreValid = false;
+
+    // check for drude particles and build the Nose-Hoover Chains
+    const System& system = context->getSystem();
+    for (auto& thermostat: thermostatData){
+        // if there are no thermostated particles or pairs in the lists this is a regular thermostat for the whole (non-Drude) system
+        if ( (thermostat.thermostatedParticles.size() == 0) && (thermostat.thermostatedPairs.size() == 0) ){
+            for(int particle = 0; particle < system.getNumParticles(); ++particle) {
+                double mass = system.getParticleMass(particle);
+                if ( (mass > 0) && (mass < 0.8) ){
+                    std::cout << "Warning: Found particles with mass between 0.0 and 0.8 dalton. Did you mean to make a DrudeVelocityVerletIntegrator instead? "
+                                 "The thermostat you are about to use will not treat these particles as Drude particles!" << std::endl;
+                }
+                if(system.getParticleMass(particle) > 0) {
+                    thermostat.thermostatedParticles.push_back(particle);
+                }
+            }
+        }
+    }
+
+    createThermostats(system);
 }
 
 void NoseHooverIntegrator::cleanup() {

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

@@ -55,8 +55,6 @@ public:
      * Create a DrudeNoseHooverIntegrator.
      *
      * @param stepSize       the step size with which to integrator the system (in picoseconds)
-     * @param system the system to be thermostated.  Note: this must be setup, i.e. all
-     *        particles should have been added, before calling this function.
      * @param temperature the target temperature for the system.
      * @param drudeTemperature the target temperature for the Drude particles, relative to their parent atom.
      * @param collisionFrequency the frequency of the system's interaction with the heat bath (in 1/ps).
@@ -66,7 +64,7 @@ public:
      * @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
      *        used in the chain propagation algorithm (must be 1, 3, or 5).
      */
-    DrudeNoseHooverIntegrator(double stepSize, System &system, double temperature, double drudeTemperature,
+    DrudeNoseHooverIntegrator(double stepSize, double temperature, double drudeTemperature,
                               double collisionFrequency, double drudeCollisionFrequency,
                               int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 3);
 

plugins/drude/openmmapi/src/DrudeNoseHooverIntegrator.cpp

@@ -46,49 +46,36 @@ using namespace OpenMM;
 using std::string;
 using std::vector;
 
-DrudeNoseHooverIntegrator::DrudeNoseHooverIntegrator(double stepSize, System &system, double temperature, double drudeTemperature,
+DrudeNoseHooverIntegrator::DrudeNoseHooverIntegrator(double stepSize, double temperature, double drudeTemperature,
                                                      double collisionFrequency, double drudeCollisionFrequency,
                                                      int chainLength, int numMTS, int numYoshidaSuzuki) :
     NoseHooverIntegrator(stepSize) {
-    const DrudeForce* drudeForce = NULL;
-    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 (drudeForce == NULL)
-        throw OpenMMException("The System does not contain a DrudeForce");
-    std::set<int> realParticlesSet;
-    vector<int> realParticles, drudeParticles, drudeParents;
-    vector<std::pair<int, int>> drudePairs;
-    for (int i = 0; i < system.getNumParticles(); 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);
-        drudePairs.push_back({p,p1});
-    }
-    for(const auto &p : realParticlesSet) realParticles.push_back(p);
 
-    addSubsystemThermostat(system, realParticles, drudePairs, temperature, drudeTemperature, collisionFrequency,
-                           drudeCollisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
+    addSubsystemThermostat(std::vector<int>(), std::vector<std::pair<int, int>>(), temperature,
+                           drudeTemperature, collisionFrequency, drudeCollisionFrequency,
+                           chainLength, numMTS, numYoshidaSuzuki);
 }
 
 DrudeNoseHooverIntegrator::~DrudeNoseHooverIntegrator() { }
 
 
 void DrudeNoseHooverIntegrator::initialize(ContextImpl& contextRef) {
-    NoseHooverIntegrator::initialize(contextRef);
+
     if (owner != NULL && &contextRef.getOwner() != owner)
         throw OpenMMException("This Integrator is already bound to a context");
+
+    context = &contextRef;
+    owner = &contextRef.getOwner();
+    vvKernel = context->getPlatform().createKernel(IntegrateVelocityVerletStepKernel::Name(), contextRef);
+    vvKernel.getAs<IntegrateVelocityVerletStepKernel>().initialize(contextRef.getSystem(), *this);
+    nhcKernel = context->getPlatform().createKernel(NoseHooverChainKernel::Name(), contextRef);
+    nhcKernel.getAs<NoseHooverChainKernel>().initialize();
+    forcesAreValid = false;
+
+    // check for drude particles and build the Nose-Hoover Chains
+    const System& system = context->getSystem();
     const DrudeForce* drudeForce = NULL;
-    const System& system = contextRef.getSystem();
+
     bool hasCMMotionRemover = false;
     for (int i = 0; i < system.getNumForces(); i++){
         if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
@@ -107,8 +94,27 @@ void DrudeNoseHooverIntegrator::initialize(ContextImpl& contextRef) {
         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();
+    for (auto& thermostat: thermostatData){
+        if ( (thermostat.thermostatedParticles.size() == 0) && (thermostat.thermostatedPairs.size() == 0) ){
+            std::set<int> realParticlesSet;
+            vector<int> realParticles, drudeParticles, drudeParents;
+            vector<std::pair<int, int>> drudePairs;
+            for (int i = 0; i < system.getNumParticles(); 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);
+                thermostat.thermostatedPairs.push_back({p,p1});
+            }
+            for(const auto &p : realParticlesSet) thermostat.thermostatedParticles.push_back(p);
+        }
+    }
+
+    createThermostats(system);
 }
 
 double DrudeNoseHooverIntegrator::computeDrudeKineticEnergy() {

plugins/drude/tests/TestDrudeNoseHoover.h

@@ -120,7 +120,7 @@ void testWaterBox(Platform& platform){
     double frequency = 100.0;
     double frequencyDrude = 80.0;
     int randomSeed = 100;
-    DrudeNoseHooverIntegrator integ(0.0005, system, temperature, temperatureDrude,
+    DrudeNoseHooverIntegrator integ(0.0005, temperature, temperatureDrude,
                                     frequency, frequencyDrude, chainLength, numMTS, numYS);;
     Context context(system, integ, platform);
     context.setPositions(positions);

tests/TestNoseHooverThermostat.h

@@ -66,7 +66,7 @@ void testHarmonicOscillator() {
     harmonic_restraint->addParticle(0);
     system.addForce(harmonic_restraint);
     NoseHooverIntegrator integrator(0.001);
-    integrator.addThermostat(system, temperature, frequency, chain_length, mts, ys);
+    integrator.addThermostat(temperature, frequency, chain_length, mts, ys);
     Context context(system, integrator, platform);
     context.setPositions(positions);
     context.setVelocities(velocities);
@@ -162,10 +162,10 @@ void testDimerBox(bool constrain=true) {
     int numMTS = 3;
     int numYS = 3;
     int chainLength = 5;
-    auto integrator = simpleConstruct ? NoseHooverIntegrator(0.001, system, temperature, collisionFrequency, chainLength, numMTS, numYS)
+    auto integrator = simpleConstruct ? NoseHooverIntegrator(0.001, temperature, collisionFrequency, chainLength, numMTS, numYS)
                                       : NoseHooverIntegrator(0.001);
     if (!simpleConstruct)
-        integrator.addThermostat(system, temperature, collisionFrequency, chainLength, numMTS, numYS);
+        integrator.addThermostat(temperature, collisionFrequency, chainLength, numMTS, numYS);
     Context context(system, integrator, platform);
     context.setPositions(positions);
     context.setVelocitiesToTemperature(temperature);
@@ -225,9 +225,9 @@ void testCheckpoints() {
     double kineticEnergy = 1e6;
     double temperature=300, collisionFrequency=1, chainLength=3, numMTS=3, numYS=3;
     chainLength = 10;
-    integrator.addSubsystemThermostat(system, std::vector<int>(), std::vector<std::pair<int,int>>{{0,1}},  temperature, temperature, collisionFrequency, collisionFrequency,
+    integrator.addSubsystemThermostat(std::vector<int>(), std::vector<std::pair<int,int>>{{0,1}},  temperature, temperature, collisionFrequency, collisionFrequency,
                                       chainLength, numMTS, numYS);
-    newIntegrator.addSubsystemThermostat(system, std::vector<int>(), std::vector<std::pair<int,int>>{{0,1}},  temperature, temperature, collisionFrequency, collisionFrequency,
+    newIntegrator.addSubsystemThermostat(std::vector<int>(), std::vector<std::pair<int,int>>{{0,1}},  temperature, temperature, collisionFrequency, collisionFrequency,
                                       chainLength, numMTS, numYS);
     Context context(system, integrator, platform);
     Context newContext(system, newIntegrator, platform);