Created DrudeSCFIntegrator, including reference implementation

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

@@ -34,6 +34,7 @@
 
 #include "openmm/DrudeForce.h"
 #include "openmm/DrudeLangevinIntegrator.h"
+#include "openmm/DrudeSCFIntegrator.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
 #include "openmm/Vec3.h"
@@ -108,6 +109,37 @@ public:
     virtual double computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) = 0;
 };
 
+/**
+ * This kernel is invoked by DrudeSCFIntegrator to take one time step.
+ */
+class IntegrateDrudeSCFStepKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "IntegrateDrudeSCFStep";
+    }
+    IntegrateDrudeSCFStepKernel(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 DrudeSCFIntegrator this kernel will be used for
+     * @param force      the DrudeForce to get particle parameters from
+     */
+    virtual void initialize(const System& system, const DrudeSCFIntegrator& integrator, const DrudeForce& force) = 0;
+    /**
+     * Execute the kernel.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param integrator     the DrudeSCFIntegrator this kernel is being used for
+     */
+    virtual void execute(ContextImpl& context, const DrudeSCFIntegrator& integrator) = 0;
+    /**
+     * Compute the kinetic energy.
+     */
+    virtual double computeKineticEnergy(ContextImpl& context, const DrudeSCFIntegrator& integrator) = 0;
+};
+
 } // namespace OpenMM
 
 #endif /*DRUDE_KERNELS_H_*/

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

+#ifndef OPENMM_DRUDESCFINTEGRATOR_H_
+#define OPENMM_DRUDESCFINTEGRATOR_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 DrudeSCFIntegrator : public Integrator {
+public:
+    /**
+     * Create a DrudeSCFIntegrator.
+     *
+     * @param stepSize       the step size with which to integrator the system (in picoseconds)
+     */
+    DrudeSCFIntegrator(double stepSize);
+    /**
+     * Get the error tolerance to use when minimizing the potential energy.  This roughly corresponds
+     * to the maximum allowed force magnitude on the Drude particles after minimization.
+     *
+     * @return the error tolerance to use, measured in kJ/mol/nm
+     */
+    double getMinimizationErrorTolerance() const {
+        return tolerance;
+    }
+    /**
+     * Set the error tolerance to use when minimizing the potential energy.  This roughly corresponds
+     * to the maximum allowed force magnitude on the Drude particles after minimization.
+     *
+     * @param tol    the error tolerance to use, measured in kJ/mol/nm
+     */
+    void setMinimizationErrorTolerance(double tol) {
+        tolerance = tol;
+    }
+    /**
+     * Advance a simulation through time by taking a series of time steps.
+     *
+     * @param steps   the number of time steps to take
+     */
+    void step(int steps);
+protected:
+    /**
+     * This will be called by the Context when it is created.  It informs the Integrator
+     * of what context it will be integrating, and gives it a chance to do any necessary initialization.
+     * It will also get called again if the application calls reinitialize() on the Context.
+     */
+    void initialize(ContextImpl& context);
+    /**
+     * 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();
+    /**
+     * When the user modifies the state, we need to mark that the forces need to be recalculated.
+     */
+    void stateChanged(State::DataType changed);
+    /**
+     * 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:
+    double tolerance;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_DRUDESCFINTEGRATOR_H_*/

