Clean up NHC interface and fix bug in swig generator

openmmapi/include/openmm/NoseHooverIntegrator.h

@@ -58,15 +58,15 @@ public:
     /**
      * Create a NoseHooverIntegrator.
      *
-     * @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 system (in Kelvin).
+     * @param collisionFrequency the frequency of the interaction with the heat bath (in inverse picoseconds).
      * @param stepSize the step size with which to integrate the system (in picoseconds)
      * @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).
      */
-    explicit NoseHooverIntegrator(double temperature, int collisionFrequnency, double stepSize, 
+    explicit NoseHooverIntegrator(double temperature, double collisionFrequency, double stepSize,
                                   int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 3);
 
     virtual ~NoseHooverIntegrator();
@@ -81,7 +81,7 @@ public:
      *
      * @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 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).
@@ -94,6 +94,7 @@ public:
      * provided and the thermostat will only control members of that list.  Additionally a list of pairs of
      * 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.
+     * If both the list of thermostated particles and thermostated pairs are empty all particles will be thermostated.
      *
      * @param thermostatedParticles list of particle ids to be thermostated.
      * @param thermostatedPairs a list of pairs of connected atoms whose absolute center of mass motion
@@ -215,9 +216,9 @@ protected:
      */
     void initialize(ContextImpl& context);
     /**
-     * Goes through the list of requested thermostat data and creates the thermostat chains.
+     * Goes through the list of thermostats, sets the number of DOFs, and checks for errors in the thermostats.
      */
-    void createThermostats(const System& system);
+    void initializeThermostats(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.
@@ -232,25 +233,6 @@ 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

@@ -53,11 +53,11 @@ NoseHooverIntegrator::NoseHooverIntegrator(double stepSize):
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
 }
