First steps towards DrudeIntegrator base class

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

+#ifndef OPENMM_DRUDEINTEGRATOR_H_
+#define OPENMM_DRUDEINTEGRATOR_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) 2008-2013 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/Integrator.h"
+#include "openmm/Kernel.h"
+#include "openmm/internal/windowsExportDrude.h"
+
+namespace OpenMM {
+
+/**
+ * This is a leap-frog Verlet Integrator that simulates systems with Drude particles.  It uses the
+ * self-consistent field (SCF) method: at every time step, the positions of Drude particles are
+ * adjusted to minimize the potential energy.
+ * 
+ * This Integrator requires the System to include a DrudeForce, which it uses to identify the Drude
+ * particles.
+ */
+
+class OPENMM_EXPORT_DRUDE DrudeIntegrator : public Integrator {
+public:
+    /**
+     * Create a DrudeSCFIntegrator.
+     *
+     * @param stepSize       the step size with which to integrator the system (in picoseconds)
+     */
+    DrudeIntegrator(double stepSize) {};
+    /**
+     * Advance a simulation through time by taking a series of time steps.
+     *
+     * @param steps   the number of time steps to take
+     */
+    virtual void step(int steps) override {};
+    /**
+     * Get the temperature of the heat bath applied to internal coordinates of Drude particles (in Kelvin).
+     *
+     * @return the temperature of the heat bath, measured in Kelvin
+     */
+    double getDrudeTemperature() const {
+        return drudeTemperature;
+    }
+    /**
+     * Set the temperature of the heat bath applied to internal coordinates of Drude particles (in Kelvin).
+     *
+     * @param temp    the temperature of the heat bath, measured in Kelvin
+     */
+    void setDrudeTemperature(double temp) {
+        drudeTemperature = temp;
+    }
+    /**
+     * Get the maximum distance a Drude particle can ever move from its parent particle, measured in nm.  This is implemented
+     * with a hard wall constraint.  If this distance is set to 0 (the default), the hard wall constraint is omitted.
+     */
+    double getMaxDrudeDistance() const;
+    /**
+     * Set the maximum distance a Drude particle can ever move from its parent particle, measured in nm.  This is implemented
+     * with a hard wall constraint.  If this distance is set to 0 (the default), the hard wall constraint is omitted.
+     */
+    void setMaxDrudeDistance(double distance);
+    /**
+     * Set the random number seed.  The precise meaning of this parameter is undefined, and is left up
+     * to each Platform to interpret in an appropriate way.  It is guaranteed that if two simulations
+     * are run with different random number seeds, the sequence of random forces will be different.  On
+     * the other hand, no guarantees are made about the behavior of simulations that use the same seed.
+     * In particular, Platforms are permitted to use non-deterministic algorithms which produce different
+     * results on successive runs, even if those runs were initialized identically.
+     *
+     * If seed is set to 0 (which is the default value assigned), a unique seed is chosen when a Context
+     * is created from this Force. This is done to ensure that each Context receives unique random seeds
+     * without you needing to set them explicitly.
+     */
+    void setRandomNumberSeed(int seed) {
+        randomNumberSeed = seed;
+    }
+    /**
+     * Get the random number seed.  See setRandomNumberSeed() for details.
+     */
+    int getRandomNumberSeed() const {
+        return randomNumberSeed;
+    }
+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.
+     */
+    virtual void initialize(ContextImpl& context) override {};
+    /**
+     * 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.
+     */
+    virtual void cleanup() override {};
+    /**
+     * When the user modifies the state, we need to mark that the forces need to be recalculated.
+     */
+    virtual void stateChanged(State::DataType changed) override {};
+    /**
+     * Get the names of all Kernels used by this Integrator.
+     */
+    virtual std::vector<std::string> getKernelNames() override { return std::vector<std::string>(); }
+    /**
+     * Compute the kinetic energy of the system at the current time.
+     */
+    virtual double computeKineticEnergy() override { return 0; }
+    /**
+     * Return a list of velocities normally distributed around a target temperature, with the Drude
+     * temperatures assigned according to the Drude temperature assigned to the integrator.
+     *
+     * @param system the system whose velocities are to be initialized.
+     * @param temperature the target temperature in Kelvin.
+     * @param randomSeed the random number seed to use when selecting velocities 
+     */
+    virtual std::vector<Vec3> getVelocitiesForTemperature(const System &system, double temperature,
+                                                          int randomSeed) const override;
+
+    int randomNumberSeed;
+    double drudeTemperature, maxDrudeDistance;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_DRUDEINTEGRATOR_H_*/

plugins/drude/openmmapi/src/DrudeIntegrator.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-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 "sfmt/SFMT.h"
+#include "SimTKOpenMMRealType.h"
+#include "openmm/DrudeForce.h"
+#include "openmm/DrudeIntegrator.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/System.h"
+
+#include <set>
+
+using namespace OpenMM;
+
+std::vector<Vec3> DrudeIntegrator::getVelocitiesForTemperature(const System &system, double temperature, int randomSeedIn) const {
+    // Find the underlying Drude force object
+    const DrudeForce* drudeForce = NULL;
+    for (int i = 0; i < system.getNumForces(); i++)
+        if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
+            if (drudeForce == NULL)
+                drudeForce = dynamic_cast<const DrudeForce*>(&system.getForce(i));
+            else
+                throw OpenMMException("The System contains multiple DrudeForces");
+        }
+    if (drudeForce == NULL)
+        throw OpenMMException("The System does not contain a DrudeForce");
+
+    // Figure out which particles are individual and which are Drude pairs
+    std::set<int> particles;
+    std::vector<std::pair<int, int>> pairParticles;
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        particles.insert(i);
+    }
+    for (int i = 0; i < drudeForce->getNumParticles(); i++) {
+        int p, p1, p2, p3, p4;
+        double charge, polarizability, aniso12, aniso34;
+        drudeForce->getParticleParameters(i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
+        particles.erase(p);
+        particles.erase(p1);
+        pairParticles.emplace_back(p, p1);
+    }
+    std::vector<int> normalParticles(particles.begin(), particles.end());
+
+    // Generate the list of Gaussian random numbers.
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(randomSeedIn, sfmt);
+    std::vector<double> randoms;
+    while (randoms.size() < system.getNumParticles()*3) {
+        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*std::log(r2))/r2);
+        randoms.push_back(x*multiplier);
+        randoms.push_back(y*multiplier);
+    }
+
+    // Assign the velocities.
+    std::vector<Vec3> velocities(system.getNumParticles(), Vec3());
+    int nextRandom = 0;
+    // First the indivitual atoms
+    for (const auto &atom : normalParticles ) {
+        double mass = system.getParticleMass(atom);
+        if (mass != 0) {
+            double velocityScale = sqrt(BOLTZ*temperature/mass);
+            velocities[atom] = Vec3(randoms[nextRandom++], randoms[nextRandom++], randoms[nextRandom++])*velocityScale;
+        }
+    }
+    // Now the particle-Drude pairs
+    for (const auto &pair : pairParticles ) {
+        const auto atom1 = pair.first;
+        const auto atom2 = pair.second;
+        double mass1 = system.getParticleMass(atom1);
+        double mass2 = system.getParticleMass(atom2);
+        if (mass1 != 0 && mass2 != 0) {
+            double invMass = 1.0 / (mass1 + mass2);
+            double redMass = mass1 * mass2 * invMass;
+            double fracM1 = mass1 * invMass;
+            double fracM2 = mass2 * invMass;
+            Vec3 comVelocity = Vec3(randoms[nextRandom++], randoms[nextRandom++], randoms[nextRandom++])*sqrt(BOLTZ*temperature*invMass);
+            Vec3 relVelocity = Vec3(randoms[nextRandom++], randoms[nextRandom++], randoms[nextRandom++])*sqrt(BOLTZ*drudeTemperature/redMass);
+            velocities[atom1] = comVelocity - fracM2 * relVelocity;
+            velocities[atom2] = comVelocity + fracM1 * relVelocity;
+        }
+    }
+    return velocities;
+}
+
+double DrudeIntegrator::getMaxDrudeDistance() const {
+    return maxDrudeDistance;
+}
+
+void DrudeIntegrator::setMaxDrudeDistance(double distance) {
+    if (distance < 0)
+        throw OpenMMException("setMaxDrudeDistance: Distance cannot be negative");
+    maxDrudeDistance = distance;
+}