Began writing reference implementation of Gay-Berne force

openmmapi/include/openmm/GayBerneForce.h

@@ -48,7 +48,7 @@ namespace OpenMM {
  * ellipsoid along each axis.  It also has three scale factors e_x, e_y, and e_z that scale the strength
  * of the interaction along each axis.  You can think of this force as a Lennard-Jones interaction computed
  * based on the distance between the nearest points on two ellipsoids.  If two particles each have all their
- * radii set to 0 and all their scale factors set 1, the interaction simplifies to a standard Lennard-Jones
+ * radii set to sigma/2 and all their scale factors set 1, the interaction simplifies to a standard Lennard-Jones
  * force between point particles.
  *
  * The orientation of a particle's ellipsoid is determined based on the positions of two other particles.
@@ -248,8 +248,9 @@ public:
      * This method has several limitations.  The only information it updates is the parameters of particles and exceptions.
      * All other aspects of the Force (the nonbonded method, the cutoff distance, etc.) are unaffected and can only be
      * changed by reinitializing the Context.  Furthermore, only the sigma and epsilon values of an exception can be
-     * changed; the pair of particles involved in the exception cannot change.  Finally, this method cannot be used to
-     * add new particles or exceptions, only to change the parameters of existing ones.
+     * changed; the pair of particles involved in the exception cannot change.  Likewise, the xparticle and yparticle
+     * defining the orientation of an ellipse cannot be changed.  Finally, this method cannot be used to add new
+     * particles or exceptions, only to change the parameters of existing ones.
      */
     void updateParametersInContext(Context& context);
     /**

openmmapi/src/GayBerneForce.cpp

@@ -83,6 +83,12 @@ void GayBerneForce::setSwitchingDistance(double distance) {
 }
 
 int GayBerneForce::addParticle(double sigma, double epsilon, int xparticle, int yparticle, double rx, double ry, double rz, double ex, double ey, double ez) {
+    if (yparticle == -1 && (ry != rz || ey != ez))
+        throw OpenMMException("GayBerneForce: yparticle is -1 for a particle that is not axially symmetric");
+    if (xparticle == -1 && (rx != rz || ex != ez))
+        throw OpenMMException("GayBerneForce: xparticle is -1 for a particle that is not spherical");
+    if (xparticle == -1 && yparticle != -1)
+        throw OpenMMException("GayBerneForce: xparticle cannot be -1 if yparticle is not also -1");
     particles.push_back(ParticleInfo(sigma, epsilon, xparticle, yparticle, rx, ry, rz, ex, ey, ez));
     return particles.size()-1;
 }
@@ -103,6 +109,12 @@ void GayBerneForce::getParticleParameters(int index, double& sigma, double& epsi
 
 void GayBerneForce::setParticleParameters(int index, double sigma, double epsilon, int xparticle, int yparticle, double rx, double ry, double rz, double ex, double ey, double ez) {
     ASSERT_VALID_INDEX(index, particles);
+    if (yparticle == -1 && (ry != rz || ey != ez))
+        throw OpenMMException("GayBerneForce: yparticle is -1 for a particle that is not axially symmetric");
+    if (xparticle == -1 && (rx != rz || ex != ez))
+        throw OpenMMException("GayBerneForce: xparticle is -1 for a particle that is not spherical");
+    if (xparticle == -1 && yparticle != -1)
+        throw OpenMMException("GayBerneForce: xparticle cannot be -1 if yparticle is not also -1");
     particles[index].sigma = sigma;
     particles[index].epsilon = epsilon;
     particles[index].xparticle = xparticle;

openmmapi/src/GayBerneForceImpl.cpp

@@ -73,12 +73,6 @@ void GayBerneForceImpl::initialize(ContextImpl& context) {
             msg << yparticle;
             throw OpenMMException(msg.str());
         }
-        if (yparticle == -1 && (ry != rz || ey != ez))
-            throw OpenMMException("GayBerneForce: yparticle is -1 for a particle that is not axially symmetric");
-        if (xparticle == -1 && (rx != rz || ex != ez))
-            throw OpenMMException("GayBerneForce: xparticle is -1 for a particle that is not spherical");
-        if (xparticle == -1 && yparticle != -1)
-            throw OpenMMException("GayBerneForce: xparticle cannot be -1 if yparticle is not also -1");
     }
     vector<set<int> > exceptions(owner.getNumParticles());
     for (int i = 0; i < owner.getNumExceptions(); i++) {
@@ -119,6 +113,7 @@ void GayBerneForceImpl::initialize(ContextImpl& context) {
 }
 
 double GayBerneForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+    if ((groups&(1<<owner.getForceGroup())) != 0)
         return kernel.getAs<CalcGayBerneForceKernel>().execute(context, includeForces, includeEnergy);
 }
 

platforms/reference/include/ReferenceGayBerneForce.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) 2016 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+#ifndef __ReferenceGayBerneForce_H__
+#define __ReferenceGayBerneForce_H__
+
+#include "openmm/GayBerneForce.h"
+#include "RealVec.h"
+#include <set>
+#include <utility>
+
+namespace OpenMM {
+
+class OPENMM_EXPORT ReferenceGayBerneForce {
+public:
+    /**
+     * Constructor.
+     */
+    ReferenceGayBerneForce(const GayBerneForce& force);
+
+    /**
+     * Compute the interaction.
+     *
+     * @param positions     the positions of the atoms
+     * @param forces        forces will be added to this vector
+     * @param boxVectors    the periodic box vectors
+     * @return the energy of the interaction
+     */
+    RealOpenMM calculateForce(const std::vector<RealVec>& positions, std::vector<RealVec>& forces, const RealVec* boxVectors);
+
+private:
+    struct ParticleInfo;
+    struct ExceptionInfo;
+    struct Matrix;
+    std::vector<ParticleInfo> particles;
+    std::vector<ExceptionInfo> exceptions;
+    std::set<std::pair<int, int> > exclusions;
+    GayBerneForce::NonbondedMethod nonbondedMethod;
+    RealOpenMM cutoffDistance, switchingDistance;
+    bool useSwitchingFunction;
+    std::vector<RealOpenMM> s;
+    std::vector<Matrix> A, B, G;
+
+    void computeEllipsoidFrames(const std::vector<RealVec>& positions);
+
+    RealOpenMM computeOneInteraction(int particle1, int particle2, RealOpenMM sigma, RealOpenMM epsilon, const std::vector<RealVec>& positions, std::vector<RealVec>& forces, const RealVec* boxVectors);
+};
+
+struct ReferenceGayBerneForce::ParticleInfo {
+    int xparticle, yparticle;
+    RealOpenMM sigma, epsilon, rx, ry, rz, ex, ey, ez;
+    bool radiiAreZero, scalesAreZero;
+};
+
+struct ReferenceGayBerneForce::ExceptionInfo {
+    int particle1, particle2;
+    RealOpenMM sigma, epsilon;
+};
+
+struct ReferenceGayBerneForce::Matrix {
+    RealOpenMM v[3][3];
+    RealVec operator*(const RealVec& r) {
+        return RealVec(v[0][0]*r[0] + v[0][1]*r[1] + v[0][2]*r[2],
+                       v[1][0]*r[0] + v[1][1]*r[1] + v[1][2]*r[2],
+                       v[2][0]*r[0] + v[2][1]*r[1] + v[2][2]*r[2]);
+    }
+
+    Matrix operator*(const Matrix& m) {
+        Matrix result;
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++) {
+                result.v[i][j] = 0;
+                for (int k = 0; k < 3; k++)
+                    result.v[i][j] += v[i][k]*m.v[k][j];
+            }
+        return result;
+    }
+
+    Matrix operator+(const Matrix& m) {
+        Matrix result;
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++)
+                result.v[i][j] = v[i][j]+m.v[i][j];
+        return result;
+    }
+
+    RealOpenMM determinant() {
+        return (v[0][0]*v[1][1]*v[2][2] + v[0][1]*v[1][2]*v[2][0] + v[0][2]*v[1][0]*v[2][1] -
+                v[0][0]*v[1][2]*v[2][1] - v[0][1]*v[1][0]*v[2][2] - v[0][2]*v[1][1]*v[2][0]);
+    }
+    
+    Matrix inverse() {
+        RealOpenMM invDet = 1/determinant();
+        Matrix result;
+        result.v[0][0] = invDet*(v[1][1]*v[2][2] - v[1][2]*v[2][1]);
+        result.v[1][0] = -invDet*(v[1][0]*v[2][2] - v[1][2]*v[2][0]);
+        result.v[2][0] = invDet*(v[1][0]*v[2][1] - v[1][1]*v[2][0]);
+        result.v[0][1] = -invDet*(v[0][1]*v[2][2] - v[0][2]*v[2][1]);
+        result.v[1][1] = invDet*(v[0][0]*v[2][2] - v[0][2]*v[2][0]);
+        result.v[2][1] = -invDet*(v[0][0]*v[2][1] - v[0][1]*v[2][0]);
+        result.v[0][2] = invDet*(v[0][1]*v[1][2] - v[0][2]*v[1][1]);
+        result.v[1][2] = -invDet*(v[0][0]*v[1][2] - v[0][2]*v[1][0]);
+        result.v[2][2] = invDet*(v[0][0]*v[1][1] - v[0][1]*v[1][0]);
+        return result;
+    }
+};
+
+} // namespace OpenMM
+
+#endif // __ReferenceGayBerneForce_H__