plugins/drude/openmmapi/src/DrudeSCFIntegrator.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-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/DrudeSCFIntegrator.h"
+#include "openmm/Context.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/DrudeKernels.h"
+#include <cmath>
+#include <ctime>
+#include <string>
+
+using namespace OpenMM;
+using std::string;
+using std::vector;
+
+DrudeSCFIntegrator::DrudeSCFIntegrator(double stepSize) {
+    setStepSize(stepSize);
+    setMinimizationErrorTolerance(0.1);
+    setConstraintTolerance(1e-5);
+}
+
+void DrudeSCFIntegrator::initialize(ContextImpl& contextRef) {
+    if (owner != NULL && &contextRef.getOwner() != owner)
+        throw OpenMMException("This Integrator is already bound to a context");
+    const DrudeForce* force = NULL;
+    const System& system = contextRef.getSystem();
+    for (int i = 0; i < system.getNumForces(); i++)
+        if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
+            if (force == NULL)
+                force = dynamic_cast<const DrudeForce*>(&system.getForce(i));
+            else
+                throw OpenMMException("The System contains multiple DrudeForces");
+        }
+    if (force == NULL)
+        throw OpenMMException("The System does not contain a DrudeForce");
+    context = &contextRef;
+    owner = &contextRef.getOwner();
+    kernel = context->getPlatform().createKernel(IntegrateDrudeSCFStepKernel::Name(), contextRef);
+    kernel.getAs<IntegrateDrudeSCFStepKernel>().initialize(contextRef.getSystem(), *this, *force);
+}
+
+void DrudeSCFIntegrator::cleanup() {
+    kernel = Kernel();
+}
+
+void DrudeSCFIntegrator::stateChanged(State::DataType changed) {
+}
+
+vector<string> DrudeSCFIntegrator::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(IntegrateDrudeSCFStepKernel::Name());
+    return names;
+}
+
+double DrudeSCFIntegrator::computeKineticEnergy() {
+    return kernel.getAs<IntegrateDrudeSCFStepKernel>().computeKineticEnergy(*context, *this);
+}
+
+void DrudeSCFIntegrator::step(int steps) {
+    for (int i = 0; i < steps; ++i) {
+        context->updateContextState();
+        context->calcForcesAndEnergy(true, false);
+        kernel.getAs<IntegrateDrudeSCFStepKernel>().execute(*context, *this);
+    }
+}

plugins/drude/platforms/reference/src/ReferenceDrudeKernelFactory.cpp

@@ -49,6 +49,7 @@ extern "C" void initDrudeReferenceKernels() {
     ReferenceDrudeKernelFactory* factory = new ReferenceDrudeKernelFactory();
     platform.registerKernelFactory(CalcDrudeForceKernel::Name(), factory);
     platform.registerKernelFactory(IntegrateDrudeLangevinStepKernel::Name(), factory);
+    platform.registerKernelFactory(IntegrateDrudeSCFStepKernel::Name(), factory);
 }
 
 KernelImpl* ReferenceDrudeKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
@@ -57,5 +58,7 @@ KernelImpl* ReferenceDrudeKernelFactory::createKernelImpl(std::string name, cons
         return new ReferenceCalcDrudeForceKernel(name, platform);
     if (name == IntegrateDrudeLangevinStepKernel::Name())
         return new ReferenceIntegrateDrudeLangevinStepKernel(name, platform, data);
+    if (name == IntegrateDrudeSCFStepKernel::Name())
+        return new ReferenceIntegrateDrudeSCFStepKernel(name, platform, data);
     throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
 }

plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp

@@ -70,6 +70,40 @@ static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorith
     }
 }
 
+static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& inverseMasses, double timeShift, ReferenceConstraintAlgorithm* constraints) {
+    const System& system = context.getSystem();
+    int numParticles = system.getNumParticles();
+    vector<RealVec>& posData = extractPositions(context);
+    vector<RealVec>& velData = extractVelocities(context);
+    vector<RealVec>& forceData = extractForces(context);
+    
+    // Compute the shifted velocities.
+    
+    vector<RealVec> shiftedVel(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        if (inverseMasses[i] > 0)
+            shiftedVel[i] = velData[i]+forceData[i]*(timeShift*inverseMasses[i]);
+        else
+            shiftedVel[i] = velData[i];
+    }
+    
+    // Apply constraints to them.
+    
+    if (constraints != NULL) {
+        constraints->setTolerance(1e-4);
+        constraints->applyToVelocities(numParticles, posData, shiftedVel, inverseMasses);
+    }
+    
+    // Compute the kinetic energy.
+    
+    double energy = 0.0;
+    for (int i = 0; i < numParticles; ++i)
+        if (inverseMasses[i] > 0)
+            energy += (shiftedVel[i].dot(shiftedVel[i]))/inverseMasses[i];
+    return 0.5*energy;
+}
+
+
 void ReferenceCalcDrudeForceKernel::initialize(const System& system, const DrudeForce& force) {
     // Initialize particle parameters.
     
@@ -330,35 +364,156 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
 }
 
 double ReferenceIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
-    const System& system = context.getSystem();
-    int numParticles = system.getNumParticles();
-    vector<RealVec>& posData = extractPositions(context);
-    vector<RealVec>& velData = extractVelocities(context);
-    vector<RealVec>& forceData = extractForces(context);
+    return computeShiftedKineticEnergy(context, particleInvMass, 0.5*integrator.getStepSize(), constraints);
+}
 
