Initial velocity verlet, with some nose hoover implementation

olla/include/openmm/kernels.h

@@ -55,6 +55,7 @@
 #include "openmm/KernelImpl.h"
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/MonteCarloBarostat.h"
+#include "openmm/NoseHooverChain.h"
 #include "openmm/PeriodicTorsionForce.h"
 #include "openmm/RBTorsionForce.h"
 #include "openmm/RMSDForce.h"
@@ -63,6 +64,7 @@
 #include "openmm/VariableLangevinIntegrator.h"
 #include "openmm/VariableVerletIntegrator.h"
 #include "openmm/VerletIntegrator.h"
+#include "openmm/VelocityVerletIntegrator.h"
 #include <iosfwd>
 #include <set>
 #include <string>
@@ -1057,6 +1059,39 @@ public:
     virtual double computeKineticEnergy(ContextImpl& context, const VerletIntegrator& integrator) = 0;
 };
 
+/**
+ * This kernel is invoked by VelocityVerletIntegrator to take one time step.
+ */
+class IntegrateVelocityVerletStepKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "IntegrateVelocityVerletStep";
+    }
+    IntegrateVelocityVerletStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param integrator the VelocityVerletIntegrator this kernel will be used for
+     */
+    virtual void initialize(const System& system, const VelocityVerletIntegrator& integrator) = 0;
+    /**
+     * Execute the kernel.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VelocityVerletIntegrator this kernel is being used for
+     */
+    virtual void execute(ContextImpl& context, const VelocityVerletIntegrator& integrator) = 0;
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VelocityVerletIntegrator this kernel is being used for
+     */
+    virtual double computeKineticEnergy(ContextImpl& context, const VelocityVerletIntegrator& integrator) = 0;
+};
+
 /**
  * This kernel is invoked by LangevinIntegrator to take one time step.
  */
@@ -1289,6 +1324,32 @@ public:
     virtual void execute(ContextImpl& context) = 0;
 };
 
+/**
+ * This kernel is invoked by NoseHooverChainThermostat at the beginning and end of each time step
+ * to update the Nose-Hoover chain.
+ */
+class PropagateNoseHooverChainKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "PropagateNoseHooverChain";
+    }
+    PropagateNoseHooverChainKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param noseHooverChain the NoseHooverChain this kernel will be used for
+     */
+    virtual void initialize(const System& system, const NoseHooverChain& noseHooverChain) = 0;
+    /**
+     * Execute the kernel.
+     * 
+     * @param context    the context in which to execute this kernel
+     */
+    virtual void execute(ContextImpl& context, const NoseHooverChain &noseHooverChain, double kineticEnergy, double timeStep) = 0;
+};
+
 /**
  * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
  */

openmmapi/include/OpenMM.h

@@ -61,6 +61,7 @@
 #include "openmm/MonteCarloBarostat.h"
 #include "openmm/MonteCarloMembraneBarostat.h"
 #include "openmm/NonbondedForce.h"
+#include "openmm/NoseHooverChain.h"
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/PeriodicTorsionForce.h"
@@ -74,6 +75,7 @@
 #include "openmm/VariableVerletIntegrator.h"
 #include "openmm/Vec3.h"
 #include "openmm/VerletIntegrator.h"
+#include "openmm/VelocityVerletIntegrator.h"
 #include "openmm/VirtualSite.h"
 #include "openmm/Platform.h"
 #include "openmm/serialization/XmlSerializer.h"

openmmapi/include/openmm/NoseHooverChain.h