platforms/reference/include/ReferenceKernels.h

@@ -9,7 +9,7 @@
  * 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.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -48,6 +48,7 @@ class ReferenceCustomCentroidBondIxn;
 class ReferenceCustomCompoundBondIxn;
 class ReferenceCustomHbondIxn;
 class ReferenceCustomManyParticleIxn;
+class ReferenceGayBerneForce;
 class ReferenceBrownianDynamics;
 class ReferenceStochasticDynamics;
 class ReferenceConstraintAlgorithm;
@@ -940,6 +941,40 @@ private:
     NonbondedMethod nonbondedMethod;
 };
 
+/**
+ * This kernel is invoked by GayBerneForce to calculate the forces acting on the system.
+ */
+class ReferenceCalcGayBerneForceKernel : public CalcGayBerneForceKernel {
+public:
+    ReferenceCalcGayBerneForceKernel(std::string name, const Platform& platform) : CalcGayBerneForceKernel(name, platform), ixn(NULL) {
+    }
+    ~ReferenceCalcGayBerneForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the GayBerneForce this kernel will be used for
+     */
+    void initialize(const System& system, const GayBerneForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the GayBerneForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const GayBerneForce& force);
+private:
+    ReferenceGayBerneForce* ixn;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/reference/src/ReferenceKernelFactory.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -56,8 +56,6 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcCustomBondForceKernel(name, platform);
     if (name == CalcHarmonicAngleForceKernel::Name())
         return new ReferenceCalcHarmonicAngleForceKernel(name, platform);
-    if (name == CalcHarmonicAngleForceKernel::Name())
-        return new ReferenceCalcHarmonicAngleForceKernel(name, platform);
     if (name == CalcCustomAngleForceKernel::Name())
         return new ReferenceCalcCustomAngleForceKernel(name, platform);
     if (name == CalcPeriodicTorsionForceKernel::Name())
@@ -82,6 +80,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcCustomCompoundBondForceKernel(name, platform);
     if (name == CalcCustomManyParticleForceKernel::Name())
         return new ReferenceCalcCustomManyParticleForceKernel(name, platform);
+    if (name == CalcGayBerneForceKernel::Name())
+        return new ReferenceCalcGayBerneForceKernel(name, platform);
     if (name == IntegrateVerletStepKernel::Name())
         return new ReferenceIntegrateVerletStepKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/reference/src/ReferenceKernels.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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -49,6 +49,7 @@
 #include "ReferenceCustomNonbondedIxn.h"
 #include "ReferenceCustomManyParticleIxn.h"
 #include "ReferenceCustomTorsionIxn.h"
+#include "ReferenceGayBerneForce.h"
 #include "ReferenceHarmonicBondIxn.h"
 #include "ReferenceLJCoulomb14.h"
 #include "ReferenceLJCoulombIxn.h"
@@ -1860,6 +1861,25 @@ void ReferenceCalcCustomManyParticleForceKernel::copyParametersToContext(Context
     }
 }
 