-    // Compute the shifted velocities.
+ReferenceIntegrateDrudeSCFStepKernel::~ReferenceIntegrateDrudeSCFStepKernel() {
+    if (constraints != NULL)
+        delete constraints;
+    if (minimizerPos != NULL)
+        lbfgs_free(minimizerPos);
+}
 
-    vector<RealVec> shiftedVel(numParticles);
-    double timeShift = 0.5*integrator.getStepSize();
-    for (int i = 0; i < numParticles; ++i) {
-        if (particleInvMass[i] > 0)
-            shiftedVel[i] = velData[i]+forceData[i]*(timeShift*particleInvMass[i]);
-        else
-            shiftedVel[i] = velData[i];
+void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, const DrudeSCFIntegrator& integrator, const DrudeForce& force) {
+    // Identify Drude particles.
+    
+    for (int i = 0; i < force.getNumParticles(); i++) {
+        int p, p1, p2, p3, p4;
+        double charge, polarizability, aniso12, aniso34;
+        force.getParticleParameters(i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
+        drudeParticles.push_back(p);
     }
 
-    // Apply constraints to them.
+    // Record particle masses.
 
-    if (constraints != NULL) {
-        constraints->setTolerance(1e-4);
-        constraints->applyToVelocities(numParticles, posData, shiftedVel, particleInvMass);
+    vector<RealOpenMM> particleMass;
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        double mass = system.getParticleMass(i);
+        particleMass.push_back(mass);
+        particleInvMass.push_back(mass == 0.0 ? 0.0 : 1.0/mass);
     }
     
-    // Compute the kinetic energy.
+    // Prepare constraints.
     
-    double energy = 0.0;
-    for (int i = 0; i < numParticles; ++i)
-        if (particleInvMass[i] > 0)
-            energy += (shiftedVel[i].dot(shiftedVel[i]))/particleInvMass[i];
-    return 0.5*energy;
+    int numConstraints = system.getNumConstraints();
+    if (numConstraints > 0) {
+        vector<pair<int, int> > constraintIndices(numConstraints);
+        vector<RealOpenMM> constraintDistances(numConstraints);
+        for (int i = 0; i < numConstraints; ++i) {
+            int particle1, particle2;
+            double distance;
+            system.getConstraintParameters(i, particle1, particle2, distance);
+            constraintIndices[i].first = particle1;
+            constraintIndices[i].second = particle2;
+            constraintDistances[i] = static_cast<RealOpenMM>(distance);
+        }
+        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
+        findAnglesForCCMA(system, angles);
+        constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, particleMass, angles, (RealOpenMM)integrator.getConstraintTolerance());
+    }
+    
+    // Initialize the energy minimizer.
+    
+    minimizerPos = lbfgs_malloc(drudeParticles.size()*3);
+    if (minimizerPos == NULL)
+        throw OpenMMException("DrudeSCFIntegrator: Failed to allocate memory");
+    lbfgs_parameter_init(&minimizerParams);
+    minimizerParams.max_iterations = 0;
+    minimizerParams.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE;
+}
+
+void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
+    vector<RealVec>& pos = extractPositions(context);
+    vector<RealVec>& vel = extractVelocities(context);
+    vector<RealVec>& force = extractForces(context);
+    
+    // Update the positions and velocities.
+    
+    int numParticles = particleInvMass.size();
+    vector<RealVec> xPrime(numParticles);
+    RealOpenMM dt = integrator.getStepSize();
+    for (int i = 0; i < numParticles; i++) {
+        if (particleInvMass[i] != 0.0) {
+            vel[i] += force[i]*particleInvMass[i]*dt;
+            xPrime[i] = pos[i]+vel[i]*dt;
+        }
+    }
+        
+    // Apply constraints.
+    
+    if (constraints != NULL)
+        constraints->apply(numParticles, pos, xPrime, particleInvMass);
+    
+    // Record the constrained positions and velocities.
+    
+    RealOpenMM dtInv = 1.0/dt;
+    for (int i = 0; i < numParticles; i++) {
+        if (particleInvMass[i] != 0.0) {
+            vel[i] = (xPrime[i]-pos[i])*dtInv;
+            pos[i] = xPrime[i];
+        }
+    }
+    
+    // Update the positions of virtual sites and Drude particles.
+    
+    ReferenceVirtualSites::computePositions(context.getSystem(), pos);
+    minimize(context, integrator.getMinimizationErrorTolerance());
+    data.time += integrator.getStepSize();
+    data.stepCount++;
+}
+
+double ReferenceIntegrateDrudeSCFStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
+    return computeShiftedKineticEnergy(context, particleInvMass, 0.5*integrator.getStepSize(), constraints);
+}
+
+struct MinimizerData {
+    ContextImpl& context;
+    vector<int>& drudeParticles;
+    MinimizerData(ContextImpl& context, vector<int>& drudeParticles) : context(context), drudeParticles(drudeParticles) {}
+};
+
+static lbfgsfloatval_t evaluate(void *instance, const lbfgsfloatval_t *x, lbfgsfloatval_t *g, const int n, const lbfgsfloatval_t step) {
+    MinimizerData* data = reinterpret_cast<MinimizerData*>(instance);
+    ContextImpl& context = data->context;
+    vector<int>& drudeParticles = data->drudeParticles;
+    int numDrudeParticles = drudeParticles.size();
+
+    // Compute the force and energy for this configuration.
+
+    vector<RealVec>& pos = extractPositions(context);
+    vector<RealVec>& force = extractForces(context);
+    for (int i = 0; i < numDrudeParticles; i++)
+        pos[drudeParticles[i]] = RealVec(x[3*i], x[3*i+1], x[3*i+2]);
+    double energy = context.calcForcesAndEnergy(true, true);
+    for (int i = 0; i < numDrudeParticles; i++) {
+        RealVec f = force[drudeParticles[i]];
+        g[3*i] = -f[0];
+        g[3*i+1] = -f[1];
+        g[3*i+2] = -f[2];
+    }
+    return energy;
+}
+
+void ReferenceIntegrateDrudeSCFStepKernel::minimize(ContextImpl& context, double tolerance) {
+    // Record the initial positions and determine a normalization constant for scaling the tolerance.
+
+    vector<RealVec>& pos = extractPositions(context);
+    int numDrudeParticles = drudeParticles.size();
+    double norm = 0.0;
+    for (int i = 0; i < numDrudeParticles; i++) {
+        RealVec p = pos[drudeParticles[i]];
+        minimizerPos[3*i] = p[0];
+        minimizerPos[3*i+1] = p[1];
+        minimizerPos[3*i+2] = p[2];
+        norm += p.dot(p);
+    }
+    norm /= numDrudeParticles;
+    norm = (norm < 1 ? 1 : sqrt(norm));
+    minimizerParams.epsilon = tolerance/norm;
+    
+    // Perform the minimization.
+
+    lbfgsfloatval_t fx;
+    MinimizerData data(context, drudeParticles);
+    lbfgs(numDrudeParticles*3, minimizerPos, &fx, evaluate, NULL, &data, &minimizerParams);
 }
