Initial infrastructure for DrudeVelocityVerletIntegrator

openmmapi/include/openmm/NoseHooverChain.h

@@ -260,13 +260,26 @@ public:
 protected:
     ForceImpl* createImpl() const;
 private:
-    double defaultTemp, defaultFreq, defaultTimeStep;
+    double defaultTemp, defaultFreq; //, defaultTimeStep;
     int defaultNumDOFs, defaultChainLength, defaultNumMTS, defaultNumYS;
     // The suffix used to distinguish NH chains, e.g. for Drude particles vs. regular particles.
     int defaultChainID;
     std::vector<int> thermostatedAtoms, parentAtoms;
 };
 
+/**
+ * Check if two Nose-Hoover chains are identical (checks all member variables but the chain id).
+ */
+//inline bool operator==(const NoseHooverChain& lhs, const NoseHooverChain& rhs){
+//    if (lhs.getDefaultTemperature() != rhs.getDefaultTemperature()) return false;
+//    if (lhs.getDefaultCollisionFrequency() != rhs.getDefaultCollisionFrequency()) return false;
+//    if (lhs.getDefaultNumDegreesOfFreedom() != rhs.getDefaultNumDegreesOfFreedom()) return false;
+//    if (lhs.getDefaultNumMultiTimeSteps() != rhs.getDefaultNumMultiTimeSteps()) return false;
+//    if (lhs.getDefaultNumYoshidaSuzukiTimeSteps() != rhs.getDefaultNumYoshidaSuzukiTimeSteps()) return false;
+//    if (lhs.getThermostatedAtoms() != rhs.getThermostatedAtoms()) return false;
+//    if (lhs.getParentAtoms() != rhs.getParentAtoms()) return false;
+//}
+
 } // namespace OpenMM
 
 #endif /*OPENMM_NOSEHOOVERCHAIN_H_*/

openmmapi/include/openmm/VelocityVerletIntegrator.h

@@ -54,6 +54,8 @@ public:
      * @param stepSize the step size with which to integrate the system (in picoseconds)
      */
     explicit VelocityVerletIntegrator(double stepSize);
+
+    virtual ~VelocityVerletIntegrator(); 
    /**
      * Advance a simulation through time by taking a series of time steps.
      * 
@@ -68,7 +70,7 @@ public:
      */
     double propagateChain(double kineticEnergy, int chainID=0);
     /**
-     * Add a Nose-Hoover Chain thermostat to control the temperature of the systeml
+     * Add a Nose-Hoover Chain thermostat to control the temperature of the 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.
@@ -81,6 +83,54 @@ public:
      */
      int addNoseHooverChainThermostat(System& system, double temperature, double collisionFrequency,
                                              int chainLength, int numMTS, int numYoshidaSuzuki);
+    /**
+     * Add a Nose-Hoover Chain thermostat to control the temperature of the 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 mask list of particle ids to be thermostated.
+     * @param parents either an empty list of a list describing the parent atoms that each thermostated
+     *        atom is connected to.
+     * @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.
+     * @param numMTS the number of step in the  multiple time step chain propagation algorithm.
+     * @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 addMaskedNoseHooverChainThermostat(System& system, const std::vector<int>& mask, const std::vector<int>& parents,
+                                             double temperature, double collisionFrequency,
+                                             int chainLength, int numMTS, int numYoshidaSuzuki);
+    /**
+     * Get the temperature of the i-th chain (in Kelvin).
+     * 
+     * @param chainID the index of the Nose-Hoover chain (default=0).
+     * 
+     * @return the temperature.
+     */
+    double getTemperature(int chainID=0) const;
+    /**
+     * set the temperature of the i-th chain.
+     *
+     * @param temperature the temperature for the Nose-Hoover chain thermostat (in Kelvin).
+     * @param chainID The id of the Nose-Hoover chain for which the temperature is set (default=0).
+     */
+    void setTemperature(double temperature, int chainID=0);
+    /**
+     * Get the collision frequency of the i-th chain (in 1/picosecond).
+     * 
+     * @param chainID the index of the Nose-Hoover chain (default=0).
+     *
+     * @return the collision frequency.
+     */
+    double getCollisionFrequency(int chainID=0) const;
+    /**
+     * Set the collision frequency of the i-th chain.
+     *
+     * @param frequency the collision frequency in picosecond.
+     * @param chainID the index of the Nose-Hoover chain (default=0).
+     */
+    void setCollisionFrequency(double frequency, int chainID=0);
     /**
      * Compute the total (potential + kinetic) heat bath energy for all heat baths
      * associated with this integrator, at the current time.
@@ -107,9 +157,9 @@ protected:
      * Compute the kinetic energy of the system at the current time.
      */
     double computeKineticEnergy();