+#ifndef OPENMM_NOSEHOOVERCHAIN_H_
+#define OPENMM_NOSEHOOVERCHAIN_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 "Force.h"
+#include <string>
+#include <vector>
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This class defines a chain of Nose-Hoover particles to be used as a heat bath to
+ * scale the velocities of a collection of particles subject to thermostating.  The
+ * heat bath is propagated using the multi time step approach detailed in
+ *
+ * G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Klein, Mol. Phys. 87, 1117 (1996).
+ *
+ * where the total number of timesteps used to propagate the chain in each step is
+ * the number of MTS steps multiplied by the number of terms in the Yoshida-Suzuki decomposition.
+ */
+
+class OPENMM_EXPORT NoseHooverChain : public Force {
+public:
+    /**
+     * This is the name of the parameter that stores the current temperature of the
+     * heat bath (in Kelvin).
+     */
+    std::string Temperature() const {
+        return defaultLabel + "NoseHooverChainTemperature";
+    }
+    /**
+     * This is the name of the parameter that stores the current collision frequency (in 1/ps).
+     */
+    std::string CollisionFrequency() const {
+        return defaultLabel + "NoseHooverChainCollisionFrequency";
+    }
+    /**
+     * This is the name of the parameter that stores the current number of degrees of freedom.
+     */
+    std::string NumDegreesOfFreedom() const {
+        return defaultLabel + "NoseHooverChainNumDegreesOfFreedom";
+    }
+    /**
+     * This is the name of the parameter that stores the current chain length
+     */
+    std::string ChainLength() const {
+        return defaultLabel + "NoseHooverChainLength";
+    }
+    /**
+     * This is the name of the parameter that stores the current number of multi time steps
+     */
+    std::string NumMultiTimeSteps() const {
+        return defaultLabel + "NoseHooverChainNumMultiTimeSteps";
+    }
+    /**
+     * This is the name of the parameter that stores the current number of Yoshida Suzuki time steps
+     */
+    std::string NumYoshidaSuzukiTimeSteps() const {
+        return defaultLabel + "NoseHooverChainNumYoshidaSuzukiTimeSteps";
+    }
+    /**
+     * This is the name of the parameter that stores the force acting on the ith bead
+     */
+    std::string Force(int i) const {
+        return defaultLabel + "NoseHooverChainForce" + std::to_string(i);
+    }
+    /**
+     * This is the name of the parameter that stores the mass of the ith bead
+     */
+    std::string Mass(int i) const {
+        return defaultLabel + "NoseHooverChainMass" + std::to_string(i);
+    }
+    /**
+     * This is the name of the parameter that stores the position of the ith bead
+     */
+    std::string Position(int i) const {
+        return defaultLabel + "NoseHooverChainPosition" + std::to_string(i);
+    }
+    /**
+     * This is the name of the parameter that stores the velocity of the ith bead
+     */
+    std::string Velocity(int i) const {
+        return defaultLabel + "NoseHooverChainVelocity" + std::to_string(i);
+    }
+    /**
+     * Create a NoseHooverChain.
+     * 
+     * @param defaultTemperature        the default temperature of the heat bath (in Kelvin)
+     * @param defaultCollisionFrequency the default collision frequency (in 1/ps)
+     * @param defaultNumDOFs            the default number of degrees of freedom in the particles that
+     *                                  interact with this chain
+     * @param defaultChainLength        the default length of (number of particles in) this heat bath
+     * @param defaultNumMTS             the default number of multi time steps used to propagate this chain
+     * @param defaultNumYoshidaSuzuki   the default number of Yoshida Suzuki steps used to propagate this chain (1, 3, or 5).
+     * @param defaultLabel              the default label used to distinguish this Nose-Hoover chain from others that may
+     *                                  be used to control a different set of particles, e.g. for Drude oscillators
+     */
+    NoseHooverChain(double defaultTemperature, double defaultCollisionFrequency, int defaultNumDOFs, int defaultChainLength,
+                    int defaultNumMTS, int defaultNumYoshidaSuzuki, const std::string &defaultLabel);
+    /**
+     * Get the default temperature of the heat bath (in Kelvin).
+     *
+     * @return the default temperature of the heat bath, measured in Kelvin.
+     */
+    double getDefaultTemperature() const {
+        return defaultTemp;
+    }
+    /**
+     * Set the default temperature of the heat bath.  This will affect any new Contexts you create,
+     * but not ones that already exist.
+     *
+     * @param temperature the default temperature of the heat bath (in Kelvin)
+     */
+    void setDefaultTemperature(double temperature) {
+        defaultTemp = temperature;
+    }
+    /**
+     * Get the default collision frequency (in 1/ps).
+     *
+     * @return the default collision frequency, measured in 1/ps.
+     */
+    double getDefaultCollisionFrequency() const {
+        return defaultFreq;
+    }
+    /**
+     * Set the default collision frequency.  This will affect any new Contexts you create,
+     * but not those that already exist.
+     *
+     * @param frequency the default collision frequency (in 1/ps)
+     */
+    void setDefaultCollisionFrequency(double frequency) {
+        defaultFreq = frequency;
+    }
+    /**
+     * Get the default number of degrees of freedom in the particles controled by this heat bath.
+     *
+     * @return the default number of degrees of freedom.
+     */
+    int getDefaultNumDegreesOfFreedom() const {
+        return defaultNumDOFs;
+    }
+    /**
+     * Set the default number of degrees of freedom in the particles controled by this heat bath.
+     * This will affect any new Contexts you create, but not those that already exist.
+     *
+     * @param numDOFs 
+     */
+    void setDefaultNumDegreesOfFreedom(int numDOFs) {
+        defaultNumDOFs = numDOFs;
+    }
+    /**
+     * Get the default chain length of this heat bath.
+     *
+     * @return the default chain length.
+     */
+    int getDefaultChainLength() const {
+        return defaultChainLength;
+    }
+    /**
+     * Set the default chain length of this heat bath.
+     * This will affect any new Contexts you create, but not those that already exist.
+     *
+     * @param numDOFs 
+     */
+    void setDefaultChainLength(int chainLength) {
+        defaultChainLength = chainLength;
+    }
+    /**
+     * Get the default number of steps used in the multi time step propagation.
+     *
+     * @returns the default number of multi time steps.
+     */
+    int getDefaultNumMultiTimeSteps() const {
+        return defaultNumMTS;
+    }
+    /**
+     * Set the default number of steps used in the multi time step propagation.
+     * This will affect any new Contexts you create, but not those that already exist.
+     *
+     * @param numSteps 
+     */
+    void setDefaultNumMultiTimeSteps(int numSteps) {
+        defaultNumMTS = numSteps;
+    }
+    /**
+     * Get the default number of steps used in the Yoshida-Suzuki decomposition for
+     * multi time step propagation.
+     *
+     * @returns the default number of multi time steps in the Yoshida-Suzuki decomposition.
+     */
+    int getDefaultNumYoshidaSuzukiTimeSteps() const {
+        return defaultNumYS;
+    }
+    /**
+     * Set the default number of steps used in the Yoshida-Suzuki decomposition for
+     * multi time step propagation. This will affect any new Contexts you create,
+     * but not those that already exist.
+     * 
+     * @param the default number of multi time steps in the Yoshida-Suzuki decomposition.
+     */
+    void setDefaultNumYoshidaSuzukiTimeSteps(int numSteps) {
+        defaultNumYS = numSteps;
+    }
+    /**
+     * Get the default label used to identify this chain
+     *
+     * @returns the default label
+     */
+    std::string getDefaultLabel() const {
+        return defaultLabel;
+    }
+    /**
+     * Set the default label used to identify this chain
+     * This will affect any new Contexts you create, but not those that already exist.
+     *
+     * @param the default label
+     */
+    void setDefaultLabel(const std::string& label) {
+        defaultLabel = label;
+    }
+    /**
+     * Get the default weights used in the Yoshida Suzuki multi time step decomposition (dimensionless) 
+     *
+     * @returns the weights for the Yoshida-Suzuki integration
+     */
+    std::vector<double> getDefaultYoshidaSuzukiWeights() const;
+    /**
+     * Returns whether or not this force makes use of periodic boundary
+     * conditions.
+     *
+     * @returns true if force uses PBC and false otherwise
+     */
+    bool usesPeriodicBoundaryConditions() const {
+        return false;
+    }
+protected:
+    ForceImpl* createImpl() const;
+private:
+    double defaultTemp, defaultFreq, defaultTimeStep;
+    int defaultNumDOFs, defaultChainLength, defaultNumMTS, defaultNumYS;
+    // The suffix used to distinguish NH chains, e.g. for Drude particles vs. regular particles.
+    std::string defaultLabel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_NOSEHOOVERCHAIN_H_*/

openmmapi/include/openmm/VelocityVerletIntegrator.h

+#ifndef OPENMM_VELOCITYVERLETINTEGRATOR_H_
+#define OPENMM_VELOCITYVERLETINTEGRATOR_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 "Integrator.h"
+#include "openmm/Kernel.h"
+#include "openmm/NoseHooverChain.h"
+#include "internal/windowsExport.h"
+
+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.
+ */
+
+class OPENMM_EXPORT VelocityVerletIntegrator : public Integrator {
+public:
+    /**
+     * Create a VelocityVerletIntegrator.
+     * 
+     * @param stepSize the step size with which to integrate the system (in picoseconds)
+     */
+    explicit VelocityVerletIntegrator(double stepSize);
+   /**
+     * 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);
+   /**
+     * Advance a simulation through time by taking a series of time steps.
+     * 
+     * @param kineticEnergy  the kinetic energy of the system that the chain is thermostating
+     * @param chainID        id of the Nose-Hoover-Chain
+     */
+    double propagateChain(double kineticEnergy, int chainID=0);
+    int addNoseHooverChainThermostat(System& system, double temperature, double collisionFrequency,
+                                     int chainLength, int numMTS, int numYoshidaSuzuki);
+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);
+    /**
+     * 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.
+     */
+    void cleanup();
+    /**
+     * Get the names of all Kernels used by this Integrator.
+     */
+    std::vector<std::string> getKernelNames();
+    /**
+     * Compute the kinetic energy of the system at the current time.
+     */
+    double computeKineticEnergy();
+private:
+    Kernel vvKernel;
+    std::vector<NoseHooverChain*> noseHooverChains;
+    std::vector<Kernel> nhcKernels;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_VELOCITYVERLETINTEGRATOR_H_*/

openmmapi/include/openmm/internal/NoseHooverChainImpl.h

+#ifndef OPENMM_NOSEHOOVERCHAINIMPL_H_
+#define OPENMM_NOSEHOOVERCHAINIMPL_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 "ForceImpl.h"
+#include "openmm/NoseHooverChain.h"
+#include "openmm/Kernel.h"
+#include <string>
+
+namespace OpenMM {
+
+class System;
+
+/**
+ * This is the internal implementation of NoseHooverChain.
+ */
+
+class OPENMM_EXPORT NoseHooverChainImpl : public ForceImpl {
+public:
+    NoseHooverChainImpl(const NoseHooverChain& owner);
+    void initialize(ContextImpl& context);
+    const NoseHooverChain& getOwner() const {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context, bool& forcesInvalid);
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
+    std::map<std::string, double> getDefaultParameters();
+    std::vector<std::string> getKernelNames();
+private:
+    const NoseHooverChain& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_NOSEHOOVERCHAINIMPL_H_*/

openmmapi/src/NoseHooverChain.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/NoseHooverChain.h"
+#include "openmm/internal/NoseHooverChainImpl.h"
+#include "openmm/OpenMMException.h"
+
+using namespace OpenMM;
+
+NoseHooverChain::NoseHooverChain(double defaultTemperature, double defaultCollisionFrequency, 
+                                 int defaultNumDOFs, int defaultChainLength, int defaultNumMTS,
+                                 int defaultNumYoshidaSuzuki, const std::string &defaultLabel) :
+        defaultTemp(defaultTemperature), defaultFreq(defaultCollisionFrequency), defaultNumDOFs(defaultNumDOFs),
+        defaultChainLength(defaultChainLength), defaultNumMTS(defaultNumMTS), defaultNumYS(defaultNumYS),
+        defaultLabel(defaultLabel) {}
+
+ForceImpl* NoseHooverChain::createImpl() const {
+    return new NoseHooverChainImpl(*this);
+}
+
+std::vector<double> NoseHooverChain::getDefaultYoshidaSuzukiWeights() const {
+    switch (defaultNumYS) {
+        case 1:
+            return {1};
+        case 3:
+            return {0.828981543588751, -0.657963087177502, 0.828981543588751};
+        case 5:
+            return {0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065,
+                    0.2967324292201065};
+        default:
+            throw OpenMMException("The number of Yoshida-Suzuki weights must be 1,3, or 5.");
+    }
+}

openmmapi/src/NoseHooverChainImpl.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) 2008-2010 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/internal/NoseHooverChainImpl.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/Integrator.h"
+#include "openmm/System.h"
+#include "openmm/kernels.h"
+#include "SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
+#include <string>
+#include <vector>
+
+using namespace OpenMM;
+using std::vector;
+
+NoseHooverChainImpl::NoseHooverChainImpl(const NoseHooverChain& owner) : owner(owner) {
+}
+
+void NoseHooverChainImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(PropagateNoseHooverChainKernel::Name(), context);
+    kernel.getAs<PropagateNoseHooverChainKernel>().initialize(context.getSystem(), owner);
+}
+
+void NoseHooverChainImpl::updateContextState(ContextImpl& context, bool& forcesInvalid) {
+    // This kernel updates the Nose-Hoover particles when invoked explicitly
+    // via its propagate(), which is called from the integrator.
+}
+
+double NoseHooverChainImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+    // This kernel doesn't compute an energy or force.  Instead it is invoked
+    // directly by the integrator via it's propagate() method.
+    return 0;
+}
+
+std::map<std::string, double> NoseHooverChainImpl::getDefaultParameters() {
+    std::map<std::string, double> parameters;
+    const auto &owner = getOwner();
+    double frequency = owner.getDefaultCollisionFrequency();
+    int chainLength = owner.getDefaultChainLength();
+    double T = owner.getDefaultTemperature();
+    double kT = BOLTZ * T;
+    int DOFs = owner.getDefaultNumDegreesOfFreedom();
+    parameters[owner.Temperature()] = T;
+    parameters[owner.CollisionFrequency()] = frequency;
+    parameters[owner.ChainLength()] = chainLength;
+    parameters[owner.NumYoshidaSuzukiTimeSteps()] = owner.getDefaultNumYoshidaSuzukiTimeSteps();
+    parameters[owner.NumMultiTimeSteps()] = owner.getDefaultNumMultiTimeSteps();
+    parameters[owner.NumDegreesOfFreedom()] = DOFs;
+    for(int i = 0; i < chainLength; ++i) {
+        parameters[owner.Force(i)] = 0;
+        parameters[owner.Position(i)] = 0;
+        parameters[owner.Mass(i)] = kT / (frequency * frequency);
+    }
+    parameters[owner.Mass(0)] *= DOFs;
+
+    // Set the velocities to the appropriate Boltzmann distribution;
+    // this is copied from Context::setVelocitiesToTemperature()
+    OpenMM_SFMT::SFMT sfmt;
+    int randomSeed = 0; //TODO figure out where / how this should be handled
+    init_gen_rand(randomSeed, sfmt);
+    vector<double> randoms;
+    while (randoms.size() < chainLength) {
+        double x, y, r2;
+        do {
+            x = 2.0*genrand_real2(sfmt)-1.0;
+            y = 2.0*genrand_real2(sfmt)-1.0;
+            r2 = x*x + y*y;
+        } while (r2 >= 1.0 || r2 == 0.0);
+        double multiplier = sqrt((-2.0*log(r2))/r2);
+        randoms.push_back(x*multiplier);
+        randoms.push_back(y*multiplier);
+    }
+    int nextRandom = 0;
+    for (int i = 0; i < chainLength; i++) {
+        double velocityScale = 1 / frequency; // = sqrt(kT / Q) = sqrt(kT / [ kT frequency^-2])
+        parameters[owner.Velocity(i)] = velocityScale * randoms[nextRandom++];
+    }
+    // N.B. the zeroth entry is computed as a function of the instantaneous KE at the start of propagate
+    for(int i = 1; i < chainLength-1; ++i) {
+        const double & v = parameters[owner.Velocity(i)];
+        parameters[owner.Force(i+1)] = (parameters[owner.Mass(i)] * v * v - kT) / parameters[owner.Mass(i+1)];
+    }
+
+
+    return parameters;
+}
+
+std::vector<std::string> NoseHooverChainImpl::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(PropagateNoseHooverChainKernel::Name());
+    return names;
+}