+ReferenceCalcGayBerneForceKernel::~ReferenceCalcGayBerneForceKernel() {
+    if (ixn != NULL)
+        delete ixn;
+}
+
+void ReferenceCalcGayBerneForceKernel::initialize(const System& system, const GayBerneForce& force) {
+    ixn = new ReferenceGayBerneForce(force);
+}
+
+double ReferenceCalcGayBerneForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    return ixn->calculateForce(extractPositions(context), extractForces(context), extractBoxVectors(context));
+}
+
+void ReferenceCalcGayBerneForceKernel::copyParametersToContext(ContextImpl& context, const GayBerneForce& force) {
+    delete ixn;
+    ixn = NULL;
+    ixn = new ReferenceGayBerneForce(force);
+}
+
 ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/reference/src/ReferencePlatform.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -64,6 +64,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcCustomCentroidBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
+    registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceGayBerneForce.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) 2016 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 "ReferenceGayBerneForce.h"
+#include "ReferenceForce.h"
+#include "openmm/OpenMMException.h"
+#include <cmath>
+
+using namespace OpenMM;
+using namespace std;
+
+ReferenceGayBerneForce::ReferenceGayBerneForce(const GayBerneForce& force) {
+    // Record the force parameters.
+
+    int numParticles = force.getNumParticles();
+    particles.resize(numParticles);
+    for (int i = 0; i < numParticles; i++) {
+        ParticleInfo& p = particles[i];
+        force.getParticleParameters(i, p.sigma, p.epsilon, p.xparticle, p.yparticle, p.rx, p.ry, p.rz, p.ex, p.ey, p.ez);
+        p.radiiAreZero = (p.rx == 0 && p.ry == 0 && p.rz == 0);
+        p.scalesAreZero = (p.ex == 0 && p.ey == 0 && p.ez == 0);
+    }
+    int numExceptions = force.getNumExceptions();
+    exceptions.resize(numExceptions);
+    for (int i = 0; i < numExceptions; i++) {
+        ExceptionInfo& e = exceptions[i];
+        force.getExceptionParameters(i, e.particle1, e.particle2, e.sigma, e.epsilon);
+        exclusions.insert(make_pair(min(e.particle1, e.particle2), max(e.particle1, e.particle2)));
+    }
+    nonbondedMethod = force.getNonbondedMethod();
+    cutoffDistance = force.getCutoffDistance();
+    switchingDistance = force.getSwitchingDistance();
+    useSwitchingFunction = force.getUseSwitchingFunction();
+
+    // Allocate workspace for calculations.
+
+    s.resize(numParticles);
+    A.resize(numParticles);
+    B.resize(numParticles);
+    G.resize(numParticles);
+
+    // We can precompute the shape factors.
+
+    for (int i = 0; i < numParticles; i++) {
+        ParticleInfo& p = particles[i];
+        s[i] = (p.rx*p.ry + p.rz*p.rz)*SQRT(p.rx*p.ry);
+    }
+}
+
+RealOpenMM ReferenceGayBerneForce::calculateForce(const vector<RealVec>& positions, vector<RealVec>& forces, const RealVec* boxVectors) {
+    if (nonbondedMethod == GayBerneForce::CutoffPeriodic) {
+        double minAllowedSize = 1.999999*cutoffDistance;
+        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
+            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
+    }
+
+    // First find the orientations of the particles and compute the matrices we'll be needing.
+
+    computeEllipsoidFrames(positions);
+
+    // Compute standard interactions.
+
+    RealOpenMM energy = 0;
+    int numParticles = particles.size();
+    for (int i = 1; i < numParticles; i++)
+        for (int j = 0; j < i; j++) {
+            if (exclusions.find(make_pair(j, i)) != exclusions.end())
+                continue; // This interaction will be handled by an exception.
+            RealOpenMM sigma = 0.5*(particles[i].sigma+particles[j].sigma);
+            RealOpenMM epsilon = SQRT(particles[i].epsilon*particles[j].epsilon);
+            energy += computeOneInteraction(i, j, sigma, epsilon, positions, forces, boxVectors);
+        }
+
+    // Compute exceptions.
+
+    int numExceptions = exceptions.size();
+    for (int i = 0; i < numExceptions; i++) {
+        ExceptionInfo& e = exceptions[i];
+        energy += computeOneInteraction(e.particle1, e.particle2, e.sigma, e.epsilon, positions, forces, boxVectors);
+    }
+    return energy;
+}
+
+void ReferenceGayBerneForce::computeEllipsoidFrames(const vector<RealVec>& positions) {
+    int numParticles = particles.size();
+    for (int particle = 0; particle < numParticles; particle++) {
+        ParticleInfo& p = particles[particle];
+
+        // Compute the local coordinate system of the ellipsoid;
+
+        RealVec xdir, ydir, zdir;
+        if (p.xparticle == -1) {
+            xdir = RealVec(1, 0, 0);
+            ydir = RealVec(0, 1, 0);
+        }
+        else {
+            xdir = positions[particle]-positions[p.xparticle];
+            xdir /= SQRT(xdir.dot(xdir));
+            if (p.yparticle == -1) {
+                if (xdir[1] > -0.5 && xdir[1] < 0.5)
+                    ydir = RealVec(0, 1, 0);
+                else
+                    ydir = RealVec(1, 0, 0);
+            }
+            else
+                ydir = positions[particle]-positions[p.yparticle];
+            ydir -= xdir*(xdir.dot(ydir));
+            ydir /= SQRT(ydir.dot(ydir));
+        }
+        zdir = xdir.cross(ydir);
+
+        // Compute matrices we will need later.
+
+        RealOpenMM (&a)[3][3] = A[particle].v;
+        RealOpenMM (&b)[3][3] = B[particle].v;
+        RealOpenMM (&g)[3][3] = G[particle].v;
+        a[0][0] = xdir[0];
+        a[0][1] = xdir[1];
+        a[0][2] = xdir[2];
+        a[1][0] = ydir[0];
+        a[1][1] = ydir[1];
+        a[1][2] = ydir[2];
+        a[2][0] = zdir[0];
+        a[2][1] = zdir[1];
+        a[2][2] = zdir[2];
+        RealVec r2(p.rx*p.rx, p.ry*p.ry, p.rz*p.rz);
+        RealVec e2(p.ex*p.ex, p.ey*p.ey, p.ez*p.ez);
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++) {
+                b[i][j] = 0;
+                g[i][j] = 0;
+                for (int k = 0; k < 3; k++) {
+                    b[i][j] += a[k][i]*e2[k]*a[k][j];
+                    g[i][j] += a[k][i]*r2[k]*a[k][j];
+                }
+            }
+    }
+}
+
+RealOpenMM ReferenceGayBerneForce::computeOneInteraction(int particle1, int particle2, RealOpenMM sigma, RealOpenMM epsilon, const vector<RealVec>& positions, vector<RealVec>& forces, const RealVec* boxVectors) {
+    // Compute the displacement and check against the cutoff.
+
+    RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+    if (nonbondedMethod == GayBerneForce::CutoffPeriodic)
+        ReferenceForce::getDeltaRPeriodic(positions[particle2], positions[particle1], boxVectors, deltaR);
+    else
+        ReferenceForce::getDeltaR(positions[particle2], positions[particle1], deltaR);
+    RealOpenMM dist = deltaR[ReferenceForce::RIndex];
+    if (nonbondedMethod != GayBerneForce::NoCutoff && dist >= cutoffDistance)
+        return 0;
+
+    // Compute vectors and matrices we'll be needing.
+
+    RealVec dr(deltaR[ReferenceForce::XIndex], deltaR[ReferenceForce::YIndex], deltaR[ReferenceForce::ZIndex]);
+    RealVec drUnit = dr/dist;
+    Matrix B12 = B[particle1]+B[particle2];
+    Matrix G12 = G[particle1]+G[particle2];
+    Matrix B12inv = B12.inverse();
+    Matrix G12inv = G12.inverse();
+
+    // Estimate the distance between the ellipsoids and compute the first term in the energy.
+
+    ParticleInfo& p1 = particles[particle1];
+    ParticleInfo& p2 = particles[particle2];
+    RealOpenMM h12 = dist;
+    if (!p1.radiiAreZero || !p2.radiiAreZero)
+        h12 -= 1/SQRT(0.5*drUnit.dot(G12inv*drUnit));
+    RealOpenMM rho = sigma/(h12+sigma);
+    RealOpenMM rho2 = rho*rho;
+    RealOpenMM rho6 = rho2*rho2*rho2;
+    RealOpenMM u = 4*epsilon*(rho6*rho6-rho6);
+
+    // Compute the second term in the energy.
+
+    RealOpenMM eta = SQRT(2*s[particle1]*s[particle2]/G12.determinant());
+
+    // Compute the third term in the energy.
+
+    RealOpenMM chi = 2*drUnit.dot(B12inv*drUnit);
+    chi *= chi;
+    return u*eta*chi;
+}