+    std::map<int, NoseHooverChain*> noseHooverChains;
 private:
     Kernel vvKernel, nhcKernel;
-    std::vector<NoseHooverChain*> noseHooverChains;
 };
 
 } // namespace OpenMM

openmmapi/src/VelocityVerletIntegrator.cpp

@@ -39,6 +39,7 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/kernels.h"
 #include <string>
+#include <algorithm>
 #include <iostream>
 
 using namespace OpenMM;
@@ -50,51 +51,151 @@ VelocityVerletIntegrator::VelocityVerletIntegrator(double stepSize) {
     setConstraintTolerance(1e-5);
 }
 
+VelocityVerletIntegrator::~VelocityVerletIntegrator() {}
+
 double VelocityVerletIntegrator::propagateChain(double kineticEnergy, int chainID) {
-    return nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, *noseHooverChains[chainID], kineticEnergy, getStepSize());
+    return nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, *noseHooverChains.at(chainID), kineticEnergy, getStepSize());
 }
 
-/*int VelocityVerletIntegrator::addMaskedNoseHooverChain(System& system, std::vector<int> mask, std::vector<int> parents, double temperature,
-                                                 double frequency, int numDOFs, int numMTS, int numYoshidaSuzuki){
-
-}*/
-
 int VelocityVerletIntegrator::addNoseHooverChainThermostat(System& system, double temperature, double collisionFrequency,
                                                            int chainLength, int numMTS, int numYoshidaSuzuki) {
-    if (context) {
-        throw OpenMMException("addNoseHooverChainThermostat cannot be called after binding this integrator to a context.");
-    }
-    int nDOF = 0;
     int numForces = system.getNumForces();
     std::vector<int> thermostatedParticles;
     std::vector<int> parentParticles;
     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) {
-            nDOF += 3;
             thermostatedParticles.push_back(particle);
         }
     }
-    nDOF -= system.getNumConstraints();
+
+    return addMaskedNoseHooverChainThermostat(system, thermostatedParticles, parentParticles, temperature, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
+}
+
+int VelocityVerletIntegrator::addMaskedNoseHooverChainThermostat(System& system, const std::vector<int>& thermostatedParticles, const std::vector<int>& parentParticles,
+                                         double temperature, double collisionFrequency,
+                                         int chainLength, int numMTS, int numYoshidaSuzuki){
+    const auto & mask = thermostatedParticles; // just an alias
+    if (context) {
+        throw OpenMMException("Nose-Hoover chains cannot be added after binding this integrator to a context.");
+    }
+    if ((parentParticles.size() != mask.size()) && (parentParticles.size() != 0)) {
+        throw OpenMMException("The number of parent particles has to be either equal to the number of thermostated particles (to thermostat relative motion)"
+                              " or zero (to thermostat absolute motion).");
+    }
+    // figure out the number of DOFs
+    int nDOF = 3*mask.size();
+    for (int constraintNum = 0; constraintNum < system.getNumConstraints(); constraintNum++) {
+        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());
+        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)
+                                        + " to NoseHooverChain, because they are connected by a constraint.");
+            }
+            if (particle1_in_mask && particle2_in_mask){
+                nDOF -= 1;
+            }
+        }
+    }
+    int numForces = system.getNumForces();
+    if (parentParticles.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;
         }
+    }
 
+    // figure out the chain ID
+    std::vector<int> takenIDs;
+    for (int forceNum = 0; forceNum < numForces; ++forceNum){
+        const auto & other_nhc = dynamic_cast<NoseHooverChain*>(&system.getForce(forceNum));
+        if (other_nhc){
+            takenIDs.push_back(other_nhc->getDefaultChainID());
+        }
+    }
+    int chainID = 0;
+    while (true) {
+        bool isChainIDTaken = (std::find(takenIDs.begin(), takenIDs.end(), chainID) != takenIDs.end());
+        bool isChainIDRegistered = (noseHooverChains.find(chainID) != noseHooverChains.end());
+        if (isChainIDTaken || isChainIDRegistered) {
+            chainID ++; 
+        } else {
+            break;
+        }
+    }
     auto nhcForce = new NoseHooverChain(temperature, collisionFrequency, nDOF, chainLength,
-                                        numMTS, numYoshidaSuzuki, noseHooverChains.size(),
+                                        numMTS, numYoshidaSuzuki, chainID,
                                         thermostatedParticles, parentParticles);