openmmapi/src/VelocityVerletIntegrator.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/VelocityVerletIntegrator.h"
+#include "openmm/Context.h"
+#include "openmm/Force.h"
+#include "openmm/System.h"
+#include "openmm/NoseHooverChain.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/CMMotionRemover.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/kernels.h"
+#include <string>
+
+using namespace OpenMM;
+using std::string;
+using std::vector;
+
+VelocityVerletIntegrator::VelocityVerletIntegrator(double stepSize) {
+    setStepSize(stepSize);
+    setConstraintTolerance(1e-5);
+}
+
+double VelocityVerletIntegrator::propagateChain(double kineticEnergy, int chainID) {
+}
+
+/*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) {
+    std::vector<int> mask, parents;
+    int nDOF = 0;
+    int numForces = system.getNumForces();
+    for(int particle = 0; particle < system.getNumParticles(); ++particle) {
+        if(system.getParticleMass(particle) > 0) nDOF += 3;
+    }
+    nDOF -= system.getNumConstraints();
+    for (int forceNum = 0; forceNum < numForces; ++forceNum) {
+        if (dynamic_cast<CMMotionRemover*>(&system.getForce(forceNum))) nDOF -= 3;
+    }
+
+    auto nhcForce = new NoseHooverChain(temperature, collisionFrequency, nDOF, chainLength,
+                                        numMTS, numYoshidaSuzuki, std::to_string(noseHooverChains.size()));
+    system.addForce(nhcForce);
+    noseHooverChains.push_back(nhcForce);
+}
+
+double VelocityVerletIntegrator::computeKineticEnergy() {
+    return vvKernel.getAs<IntegrateVelocityVerletStepKernel>().computeKineticEnergy(*context, *this);
+}
+
+void VelocityVerletIntegrator::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);
+    vvKernel.getAs<IntegrateVelocityVerletStepKernel>().initialize(contextRef.getSystem(), *this);
+    for (auto * nhc: noseHooverChains){
+        nhcKernels.push_back(context->getPlatform().createKernel(PropagateNoseHooverChainKernel::Name(), contextRef));
+        nhcKernels.back().getAs<PropagateNoseHooverChainKernel>().initialize(contextRef.getSystem(), *nhc);
+    }
+}
+
+void VelocityVerletIntegrator::cleanup() {
+    vvKernel = Kernel();
+    nhcKernels.clear();
+}
+
+vector<string> VelocityVerletIntegrator::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(PropagateNoseHooverChainKernel::Name());
+    names.push_back(IntegrateVelocityVerletStepKernel::Name());
+    return names;
+}
+
+void VelocityVerletIntegrator::step(int steps) {
+    if (context == NULL)
+        throw OpenMMException("This Integrator is not bound to a context!");
+    for (int i = 0; i < steps; ++i) {
+        context->updateContextState();
+        vvKernel.getAs<IntegrateVelocityVerletStepKernel>().execute(*context, *this);
+    }
+}

platforms/reference/include/ReferenceConstraints.h

@@ -61,7 +61,7 @@ public:
      * Apply the constraint algorithm to velocities.
      * 
      * @param atomCoordinates  the atom coordinates
-     * @param atomCoordinatesP the velocities to modify
+     * @param velocities       the velocities to modify
      * @param inverseMasses    1/mass
      * @param tolerance        the constraint tolerance
      */