platforms/reference/tests/TestReferenceGayBerneForce.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) 2016 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 "TestGayBerneForce.h"
+
+void runPlatformTests() {
+}

tests/TestGayBerneForce.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) 2016 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+#ifdef WIN32
+  #define _USE_MATH_DEFINES // Needed to get M_PI
+#endif
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "openmm/GayBerneForce.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testPointParticles() {
+    // For point particles, it should be identical to a standard Lennard-Jones force.
+
+    const int numParticles = 10;
+    const double sigma = 0.5;
+    const double epsilon = 1.5;
+    System system;
+    GayBerneForce* gb = new GayBerneForce();
+    NonbondedForce* nb = new NonbondedForce();
+    system.addForce(gb);
+    system.addForce(nb);
+    gb->setForceGroup(1);
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        gb->addParticle(sigma, epsilon, -1, -1, sigma/2, sigma/2, sigma/2, 1, 1, 1);
+        nb->addParticle(0, sigma, epsilon);
+        positions.push_back(Vec3(2.0*genrand_real2(sfmt), 2.0*genrand_real2(sfmt), 2.0*genrand_real2(sfmt)));
+    }
+    VerletIntegrator integ(0.001);
+
+    // Compute forces and energy with each one and compare them.
+
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+    State state1 = context.getState(State::Forces | State::Energy, false, 1);
+    State state2 = context.getState(State::Forces | State::Energy, false, 2);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
+}
+
+void runPlatformTests();
+
+int main(int argc, char* argv[]) {
+    try {
+        initializeTests(argc, argv);
+        testPointParticles();
+        runPlatformTests();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}