-NoseHooverIntegrator::NoseHooverIntegrator(double temperature, int collisionFrequnency, double stepSize,
+NoseHooverIntegrator::NoseHooverIntegrator(double temperature, double collisionFrequency, double stepSize,
                                            int chainLength, int numMTS, int numYoshidaSuzuki) : forcesAreValid(false) {
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
-    addThermostat(temperature, collisionFrequnency, chainLength, numMTS, numYoshidaSuzuki);
+    addThermostat(temperature, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
 }
 
 NoseHooverIntegrator::~NoseHooverIntegrator() {}
@@ -79,40 +79,45 @@ int NoseHooverIntegrator::addSubsystemThermostat(const std::vector<int>& thermos
                                                  double temperature, double collisionFrequency,
                                                  double relativeTemperature, double relativeCollisionFrequency,
                                                  int chainLength, int numMTS, int numYoshidaSuzuki) {
-    auto data = ThermostatData(thermostatedParticles, thermostatedPairs, temperature,
-                               relativeTemperature, collisionFrequency,
-                               relativeCollisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
-    thermostatData.push_back(data);
-    return thermostatData.size() - 1;
+    int chainID = noseHooverChains.size();
+    int nDOF = 0; // is set in initializeThermostats()
+    noseHooverChains.emplace_back(temperature, relativeTemperature,
+                                 collisionFrequency, relativeCollisionFrequency,
+                                 nDOF, chainLength, numMTS,
+                                 numYoshidaSuzuki, chainID,
+                                 thermostatedParticles, thermostatedPairs);
+    return chainID;
 }
 
 
-void NoseHooverIntegrator::createThermostats(const System &system) {
+void NoseHooverIntegrator::initializeThermostats(const System &system) {
 
-    for (const auto &thermostat : thermostatData) {
+    for (auto &thermostat : noseHooverChains) {
+        const auto &thermostatedParticles = thermostat.getThermostatedAtoms();
+        const auto &thermostatedPairs = thermostat.getThermostatedPairs();
 
         // figure out the number of DOFs
-        int nDOF = 3*(thermostat.thermostatedParticles.size() + thermostat.thermostatedPairs.size());
+        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(thermostat.thermostatedParticles.begin(),
-                                                                  thermostat.thermostatedParticles.end(), particle1)
-                                                                    != thermostat.thermostatedParticles.end())) ||
-                                                      (std::find_if(thermostat.thermostatedPairs.begin(),
-                                                                    thermostat.thermostatedPairs.end(),
+            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;})
-                                                                      != 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(),
+                                                                      != 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;})
-                                                                      != thermostat.thermostatedPairs.end());
+                                                                      != thermostatedPairs.end());
             if ((system.getParticleMass(particle1) > 0) && (system.getParticleMass(particle2) > 0)){
                 if ((particle1_in_thermostatedParticles && !particle2_in_thermostatedParticles) ||
                      (!particle1_in_thermostatedParticles && particle2_in_thermostatedParticles)){
@@ -127,21 +132,14 @@ void NoseHooverIntegrator::createThermostats(const System &system) {
 
         // remove 3 degrees of freedom from thermostats that act on absolute motions
         int numForces = system.getNumForces();
-        if (thermostat.thermostatedPairs.size() == 0){
+        if (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);
-
+        // set number of DoFs for chain 
+        thermostat.setDefaultNumDegreesOfFreedom(nDOF);
     }
 
 
@@ -271,9 +269,10 @@ void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
 
     // check for drude particles and build the Nose-Hoover Chains
     const System& system = context->getSystem();
-    for (auto& thermostat: thermostatData){
+    for (auto& thermostat: noseHooverChains){
         // 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) ){
+        if ( (thermostat.getThermostatedAtoms().size() == 0) && (thermostat.getThermostatedPairs().size() == 0) ){
+            std::vector<int> thermostatedParticles;
             for(int particle = 0; particle < system.getNumParticles(); ++particle) {
                 double mass = system.getParticleMass(particle);
                 if ( (mass > 0) && (mass < 0.8) ){
@@ -281,13 +280,14 @@ void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
                                  "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);
+                   thermostatedParticles.push_back(particle);
                 }
             }
+            thermostat.setThermostatedAtoms(thermostatedParticles);
         }
     }
 
-    createThermostats(system);
+    initializeThermostats(system);
 }
 
 void NoseHooverIntegrator::cleanup() {

plugins/drude/openmmapi/include/OpenMMDrude.h

@@ -35,6 +35,6 @@
 #include "openmm/DrudeForce.h"
 #include "openmm/DrudeLangevinIntegrator.h"
 #include "openmm/DrudeSCFIntegrator.h"
-#include "openmm/DrudeNIntegrator.h"
+#include "openmm/DrudeNoseHooverIntegrator.h"
 
 #endif /*OPENMM_DRUDE_H_*/

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

@@ -54,10 +54,10 @@ public:
     /**
      * Create a DrudeNoseHooverIntegrator.
      *
-     * @param temperature the target temperature for the system.
-     * @param collisionFrequency the frequency of the system's interaction with the heat bath (in 1/ps).
-     * @param drudeTemperature the target temperature for the Drude particles, relative to their parent atom.
-     * @param drudeCollisionFrequency the frequency of the drude particles' interaction with the heat bath (in 1/ps).
+     * @param temperature the target temperature for the system (in Kelvin).
+     * @param collisionFrequency the frequency of the system's interaction with the heat bath (in inverse picoseconds).
+     * @param drudeTemperature the target temperature for the Drude particles, relative to their parent atom (in Kelvin).
+     * @param drudeCollisionFrequency the frequency of the drude particles' interaction with the heat bath (in inverse picoseconds).
      * @param stepSize       the step size with which to integrator the system (in picoseconds)
      * @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.

plugins/drude/openmmapi/src/DrudeNoseHooverIntegrator.cpp

@@ -94,8 +94,8 @@ 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;
     }
-    for (auto& thermostat: thermostatData){
-        if ( (thermostat.thermostatedParticles.size() == 0) && (thermostat.thermostatedPairs.size() == 0) ){
+    for (auto& thermostat: noseHooverChains){
+        if ( (thermostat.getThermostatedAtoms().size() == 0) && (thermostat.getThermostatedPairs().size() == 0) ){
             std::set<int> realParticlesSet;
             vector<int> realParticles, drudeParticles, drudeParents;
             vector<std::pair<int, int>> drudePairs;
@@ -108,13 +108,14 @@ void DrudeNoseHooverIntegrator::initialize(ContextImpl& contextRef) {
                 drudeForce->getParticleParameters(i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
                 realParticlesSet.erase(p);
                 realParticlesSet.erase(p1);
-                thermostat.thermostatedPairs.push_back({p,p1});
+                drudePairs.push_back({p,p1});
             }
-            for(const auto &p : realParticlesSet) thermostat.thermostatedParticles.push_back(p);
+            thermostat.setThermostatedPairs(drudePairs);
+            thermostat.setThermostatedAtoms(std::vector<int>(realParticlesSet.begin(), realParticlesSet.end()));
         }
     }
 
-    createThermostats(system);
+    initializeThermostats(system);
 }
 
 double DrudeNoseHooverIntegrator::computeDrudeKineticEnergy() {

wrappers/python/src/swig_doxygen/swigInputBuilder.py

@@ -21,6 +21,9 @@ except ImportError:
 INDENT = "   "
 docTags = {'emphasis':'i', 'bold':'b', 'itemizedlist':'ul', 'listitem':'li', 'preformatted':'pre', 'computeroutput':'tt', 'subscript':'sub'}
 
+def is_method_abstract(argstring):
+    return argstring.split(")")[-1].find("=0") >= 0
+
 def striphtmltags(s):
     """Strip a couple html tags used inside docstrings in the C++ source
     to produce something more easily read as plain text.
@@ -511,7 +514,7 @@ class SwigInputBuilder:
                 mArgsstring = getText("argsstring", memberNode)
             if self.fOutPythonprepend and \
                len(paramList) and \
-               mArgsstring.find('=0') < 0:
+               not is_method_abstract(mArgsstring):
                 text = '''
 %pythonprepend OpenMM::{shortClassName}::{methName}{mArgsstring} %{{{{{{0}}
 %}}}}'''.format(shortClassName=shortClassName, methName=methName, mArgsstring=mArgsstring)
@@ -529,6 +532,7 @@ class SwigInputBuilder:
              % self.__class__.__name__)
         raise Exception(s)'''.format(argName=argName)
 
+
                 # Convert input arguments to the proper units, if specified.
                 if key not in methodsWithOutputArgs:
                     if key in self.configModule.UNITS:
@@ -563,9 +567,10 @@ class SwigInputBuilder:
 
             # write pythonappend blocks
             if self.fOutPythonappend \
-               and mArgsstring.find('=0') < 0:
+               and not is_method_abstract(mArgsstring):
                 key = (shortClassName, methName)
-                #print "key %s %s \n" % (shortClassName, methName)
+                #sys.stdout.write("key %s %s \n" % (shortClassName, methName))
+
                 addText=''
                 returnType = getText("type", memberNode)
 
@@ -619,8 +624,15 @@ class SwigInputBuilder:
                             pType = getText('type', pNode)
                         except IndexError:
                             pType = getText('type/ref', pNode)
+                        # parse default arguments
+                        try:
+                            defaultValue = getText('defval', pNode)
+                        except:
+                            defaultValue = ""
+                        if defaultValue != "":
+                            defaultValue = "=%s" %defaultValue
                         pName = getText('declname', pNode)
-                        self.fOutPythonappend.write("%s%s %s" % (sepChar, pType, pName))
+                        self.fOutPythonappend.write("%s%s %s%s" % (sepChar, pType, pName, defaultValue))
                         sepChar=', '
 
                         if pType.find('&')>=0 and \

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -217,6 +217,9 @@ UNITS = {
 ("*", "getTabulatedFunction") : (None, ()),
 ("*", "getUseDispersionCorrection") : (None, ()),
 ("*", "getTemperature") : ("unit.kelvin", ()),
+("*", "getRelativeTemperature") : ("unit.kelvin", ()),
+("*", "getCollisionFrequency") : ( "1/unit.picosecond", ()),
+("*", "getRelativeCollisionFrequency") : ( "1/unit.picosecond", ()),
 ("*", "getUseDispersionCorrection") : (None, ()),
 ("*", "getWeight") : (None, ()),
 ("*", "getWeight12") : (None, ()),
@@ -343,9 +346,11 @@ UNITS = {
 ("HippoNonbondedForce",                 "getInducedDipoles")                             :  ( None, ()),
 ("HippoNonbondedForce",                 "getLabFramePermanentDipoles")                   :  ( None, ()),
 
+
 ("Context", "getParameter") : (None, ()),
 ("Context", "getParameters") : (None, ()),
 ("Context", "getMolecules") : (None, ()),
+("Context", "getState") : (None, (None, None, None)),
 ("CMAPTorsionForce", "getMapParameters") : (None, (None, 'unit.kilojoule_per_mole')),
 ("CMAPTorsionForce", "getTorsionParameters") : (None, ()),
 ("CMMotionRemover", "getFrequency") : (None, ()),
@@ -466,17 +471,28 @@ UNITS = {
 ("System", "getVirtualSite") : (None, ()),
 ("DrudeLangevinIntegrator", "getDrudeTemperature") : ("unit.kelvin", ()),
 ("DrudeLangevinIntegrator", "getMaxDrudeDistance") : ("unit.nanometer", ()),
+
 ("MonteCarloMembraneBarostat", "MonteCarloMembraneBarostat") : (None, ("unit.bar", "unit.bar*unit.nanometer", "unit.kelvin", None, None, None)),
 ("MonteCarloMembraneBarostat", "getXYMode") : (None, ()),
 ("MonteCarloMembraneBarostat", "getZMode") : (None, ()),
 ("DrudeLangevinIntegrator", "getDrudeFriction") : ("1/unit.picosecond", ()),
-("NoseHooverChain", "getThermostatedPairs") : (None, ()),
-("NoseHooverChain", "getThermostatedAtoms") : (None, ()),
-("NoseHooverChain", "getDefaultYoshidaSuzukiWeights") : (None, ()),
 ("DrudeSCFIntegrator", "getMinimizationErrorTolerance") : ("unit.kilojoules_per_mole/unit.nanometer", ()),
 ("RPMDIntegrator", "getContractions") : (None, ()),
 ("RPMDIntegrator", "getTotalEnergy") : ("unit.kilojoules_per_mole", ()),
+("RPMDIntegrator", "getState"): (None,(None, None, None, None)),
 ("RMSDForce", "getReferencePositions") : ("unit.nanometer", ()),
 ("RMSDForce", "getParticles") : (None, ()),
+#("NoseHooverChain", "getThermostatedPairs") : (None, ()),
+#("NoseHooverChain", "getThermostatedAtoms") : (None, ()),
+#("NoseHooverChain", "getDefaultYoshidaSuzukiWeights") : (None, ()),
+#("NoseHooverIntegrator", "NoseHooverIntegrator") : (None, ("unit.kelvin", "1/unit.picosecond", "unit,picosecond", None, None, None)),
+("NoseHooverIntegrator", "setTemperature") : (None, ("unit.kelvin", None)),
+("NoseHooverIntegrator", "setRelativeTemperature") : (None, ("unit.kelvin", None) ),
+("NoseHooverIntegrator", "setCollisionFrequency") : (None, ("1/unit.picosecond", None)),
+("NoseHooverIntegrator", "setRelativeCollisionFrequency") : (None, ("1/unit.picosecond", None)),
+("NoseHooverIntegrator", "computeHeatBathEnergy") : ( "unit.kilojoules_per_mole", ()),
+("NoseHooverIntegrator", "addThermostat"): (None, ("unit.kelvin", "1/unit.picosecond", None, None, None)),
+("NoseHooverIntegrator", "addSubsystemThermostat"):
+    (None, (None, None, "unit.kelvin", "1/unit.picosecond", "unit.kelvin", "1/unit.picosecond", None, None, None))
 }
 

wrappers/python/tests/TestAPIUnits.py

@@ -1165,6 +1165,21 @@ class TestAPIUnits(unittest.TestCase):
         self.assertAlmostEqualUnit(integrator.getFriction(), 0.1/microsecond)
         self.assertEqual(integrator.getStepSize(), 1*femtosecond)
 
+    def testNoseHooverIntegrator(self):
+        """ Tests the NoseHooverIntegrator API features """
+        integrator = NoseHooverIntegrator(300, 0.1, 1)
+        self.assertEqual(integrator.getTemperature(0), 300*kelvin)
+        self.assertEqual(integrator.getCollisionFrequency(), 0.1/picosecond)
+        self.assertEqual(integrator.getStepSize(), 1*picosecond)
+
+        integrator = NoseHooverIntegrator(300*kelvin, 0.1/microsecond, 1*femtosecond)
+        self.assertEqual(integrator.getTemperature(), 300*kelvin)
+        self.assertAlmostEqualUnit(integrator.getCollisionFrequency(), 0.1/microsecond)
+        self.assertEqual(integrator.getStepSize(), 1*femtosecond)
+
+        integrator = NoseHooverIntegrator(1*femtosecond)
+        self.assertEqual(integrator.getStepSize(), 1*femtosecond)
+
     def testAndersenThermostat(self):
         """ Tests the AndersenThermostat API features """
         force = AndersenThermostat(300*kelvin, 1/microsecond)