platforms/reference/include/ReferenceKernels.h

@@ -58,7 +58,9 @@ class ReferenceGayBerneForce;
 class ReferenceBrownianDynamics;
 class ReferenceStochasticDynamics;
 class ReferenceConstraintAlgorithm;
+class ReferenceNoseHooverChain;
 class ReferenceMonteCarloBarostat;
+class ReferenceVelocityVerletDynamics;
 class ReferenceVariableStochasticDynamics;
 class ReferenceVariableVerletDynamics;
 class ReferenceVerletDynamics;
@@ -1131,6 +1133,43 @@ private:
     ReferenceVerletDynamics* dynamics;
     std::vector<double> masses;
     double prevStepSize;
+}
+;
+/**
+ * This kernel is invoked by VelocityVerletIntegrator to take one time step.
+ */
+class ReferenceIntegrateVelocityVerletStepKernel : public IntegrateVelocityVerletStepKernel {
+public:
+    ReferenceIntegrateVelocityVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVelocityVerletStepKernel(name, platform),
+        data(data), dynamics(0) {
+    }
+    ~ReferenceIntegrateVelocityVerletStepKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param integrator the VelocityVerletIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const VelocityVerletIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VerletIntegrator this kernel is being used for
+     */
+    void execute(ContextImpl& context, const VelocityVerletIntegrator& integrator);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param integrator the VelocityVerletIntegrator this kernel is being used for
+     */
+    double computeKineticEnergy(ContextImpl& context, const VelocityVerletIntegrator& integrator);
+private:
+    ReferencePlatform::PlatformData& data;
+    ReferenceVelocityVerletDynamics* dynamics;
+    std::vector<double> masses;
+    double prevStepSize;
 };
 
 /**
@@ -1387,6 +1426,32 @@ private:
     std::vector<double> masses;
 };
 
+/**
+ * This kernel is invoked by NoseHooverChain at the start of each time step to adjust the thermostat
+ * and update the associated particle velocities.
+ */
+class ReferencePropagateNoseHooverChainKernel : public PropagateNoseHooverChainKernel {
+public:
+    ReferencePropagateNoseHooverChainKernel(std::string name, const Platform& platform) : PropagateNoseHooverChainKernel(name, platform), chainPropagator(0) {
+    }
+    ~ReferencePropagateNoseHooverChainKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param thermostat the NoseHooverChain this kernel will be used for
+     */
+    void initialize(const System& system, const NoseHooverChain& chain);
+    /**
+     * Execute the kernel.
+     * 
+     * @param context    the context in which to execute this kernel
+     */
+    void execute(ContextImpl& context, const NoseHooverChain &nhc, double kineticEnergy, double timeStep);
+private:
+    ReferenceNoseHooverChain* chainPropagator;
+};
+
 /**
  * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
  */

