Merge pull request #5 from peastman/master

Fixes to a couple of bugs (in CustomIntegrator and DrudeLangevinIntegrator)

platforms/cuda/src/CudaKernels.cpp

@@ -4554,7 +4554,7 @@ void CudaIntegrateCustomStepKernel::initialize(const System& system, const Custo
     numGlobalVariables = integrator.getNumGlobalVariables();
     int elementSize = (cu.getUseDoublePrecision() || cu.getUseMixedPrecision() ? sizeof(double) : sizeof(float));
     globalValues = new CudaArray(cu, max(1, numGlobalVariables), elementSize, "globalVariables");
-    sumBuffer = new CudaArray(cu, 3*system.getNumParticles(), elementSize, "sumBuffer");
+    sumBuffer = new CudaArray(cu, ((3*system.getNumParticles()+3)/4)*4, elementSize, "sumBuffer");
     potentialEnergy = new CudaArray(cu, 1, cu.getEnergyBuffer().getElementSize(), "potentialEnergy");
     kineticEnergy = new CudaArray(cu, 1, elementSize, "kineticEnergy");
     perDofValues = new CudaParameterSet(cu, integrator.getNumPerDofVariables(), 3*system.getNumParticles(), "perDofVariables", false, cu.getUseDoublePrecision() || cu.getUseMixedPrecision());
@@ -5103,6 +5103,7 @@ void CudaIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegrat
             kernelArgs[i][0][10] = &randomIndex;
             if (requiredUniform[i] > 0)
                 cu.executeKernel(randomKernel, &randomArgs[0], numAtoms);
+            cu.clearBuffer(*sumBuffer);
             cu.executeKernel(kernels[i][0], &kernelArgs[i][0][0], numAtoms);
             cu.executeKernel(kernels[i][1], &kernelArgs[i][1][0], CudaContext::ThreadBlockSize, CudaContext::ThreadBlockSize);
         }