+
+    // make sure that thermostats do not overlap 
+    for (auto &other_nhc: noseHooverChains) {
+        for (auto &particle: mask){
+            bool isParticleInOtherChain = (std::find(other_nhc.second->getThermostatedAtoms().begin(), 
+                                                    other_nhc.second->getThermostatedAtoms().end(),
+                                                    particle) == other_nhc.second->getThermostatedAtoms().begin());
+            if (isParticleInOtherChain){
+                throw OpenMMException("Found particle " + std::to_string(particle) + "in a different NoseHooverChain, "
+                                      "but particles can only be thermostated by one thermostat.");
+            }
+        }
+    }
     system.addForce(nhcForce);
-    noseHooverChains.push_back(nhcForce);
-    return noseHooverChains.size() - 1;
+    noseHooverChains[chainID] = nhcForce;
+    return chainID;
+}
+
+double VelocityVerletIntegrator::getTemperature(int chainID) const {
+    if (chainID >= noseHooverChains.size()) {
+        throw OpenMMException("Cannot get temperature for chainID " + std::to_string(chainID)
+                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
+        );
+    }
+    return noseHooverChains.at(chainID)->getDefaultTemperature();
+}
+
+void VelocityVerletIntegrator::setTemperature(double temperature, int chainID){
+    if (chainID >= noseHooverChains.size()) {
+        throw OpenMMException("Cannot set temperature for chainID " + std::to_string(chainID)
+                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
+        );
+    }
+    noseHooverChains.at(chainID)->setDefaultTemperature(temperature);
+
+}
+
+double VelocityVerletIntegrator::getCollisionFrequency(int chainID) const {
+    if (chainID >= noseHooverChains.size()) {
+        throw OpenMMException("Cannot get collision frequency for chainID " + std::to_string(chainID)
+                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
+        );
+    }
+    return noseHooverChains.at(chainID)->getDefaultCollisionFrequency();
+} 
+void VelocityVerletIntegrator::setCollisionFrequency(double frequency, int chainID){
+    if (chainID >= noseHooverChains.size()) {
+        throw OpenMMException("Cannot set collision frequency for chainID " + std::to_string(chainID)
+                + ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
+        );
+    }
+    noseHooverChains.at(chainID)->setDefaultCollisionFrequency(frequency);
 }
 
 double VelocityVerletIntegrator::computeKineticEnergy() {
     return vvKernel.getAs<IntegrateVelocityVerletStepKernel>().computeKineticEnergy(*context, *this);
+
 }
 
 double VelocityVerletIntegrator::computeHeatBathEnergy() {
     double energy = 0;
     for(auto &nhc : noseHooverChains) {
-        energy += nhcKernel.getAs<NoseHooverChainKernel>().computeHeatBathEnergy(*context, *nhc);
+        energy += nhcKernel.getAs<NoseHooverChainKernel>().computeHeatBathEnergy(*context, *nhc.second);
     }
     return energy;
 }
@@ -129,15 +230,15 @@ 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);
-            scale = nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, *nhc, kineticEnergy, getStepSize());
-            nhcKernel.getAs<NoseHooverChainKernel>().scaleVelocities(*context, *nhc, scale);
+            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, *nhc.second);
+            scale = nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, *nhc.second, kineticEnergy, getStepSize());
+            nhcKernel.getAs<NoseHooverChainKernel>().scaleVelocities(*context, *nhc.second, scale);
         }
         vvKernel.getAs<IntegrateVelocityVerletStepKernel>().execute(*context, *this);
         for(auto &nhc : noseHooverChains) {
-            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, *nhc);
-            scale = nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, *nhc, kineticEnergy, getStepSize());
-            nhcKernel.getAs<NoseHooverChainKernel>().scaleVelocities(*context, *nhc, scale);
+            kineticEnergy = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, *nhc.second);
+            scale = nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, *nhc.second, kineticEnergy, getStepSize());
+            nhcKernel.getAs<NoseHooverChainKernel>().scaleVelocities(*context, *nhc.second, scale);
         }
     }
 }

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