\ No newline at end of file

plugins/drude/platforms/reference/src/ReferenceDrudeKernels.h

@@ -35,6 +35,7 @@
 #include "ReferencePlatform.h"
 #include "openmm/DrudeKernels.h"
 #include "SimTKUtilities/RealVec.h"
+#include "lbfgs.h"
 #include <utility>
 #include <vector>
 
@@ -120,6 +121,47 @@ private:
     ReferenceConstraintAlgorithm* constraints;
 };
 
+/**
+ * This kernel is invoked by DrudeSCFIntegrator to take one time step
+ */
+class ReferenceIntegrateDrudeSCFStepKernel : public IntegrateDrudeSCFStepKernel {
+public:
+    ReferenceIntegrateDrudeSCFStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) :
+        IntegrateDrudeSCFStepKernel(name, platform), data(data), constraints(NULL), minimizerPos(NULL) {
+    }
+    ~ReferenceIntegrateDrudeSCFStepKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param integrator the DrudeSCFIntegrator this kernel will be used for
+     * @param force      the DrudeForce to get particle parameters from
+     */
+    void initialize(const System& system, const DrudeSCFIntegrator& integrator, const DrudeForce& force);
+    /**
+     * Execute the kernel.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param integrator     the DrudeSCFIntegrator this kernel is being used for
+     */
+    void execute(ContextImpl& context, const DrudeSCFIntegrator& integrator);
+    /**
+     * Compute the kinetic energy.
+     * 
+     * @param context     the context in which to execute this kernel
+     * @param integrator  the DrudeSCFIntegrator this kernel is being used for
+     */
+    double computeKineticEnergy(ContextImpl& context, const DrudeSCFIntegrator& integrator);
+private:
+    void minimize(ContextImpl& context, double tolerance);
+    ReferencePlatform::PlatformData& data;
+    std::vector<int> drudeParticles;
+    std::vector<double> particleInvMass;
+    ReferenceConstraintAlgorithm* constraints;
+    lbfgsfloatval_t *minimizerPos;
+    lbfgs_parameter_t minimizerParams;
+};
+
 } // namespace OpenMM
 
 #endif /*REFERENCE_DRUDE_KERNELS_H_*/