platforms/reference/include/ReferenceNoseHooverChain.h

+
+/* Portions copyright (c) 2019 Stanford University and Simbios.
+ * Contributors: Andreas Krämer and Andrew C. Simmonett
+ *
+ * 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.
+ */
+
+#ifndef __ReferenceNoseHooverChain_H__
+#define __ReferenceNoseHooverChain_H__
+
+#include "openmm/Vec3.h"
+#include <vector>
+
+namespace OpenMM {
+
+using std::vector;
+
+class ReferenceNoseHooverChain {
+
+   private:
+       
+   public:
+
+      /**---------------------------------------------------------------------------------------
+      
+         Constructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceNoseHooverChain();
+
+      /**---------------------------------------------------------------------------------------
+      
+         Destructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceNoseHooverChain();
+
+      /**---------------------------------------------------------------------------------------
+
+         Propagate the Nose-Hoover chain a half timestep and find the appropriate velocity scaling
+
+         @param kineticEnergy      the instantaneous kinetic energy of the particles being thermostated
+         @param chainMasses        the "masses" assigned to each thermostat bead in ps^2 kJ / mol
+         @param chainVelocities    the velocities of the chain's beads in nm / ps
+         @param chainPositions     the positions of the chains's beads in nm
+         @param chainForces        the forces on each bead in the chain
+         @param numDOFs            the number of degrees of freedom in the system that this chain thermostats
+         @param temperature        thermostat temperature in Kelvin
+         @param collisionFrequency collision frequency for each atom in ps^-1
+         @param timeStep           full integration step size in ps (this only propagates half way)
+         @param numMTS             number of multi timestep increments used in the Trotter expansion
+         @param YSWeights          vector of weights used in the Yoshida-Suzuki multi-timestepping.
+         --------------------------------------------------------------------------------------- */
+      double propagate(double kineticEnergy, const vector<double>& chainMasses,
+                       vector<double>& chainVelocities, vector<double>& chainPositions,
+                       vector<double>& chainForces, int numDOFs,
+                       double temperature, double collisionFrequency, double timeStep,
+                       int numMTS, const vector<double>& YSWeights) const;
+};
+
+} // namespace OpenMM
+
+#endif // __ReferenceNoseHooverChain_H__