+#ifndef OPENMM_DRUDEVELOCITYVERLETINTEGRATOR_H_
+#define OPENMM_DRUDEVELOCITYVERLETINTEGRATOR_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   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) 2019 Stanford University and the Authors.           *
+ * Authors: Andreas Krämer and Andrew C. Simmonett                            *
+ * 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 "openmm/VelocityVerletIntegrator.h"
+#include "openmm/Kernel.h"
+#include "openmm/internal/windowsExportDrude.h"
+
+namespace OpenMM {
+
+/**
+ * This Integrator simulates systems that include Drude particles.  It applies two different Nose-Hoover
+ * chain thermostats to the different parts of the system.  The first is applied to ordinary particles (ones
+ * that are not part of a Drude particle pair), as well as to the center of mass of each Drude particle pair.
+ * A second thermostat, typically with a much lower temperature, is applied to the relative internal
+ * displacement of each pair.
+ *
+ * This Integrator requires the System to include a DrudeForce, which it uses to identify the Drude
+ * particles.
+ */
+
+class OPENMM_EXPORT_DRUDE DrudeVelocityVerletIntegrator : public VelocityVerletIntegrator {
+public:
+    /**
+     * Create a DrudeVelocityVerletIntegrator.
+     *
+     * @param stepSize       the step size with which to integrator the system (in picoseconds)
+     */
+    DrudeVelocityVerletIntegrator(double stepSize);
+
+    virtual ~DrudeVelocityVerletIntegrator(); 
+    /**
+     * Add a Nose-Hoover Chain thermostat to control the temperature of the 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 temperature the target temperature for the drude particles..
+     * @param collisionFrequency the frequency of the interaction of the drude particles with the heat bath (in 1/ps).
+     * @param chainLength the number of beads in the Nose-Hoover chain.
+     * @param numMTS the number of step in the  multiple time step chain propagation algorithm.
+     * @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 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.
+     * It will also get called again if the application calls reinitialize() on the Context.
+     */
+    void initialize(ContextImpl& context);
+    /**
+     * Compute the kinetic energy of the system at the current time.
+     */
+    //double computeKineticEnergy();
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_DRUDEVELOCITYVERLETINTEGRATOR_H_*/

plugins/drude/openmmapi/src/DrudeVelocityVerletIntegrator.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) 2019 Stanford University and the Authors.           *
+ * Authors: Andreas Krämer and Andrew C. Simmonett                            *
+ * 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 "openmm/DrudeVelocityVerletIntegrator.h"
+#include "openmm/Context.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/DrudeKernels.h"
+#include "openmm/DrudeForce.h"
+#include <ctime>
+#include <string>
+#include <set>
+
+using namespace OpenMM;
+using std::string;
+using std::vector;
+
+DrudeVelocityVerletIntegrator::DrudeVelocityVerletIntegrator(double stepSize) : VelocityVerletIntegrator(stepSize) { }
+
+DrudeVelocityVerletIntegrator::~DrudeVelocityVerletIntegrator() { }
+
+int DrudeVelocityVerletIntegrator::addDrudeNoseHooverChainThermostat(System& system, double temperature, double collisionFrequency,
+                                                                     double drudeTemperature, double drudeCollisionFrequency,
+                                                                     int chainLength, int numMTS, int numYoshidaSuzuki) {
+    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;
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        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);
+    }
+    for(const auto &p : realParticlesSet) realParticles.push_back(p);
+
+    addMaskedNoseHooverChainThermostat(system, realParticles, vector<int>(), temperature, collisionFrequency,
+                                      chainLength, numMTS, numYoshidaSuzuki);
+    addMaskedNoseHooverChainThermostat(system, drudeParticles, drudeParents, drudeTemperature, drudeCollisionFrequency,
+                                      chainLength, numMTS, numYoshidaSuzuki);
+    return noseHooverChains.size() - 1;
+}
+
+void DrudeVelocityVerletIntegrator::initialize(ContextImpl& 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++)
+        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");
+    context = &contextRef;
+    owner = &contextRef.getOwner();
+}