plugins/drude/platforms/reference/tests/TestReferenceDrudeSCFIntegrator.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) 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Reference implementation of DrudeSCFIntegrator.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/Platform.h"
+#include "openmm/System.h"
+#include "openmm/VirtualSite.h"
+#include "openmm/DrudeForce.h"
+#include "openmm/DrudeSCFIntegrator.h"
+#include "SimTKUtilities/SimTKOpenMMUtilities.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testWater() {
+    // Create a box of SWM4-NDP water molecules.  This involves constraints, virtual sites,
+    // and Drude particles.
+    
+    const int gridSize = 4;
+    const int numMolecules = gridSize*gridSize*gridSize;
+    const double spacing = 0.6;
+    const double boxSize = spacing*(gridSize+1);
+    System system;
+    NonbondedForce* nonbonded = new NonbondedForce();
+    DrudeForce* drude = new DrudeForce();
+    system.addForce(nonbonded);
+    system.addForce(drude);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nonbonded->setCutoffDistance(1.0);
+    for (int i = 0; i < numMolecules; i++) {
+        int startIndex = system.getNumParticles();
+        system.addParticle(15.6); // O
+        system.addParticle(0.4);  // D
+        system.addParticle(1.0);  // H1
+        system.addParticle(1.0);  // H2
+        system.addParticle(0.0);  // M
+        nonbonded->addParticle(1.71636, 0.318395, 0.21094*4.184);
+        nonbonded->addParticle(-1.71636, 1, 0);
+        nonbonded->addParticle(0.55733, 1, 0);
+        nonbonded->addParticle(0.55733, 1, 0);
+        nonbonded->addParticle(-1.11466, 1, 0);
+        for (int j = 0; j < 5; j++)
+            for (int k = 0; k < j; k++)
+                nonbonded->addException(startIndex+j, startIndex+k, 0, 1, 0);
+        system.addConstraint(startIndex, startIndex+2, 0.09572);
+        system.addConstraint(startIndex, startIndex+3, 0.09572);
+        system.addConstraint(startIndex+2, startIndex+3, 0.15139);
+        system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, 1.71636*1.71636/(100000*4.184), 1, 1);
+    }
+    vector<Vec3> positions;
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                Vec3 pos(i*spacing, j*spacing, k*spacing);
+                positions.push_back(pos);
+                positions.push_back(pos);
+                positions.push_back(pos+Vec3(0.09572, 0, 0));
+                positions.push_back(pos+Vec3(-0.023999, 0.092663, 0));
+                positions.push_back(pos);
+            }
+    
+    // Simulate it and check energy conservation and the total force on the Drude particles.
+    
+    DrudeSCFIntegrator integ(0.0005);
+    Platform& platform = Platform::getPlatformByName("Reference");
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+    context.applyConstraints(1e-5);
+    context.setVelocitiesToTemperature(300.0);
+    State state = context.getState(State::Energy);
+    double initialEnergy = state.getPotentialEnergy()+state.getKineticEnergy();
+    int numSteps = 1000;
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        state = context.getState(State::Energy | State::Forces);
+        ASSERT_EQUAL_TOL(initialEnergy, state.getPotentialEnergy()+state.getKineticEnergy(), 0.01);
+        const vector<Vec3>& force = state.getForces();
+        double norm = 0.0;
+        for (int j = 1; j < (int) force.size(); j += 5)
+            norm += sqrt(force[j].dot(force[j]));
+        norm = (norm/numMolecules);
+        ASSERT(norm < 1.0);
+    }
+}
+
+int main() {
+    try {
+        testWater();
+    }
+    catch(const std::exception& e) {
+        std::cout << "exception: " << e.what() << std::endl;
+        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
+        return 1;
+    }
+    std::cout << "Done" << std::endl;
+    return 0;
+}