platforms/reference/include/ReferenceSETTLEAlgorithm.h

@@ -75,7 +75,7 @@ public:
      * Apply the constraint algorithm to velocities.
      * 
      * @param atomCoordinates  the atom coordinates
-     * @param atomCoordinatesP the velocities to modify
+     * @param velocities       the velocities to modify
      * @param inverseMasses    1/mass
      * @param tolerance        the constraint tolerance
      */

platforms/reference/include/ReferenceVelocityVerletDynamics.h

+
+/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * 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.
+ */
+
+#ifndef __ReferenceVelocityVerletDynamics_H__
+#define __ReferenceVelocityVerletDynamics_H__
+
+#include "ReferenceDynamics.h"
+
+namespace OpenMM {
+
+class ContextImpl;
+
+class ReferenceVelocityVerletDynamics : public ReferenceDynamics {
+
+   private:
+      std::vector<OpenMM::Vec3> xPrime;
+      std::vector<double> inverseMasses;
+      
+   public:
+
+      /**---------------------------------------------------------------------------------------
+      
+         Constructor
+
+         @param numberOfAtoms  number of atoms
+         @param deltaT         delta t for dynamics
+         @param friction       friction coefficient
+         @param temperature    temperature
+      
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceVelocityVerletDynamics(int numberOfAtoms, double deltaT);
+
+      /**---------------------------------------------------------------------------------------
+      
+         Destructor
+      
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceVelocityVerletDynamics();
+
+      /**---------------------------------------------------------------------------------------
+      
+         Update
+      
+         @param system              the System to be integrated
+         @param atomCoordinates     atom coordinates
+         @param velocities          velocities
+         @param forces              forces
+         @param masses              atom masses
+         @param tolerance           the constraint tolerance
+      
+         --------------------------------------------------------------------------------------- */
+     
+      void update(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
+                  std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance);
+      
+};
+
+} // namespace OpenMM
+
+#endif // __ReferenceVelocityVerletDynamics_H__

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -88,6 +88,10 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcGayBerneForceKernel(name, platform);
     if (name == IntegrateVerletStepKernel::Name())
         return new ReferenceIntegrateVerletStepKernel(name, platform, data);
+    if (name == IntegrateVelocityVerletStepKernel::Name())
+        return new ReferenceIntegrateVelocityVerletStepKernel(name, platform, data);
+    if (name == PropagateNoseHooverChainKernel::Name())
+        return new ReferencePropagateNoseHooverChainKernel(name, platform);
     if (name == IntegrateLangevinStepKernel::Name())
         return new ReferenceIntegrateLangevinStepKernel(name, platform, data);
     if (name == IntegrateBrownianStepKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -55,6 +55,7 @@
 #include "ReferenceLJCoulomb14.h"
 #include "ReferenceLJCoulombIxn.h"
 #include "ReferenceMonteCarloBarostat.h"
+#include "ReferenceNoseHooverChain.h"
 #include "ReferenceProperDihedralBond.h"
 #include "ReferenceRbDihedralBond.h"
 #include "ReferenceRMSDForce.h"
@@ -62,6 +63,7 @@
 #include "ReferenceTabulatedFunction.h"
 #include "ReferenceVariableStochasticDynamics.h"
 #include "ReferenceVariableVerletDynamics.h"
+#include "ReferenceVelocityVerletDynamics.h"
 #include "ReferenceVerletDynamics.h"
 #include "ReferenceVirtualSites.h"
 #include "openmm/CMMotionRemover.h"
@@ -2105,6 +2107,40 @@ double ReferenceIntegrateVerletStepKernel::computeKineticEnergy(ContextImpl& con
     return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
 }
 
+ReferenceIntegrateVelocityVerletStepKernel::~ReferenceIntegrateVelocityVerletStepKernel() {
+    if (dynamics)
+        delete dynamics;
+}
+
+void ReferenceIntegrateVelocityVerletStepKernel::initialize(const System& system, const VelocityVerletIntegrator& integrator) {
+    int numParticles = system.getNumParticles();
+    masses.resize(numParticles);
+    for (int i = 0; i < numParticles; ++i)
+        masses[i] = system.getParticleMass(i);
+}
+
+void ReferenceIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const VelocityVerletIntegrator& integrator) {
+    double stepSize = integrator.getStepSize();
+    vector<Vec3>& posData = extractPositions(context);
+    vector<Vec3>& velData = extractVelocities(context);
+    vector<Vec3>& forceData = extractForces(context);
+    if (dynamics == 0 || stepSize != prevStepSize) {
+        // Recreate the computation objects with the new parameters.
+        if (dynamics)
+            delete dynamics;
+        dynamics = new ReferenceVelocityVerletDynamics(context.getSystem().getNumParticles(), stepSize);
+        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
+        prevStepSize = stepSize;
+    }
+    dynamics->update(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
+    data.time += stepSize;
+    data.stepCount++;
+}
+
+double ReferenceIntegrateVelocityVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VelocityVerletIntegrator& integrator) {
+    return computeShiftedKineticEnergy(context, masses, 0);
+}
+
 ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
     if (dynamics)
         delete dynamics;