@@ -5150,6 +5151,7 @@ double CudaIntegrateCustomStepKernel::computeKineticEnergy(ContextImpl& context,
     int randomIndex = 0;
     kineticEnergyArgs[1] = &posCorrection;
     kineticEnergyArgs[10] = &randomIndex;
+    cu.clearBuffer(*sumBuffer);
     cu.executeKernel(kineticEnergyKernel, &kineticEnergyArgs[0], cu.getNumAtoms());
     void* args[] = {&sumBuffer->getDevicePointer(), &kineticEnergy->getDevicePointer()};
     cu.executeKernel(sumKineticEnergyKernel, args, CudaContext::ThreadBlockSize, CudaContext::ThreadBlockSize);

platforms/cuda/tests/TestCudaCustomIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -345,20 +345,20 @@ void testSum() {
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
     for (int i = 0; i < numParticles; i++) {
-        system.addParticle(1.5);
-        nb->addParticle(i%2 == 0 ? 1 : -1, 0.1, 1);
+        system.addParticle(i%10 == 0 ? 0.0 : 1.5);
+        nb->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.1, 1);
         bool close = true;
         while (close) {
             positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
             close = false;
             for (int j = 0; j < i; ++j) {
                 Vec3 delta = positions[i]-positions[j];
-                if (delta.dot(delta) < 0.1)
+                if (delta.dot(delta) < 1)
                     close = true;
             }
         }
     }
-    CustomIntegrator integrator(0.01);
+    CustomIntegrator integrator(0.005);
     integrator.addGlobalVariable("ke", 0);
     integrator.addComputePerDof("v", "v+dt*f/m");
     integrator.addComputePerDof("x", "x+dt*v");
@@ -368,10 +368,8 @@ void testSum() {
     
     // See if the sum is being computed correctly.
     
-    State state = context.getState(State::Energy);
-    const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
     for (int i = 0; i < 100; ++i) {
-        state = context.getState(State::Energy);
+        State state = context.getState(State::Energy);
         ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5);
         integrator.step(1);
     }

platforms/opencl/src/OpenCLKernels.cpp

@@ -4786,7 +4786,7 @@ void OpenCLIntegrateCustomStepKernel::initialize(const System& system, const Cus
     numGlobalVariables = integrator.getNumGlobalVariables();
     int elementSize = (cl.getUseDoublePrecision() || cl.getUseMixedPrecision() ? sizeof(double) : sizeof(float));
     globalValues = new OpenCLArray(cl, max(1, numGlobalVariables), elementSize, "globalVariables");
-    sumBuffer = new OpenCLArray(cl, 3*system.getNumParticles(), elementSize, "sumBuffer");
+    sumBuffer = new OpenCLArray(cl, ((3*system.getNumParticles()+3)/4)*4, elementSize, "sumBuffer");
     potentialEnergy = new OpenCLArray(cl, 1, cl.getEnergyBuffer().getElementSize(), "potentialEnergy");
     kineticEnergy = new OpenCLArray(cl, 1, elementSize, "kineticEnergy");
     perDofValues = new OpenCLParameterSet(cl, integrator.getNumPerDofVariables(), 3*system.getNumParticles(), "perDofVariables", false, cl.getUseDoublePrecision() || cl.getUseMixedPrecision());
@@ -5329,6 +5329,7 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
             kernels[i][0].setArg<cl_uint>(10, integration.prepareRandomNumbers(requiredGaussian[i]));
             if (requiredUniform[i] > 0)
                 cl.executeKernel(randomKernel, numAtoms);
+            cl.clearBuffer(*sumBuffer);
             cl.executeKernel(kernels[i][0], numAtoms);
             cl.executeKernel(kernels[i][1], OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize);
         }
@@ -5375,6 +5376,7 @@ double OpenCLIntegrateCustomStepKernel::computeKineticEnergy(ContextImpl& contex
             cl.executeKernel(sumPotentialEnergyKernel, OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize);
         forcesAreValid = true;
     }
+    cl.clearBuffer(*sumBuffer);
     cl.executeKernel(kineticEnergyKernel, cl.getNumAtoms());
     cl.executeKernel(sumKineticEnergyKernel, OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize);
     if (cl.getUseDoublePrecision() || cl.getUseMixedPrecision()) {

platforms/opencl/tests/TestOpenCLCustomIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2011 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -345,20 +345,20 @@ void testSum() {
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
     for (int i = 0; i < numParticles; i++) {
-        system.addParticle(1.5);
-        nb->addParticle(i%2 == 0 ? 1 : -1, 0.1, 1);
+        system.addParticle(i%10 == 0 ? 0.0 : 1.5);
+        nb->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.1, 1);
         bool close = true;
         while (close) {
             positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
             close = false;
             for (int j = 0; j < i; ++j) {
                 Vec3 delta = positions[i]-positions[j];
-                if (delta.dot(delta) < 0.1)
+                if (delta.dot(delta) < 1)
                     close = true;
             }
         }
     }
-    CustomIntegrator integrator(0.01);
+    CustomIntegrator integrator(0.005);
     integrator.addGlobalVariable("ke", 0);
     integrator.addComputePerDof("v", "v+dt*f/m");
     integrator.addComputePerDof("x", "x+dt*v");
@@ -368,10 +368,8 @@ void testSum() {
     
     // See if the sum is being computed correctly.
     
-    State state = context.getState(State::Energy);
-    const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
     for (int i = 0; i < 100; ++i) {
-        state = context.getState(State::Energy);
+        State state = context.getState(State::Energy);
         ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5);
         integrator.step(1);
     }

platforms/reference/tests/TestReferenceCustomIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2011 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -336,20 +336,20 @@ void testSum() {
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
     for (int i = 0; i < numParticles; i++) {
-        system.addParticle(1.5);
-        nb->addParticle(i%2 == 0 ? 1 : -1, 0.1, 1);
+        system.addParticle(i%10 == 0 ? 0.0 : 1.5);
+        nb->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.1, 1);
         bool close = true;
         while (close) {
             positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
             close = false;
             for (int j = 0; j < i; ++j) {
                 Vec3 delta = positions[i]-positions[j];
-                if (delta.dot(delta) < 0.1)
+                if (delta.dot(delta) < 1)
                     close = true;
             }
         }
     }
-    CustomIntegrator integrator(0.01);
+    CustomIntegrator integrator(0.005);
     integrator.addGlobalVariable("ke", 0);
     integrator.addComputePerDof("v", "v+dt*f/m");
     integrator.addComputePerDof("x", "x+dt*v");
@@ -359,10 +359,8 @@ void testSum() {
     
     // See if the sum is being computed correctly.
     
-    State state = context.getState(State::Energy);
-    const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
     for (int i = 0; i < 100; ++i) {
-        state = context.getState(State::Energy);
+        State state = context.getState(State::Energy);
         ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5);
         integrator.step(1);
     }

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/opencl/src/kernels/drudeLangevin.cl

@@ -13,10 +13,12 @@ __kernel void integrateDrudeLangevinPart1(__global mixed4* restrict velm, __glob
         int index = normalParticles[i];
         mixed4 velocity = velm[index];
         if (velocity.w != 0) {
-            mixed sqrtInvMass = SQRT(velocity.w);
+            mixed sqrtInvMass = sqrt(velocity.w);
             float4 rand = random[randomIndex+index];
             real4 f = force[index];
-            velocity.xyz = vscale*velocity.xyz + fscale*velocity.w*f.xyz + noisescale*sqrtInvMass*rand.xyz;
+            velocity.x = vscale*velocity.x + fscale*velocity.w*f.x + noisescale*sqrtInvMass*rand.x;
+            velocity.y = vscale*velocity.y + fscale*velocity.w*f.y + noisescale*sqrtInvMass*rand.y;
+            velocity.z = vscale*velocity.z + fscale*velocity.w*f.z + noisescale*sqrtInvMass*rand.z;
             velm[index] = velocity;
             posDelta[index] = (mixed4) (stepSize*velocity.x, stepSize*velocity.y, stepSize*velocity.z, 0);
         }
@@ -29,24 +31,28 @@ __kernel void integrateDrudeLangevinPart1(__global mixed4* restrict velm, __glob
         int2 particles = pairParticles[i];
         mixed4 velocity1 = velm[particles.x];
         mixed4 velocity2 = velm[particles.y];
-        mixed mass1 = RECIP(velocity1.w);
-        mixed mass2 = RECIP(velocity2.w);
-        mixed invTotalMass = RECIP(mass1+mass2);
+        mixed mass1 = 1/velocity1.w;
+        mixed mass2 = 1/velocity2.w;
+        mixed invTotalMass = 1/(mass1+mass2);
         mixed invReducedMass = (mass1+mass2)*velocity1.w*velocity2.w;
         mixed mass1fract = invTotalMass*mass1;
         mixed mass2fract = invTotalMass*mass2;
-        mixed sqrtInvTotalMass = SQRT(invTotalMass);
-        mixed sqrtInvReducedMass = SQRT(invReducedMass);
+        mixed sqrtInvTotalMass = sqrt(invTotalMass);
+        mixed sqrtInvReducedMass = sqrt(invReducedMass);
         mixed4 cmVel = velocity1*mass1fract+velocity2*mass2fract;
         mixed4 relVel = velocity2-velocity1;
-        mixed4 force1 = force[particles.x];
-        mixed4 force2 = force[particles.y];
+        mixed4 force1 = convert_mixed4(force[particles.x]);
+        mixed4 force2 = convert_mixed4(force[particles.y]);
         mixed4 cmForce = force1+force2;
         mixed4 relForce = force2*mass1fract - force1*mass2fract;
         float4 rand1 = random[randomIndex+2*i];
         float4 rand2 = random[randomIndex+2*i+1];
-        cmVel.xyz = vscale*cmVel.xyz + fscale*invTotalMass*cmForce.xyz + noisescale*sqrtInvTotalMass*rand1.xyz;
-        relVel.xyz = vscaleDrude*relVel.xyz + fscaleDrude*invReducedMass*relForce.xyz + noisescaleDrude*sqrtInvReducedMass*rand2.xyz;
+        cmVel.x = vscale*cmVel.x + fscale*invTotalMass*cmForce.x + noisescale*sqrtInvTotalMass*rand1.x;
+        cmVel.y = vscale*cmVel.y + fscale*invTotalMass*cmForce.y + noisescale*sqrtInvTotalMass*rand1.y;
+        cmVel.z = vscale*cmVel.z + fscale*invTotalMass*cmForce.z + noisescale*sqrtInvTotalMass*rand1.z;
+        relVel.x = vscaleDrude*relVel.x + fscaleDrude*invReducedMass*relForce.x + noisescaleDrude*sqrtInvReducedMass*rand2.x;
+        relVel.y = vscaleDrude*relVel.y + fscaleDrude*invReducedMass*relForce.y + noisescaleDrude*sqrtInvReducedMass*rand2.y;
+        relVel.z = vscaleDrude*relVel.z + fscaleDrude*invReducedMass*relForce.z + noisescaleDrude*sqrtInvReducedMass*rand2.z;
         velocity1.xyz = cmVel.xyz-relVel.xyz*mass2fract;
         velocity2.xyz = cmVel.xyz+relVel.xyz*mass1fract;
         velm[particles.x] = velocity1;

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);
+}
+
+ReferenceIntegrateDrudeSCFStepKernel::~ReferenceIntegrateDrudeSCFStepKernel() {
+    if (constraints != NULL)
+        delete constraints;
+    if (minimizerPos != NULL)
+        lbfgs_free(minimizerPos);
+}
+
+void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, const DrudeSCFIntegrator& integrator, const DrudeForce& force) {
+    // Identify Drude particles.
     
-    // Compute the shifted velocities.
+    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);
+    }
+
+    // Record particle masses.
+
+    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);
+    }
     
-    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];
+    // Prepare constraints.
+    
+    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());
     }
     
-    // Apply constraints to them.
+    // Initialize the energy minimizer.
     
-    if (constraints != NULL) {
-        constraints->setTolerance(1e-4);
-        constraints->applyToVelocities(numParticles, posData, shiftedVel, particleInvMass);
+    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.
     
-    // Compute the kinetic energy.
+    if (constraints != NULL)
+        constraints->apply(numParticles, pos, xPrime, particleInvMass);
     
-    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;
+    // 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;
+}