@@ -2386,6 +2422,39 @@ void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
         context.getIntegrator().getStepSize());
 }
 
+ReferencePropagateNoseHooverChainKernel::~ReferencePropagateNoseHooverChainKernel() {
+    if (chainPropagator)
+        delete chainPropagator;
+}
+
+void ReferencePropagateNoseHooverChainKernel::initialize(const System& system, const NoseHooverChain& nhc) {
+    this->chainPropagator = new ReferenceNoseHooverChain();
+    //SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
+}
+
+void ReferencePropagateNoseHooverChainKernel::execute(ContextImpl& context, const NoseHooverChain &nhc, double kineticEnergy, double timeStep) {
+    int chainLength = context.getParameter(nhc.ChainLength());
+    std::vector<double> chainPositions(chainLength), chainVelocities(chainLength), chainForces(chainLength),
+                        chainMasses(chainLength);
+    double temperature = context.getParameter(nhc.Temperature());
+    double collisionFrequency = context.getParameter(nhc.CollisionFrequency());
+    int numDOFs = context.getParameter(nhc.NumDegreesOfFreedom());
+    int numMTS = context.getParameter(nhc.NumMultiTimeSteps());
+    for(int i = 0; i < chainLength; ++i) {
+        chainPositions[i] = context.getParameter(nhc.Position(i));
+        chainVelocities[i] = context.getParameter(nhc.Velocity(i));
+        chainMasses[i] = context.getParameter(nhc.Mass(i));
+        chainForces[i] = context.getParameter(nhc.Force(i));
+    }
+
+    chainPropagator->propagate(kineticEnergy, chainMasses, chainVelocities, chainPositions, chainForces, numDOFs,
+                               temperature, collisionFrequency, timeStep, numMTS, nhc.getDefaultYoshidaSuzukiWeights());
+    //chain->applyThermostat(particleGroups, velData, masses,
+    //    context.getParameter(AndersenThermostat::Temperature()),
+    //    context.getParameter(AndersenThermostat::CollisionFrequency()),
+    //    context.getIntegrator().getStepSize());
+}
+
 ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel() {
     if (barostat)
         delete barostat;

platforms/reference/src/ReferencePlatform.cpp

@@ -66,6 +66,8 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
+    registerKernelFactory(IntegrateVelocityVerletStepKernel::Name(), factory);
+    registerKernelFactory(PropagateNoseHooverChainKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
     registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceNoseHooverChain.cpp

+
+/* Portions copyright (c) 2008-2010 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * 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 <cmath>
+#include <string.h>
+#include <sstream>
+
+#include "SimTKOpenMMUtilities.h"
+#include "ReferenceNoseHooverChain.h"
+
+using std::vector;
+using namespace OpenMM;
+
+/**---------------------------------------------------------------------------------------
+
+   Constructor
+
+   --------------------------------------------------------------------------------------- */
+
+ ReferenceNoseHooverChain::ReferenceNoseHooverChain() {
+ }
+
+/**---------------------------------------------------------------------------------------
+
+   Destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ ReferenceNoseHooverChain::~ReferenceNoseHooverChain() {
+ }
+
+double ReferenceNoseHooverChain::propagate(double kineticEnergy, const vector<double>& chainMasses,
+                                           vector<double>& chainVelocities, vector<double>& chainPositions,
+                                           vector<double>& chainForces, int numDOFs,
+                                           double temperature, double collisionFrequency, double timeStep,
+                                           int numMTS, const vector<double>& YSWeights) const {
+     double scale = 1;
+     double KE2 = 2 * kineticEnergy;
+     const double kT = BOLTZ * temperature;
+     const size_t chainLength = chainMasses.size();
+     chainForces[0] = (KE2 - numDOFs * kT) / chainMasses[0];
+     for (int mts = 0; mts < numMTS; ++mts) {
+         for (const auto &ys : YSWeights) {
+             double wdt = ys * timeStep / numMTS;
+             chainVelocities.back() += 0.25 * wdt + chainForces.back();
+             for (int bead = chainLength - 2; bead >= 0; --bead) {
+                 double aa = exp(-0.125 * wdt * chainVelocities[bead + 1]);
+                 chainVelocities[bead] = aa * (chainVelocities[bead] * aa + 0.25 * wdt * chainForces[bead]);
+             }
+             // update particle velocities
+             double aa = exp(-0.5 * wdt * chainVelocities[0]);
+             scale *= aa;
+             // update the thermostat positions
+             for (int bead = 0; bead < chainLength; ++bead) {
+                 chainPositions[bead] += 0.5 * chainVelocities[bead] * wdt;
+             }
+             // update the forces
+             chainForces[0] = (scale * scale * KE2 - numDOFs * kT) / chainMasses[0];
+             // update thermostat velocities
+             for (int bead = 0; bead < chainLength - 1; ++bead) {
+                 double aa = exp(-0.125 * wdt * chainVelocities[bead + 1]);
+                 chainVelocities[bead] = aa * (aa * chainVelocities[bead] + 0.25 * wdt * chainForces[bead]);
+                 chainForces[bead + 1] = (chainMasses[bead] * chainVelocities[bead] * chainVelocities[bead] - kT) / chainMasses[bead + 1];
+             }
+             chainVelocities.back() += 0.25 * wdt * chainForces.back();
+         }  // YS loop
+     } // MTS loop
+     return scale;
+}

platforms/reference/src/SimTKReference/ReferenceVelocityVerletDynamics.cpp

+
+/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
+ * Contributors: Peter Eastman, Pande Group
+ *
+ * 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 <cstring>
+#include <sstream>
+
+#include "SimTKOpenMMUtilities.h"
+#include "openmm/internal/ContextImpl.h"
+#include "ReferenceVelocityVerletDynamics.h"
+#include "ReferenceVirtualSites.h"
+
+#include <cstdio>
+
+using std::vector;
+using namespace OpenMM;
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceVelocityVerletDynamics constructor
+
+   @param numberOfAtoms  number of atoms
+   @param deltaT         delta t for dynamics
+   @param friction       friction coefficient
+   @param temperature    temperature
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceVelocityVerletDynamics::ReferenceVelocityVerletDynamics(int numberOfAtoms, double deltaT) :
+           ReferenceDynamics(numberOfAtoms, deltaT, 0.0) {
+   xPrime.resize(numberOfAtoms);
+   inverseMasses.resize(numberOfAtoms);
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceVelocityVerletDynamics destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceVelocityVerletDynamics::~ReferenceVelocityVerletDynamics() {
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Update -- driver routine for performing Velocity Verlet dynamics update of coordinates
+   and velocities
+
+   @param system              the System to be integrated
+   @param atomCoordinates     atom coordinates
+   @param velocities          velocities
+   @param forces              forces
+   @param masses              atom masses
+
+   --------------------------------------------------------------------------------------- */
+
+void ReferenceVelocityVerletDynamics::update(OpenMM::ContextImpl &context, const OpenMM::System& system, vector<Vec3>& atomCoordinates,
+                                          vector<Vec3>& velocities,
+                                          vector<Vec3>& forces, vector<double>& masses, double tolerance) {
+   // first-time-through initialization
+
+   int numberOfAtoms = system.getNumParticles();
+   if (getTimeStep() == 0) {
+      // invert masses
+
+      for (int ii = 0; ii < numberOfAtoms; ii++) {
+         if (masses[ii] == 0.0)
+             inverseMasses[ii] = 0.0;
+         else
+             inverseMasses[ii] = 1.0/masses[ii];
+      }
+   }
+
+   // Perform the integration.
+
+   for (int i = 0; i < numberOfAtoms; ++i) {
+       if (masses[i] != 0.0)
+           for (int j = 0; j < 3; ++j) {
+               velocities[i][j] += 0.5 * inverseMasses[i]*forces[i][j]*getDeltaT();
+               xPrime[i][j] = atomCoordinates[i][j];
+               atomCoordinates[i][j] += velocities[i][j]*getDeltaT();
+           }
+   }
+
+   // 
+
+   ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
+   if (referenceConstraintAlgorithm)
+      referenceConstraintAlgorithm->apply(xPrime, atomCoordinates, inverseMasses, tolerance);
+
+   context.calcForcesAndEnergy(true, false);
+
+   for (int i = 0; i < numberOfAtoms; ++i) {
+       if (masses[i] != 0.0)
+           for (int j = 0; j < 3; ++j) {
+               xPrime[i][j] += velocities[i][j]*getDeltaT();
+           }
+   }
+
+   // Update the positions and velocities.
+
+   for (int i = 0; i < numberOfAtoms; ++i) {
+       if (masses[i] != 0.0)
+           for (int j = 0; j < 3; ++j) {
+               velocities[i][j] += 0.5*inverseMasses[i]*forces[i][j]*getDeltaT() + (atomCoordinates[i][j] - xPrime[i][j]) / getDeltaT();
+           }
+   }
+   if (referenceConstraintAlgorithm)
+      referenceConstraintAlgorithm->applyToVelocities(atomCoordinates, velocities, inverseMasses, tolerance);
+
+   ReferenceVirtualSites::computePositions(system, atomCoordinates);
+   incrementTimeStep();
+}

platforms/reference/tests/TestReferenceNoseHooverThermostat.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 "ReferenceTests.h"
+#include "TestNoseHooverThermostat.h"
+
+void runPlatformTests() {
+}

platforms/reference/tests/TestReferenceVelocityVerletIntegrator.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) 2015 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 "TestVerletIntegrator.h"
+
+void runPlatformTests() {
+}