Merge pull request #1553 from peastman/gayberne

Implemented Gay-Berne ellipsoids

docs-source/usersguide/references.bib

@@ -80,6 +80,17 @@
    type = {Journal Article}
 }
 
+@article{Everaers2003
+   author = {Everaers, R. and Ejtehadi, M. R.},
+   title = {Interaction potentials for soft and hard ellipsoids},
+   journal = {Physical Review E},
+   volume = {67},
+   issue = {4},
+   pages = {041710},
+   year = {2003},
+   type = {Journal Article}
+}
+
 @article{Hall1984
    author = {Hall, Randall W. and Berne, B. J.},
    title = {Nonergodicity in path integral molecular dynamics},

docs-source/usersguide/theory.rst

@@ -182,7 +182,7 @@ an energy term of the form
 
 
 .. math::
-   E=4\epsilon\left({\left(\frac{\sigma}{r}\right)}^{\text{12}}-{\left(\frac{\sigma}{r}\right)}^{6}\right)
+   E=4\epsilon\left({\left(\frac{\sigma}{r}\right)}^{12}-{\left(\frac{\sigma}{r}\right)}^{6}\right)
 
 
 where *r* is the distance between the two particles, :math:`\sigma` is the distance
@@ -201,7 +201,7 @@ at the cutoff distance.  When :math:`r_\mathit{switch} < r < r_\mathit{cutoff}`\
 
 where :math:`x = (r-r_\mathit{switch})/(r_\mathit{cutoff}-r_\mathit{switch})`. This function decreases smoothly from 1 at
 :math:`r = r_\mathit{switch}` to 0 at :math:`r = r_\mathit{cutoff}`, and has continuous first and
-second derivatives at both ends
+second derivatives at both ends.
 
 When an exception has been added for a pair of particles, :math:`\sigma` and :math:`\epsilon`
 are the parameters specified by the exception.  Otherwise they are determined
@@ -497,6 +497,84 @@ its atomic radius, and :math:`r_\mathit{solvent}` is the solvent radius, which i
 to be 0.14 nm.  The default value for the energy scale :math:`E_{SA}` is 2.25936 kJ/mol/nm\ :sup:`2`\ .
 
 
+GayBerneForce
+*************
+
+This is similar to the Lennard-Jones interaction described in section :ref:`lennard-jones-interaction`,
+but instead of being based on the distance between two point particles, it is based
+on the distance of closest approach between two ellipsoids.\ :cite:`Everaers2003`
+Let :math:`\mathbf{A}_1` and :math:`\mathbf{A}_2` be rotation matrices that transform
+from the lab frame to the body frames of two interacting ellipsoids.  These rotations
+are determined from the positions of other particles, as described in the API documentation.
+Let :math:`\mathbf{r}_{12}` be the vector pointing from particle 1 to particle 2, and
+:math:`\hat{\mathbf{r}}_{12}=\mathbf{r}_{12}/|\mathbf{r}_{12}|`.  Let :math:`\mathbf{S}_1`
+and :math:`\mathbf{S}_2` be diagonal matrices containing the three radii of each particle:
+
+.. math::
+   \mathbf{S}_i=\begin{bmatrix}
+   a_i & 0 & 0 \\
+   0 & b_i & 0 \\
+   0 & 0 & c_i
+   \end{bmatrix}
+
+The energy is computed as a product of three terms:
+
+.. math::
+   E=U_r(\mathbf{A}_1, \mathbf{A}_2, \mathbf{r}_{12}) \cdot \eta_{12}(\mathbf{A}_1, \mathbf{A}_2) \cdot \chi_{12}(\mathbf{A}_1, \mathbf{A}_2, \hat{\mathbf{r}}_{12}) 
+
+The first term describes the distance dependence, and is very similar in form to
+the Lennard-Jones interaction:
+
+.. math::
+   U_r=4\epsilon\left({\left(\frac{\sigma}{h_{12}+\sigma}\right)}^{12}-{\left(\frac{\sigma}{h_{12}+\sigma}\right)}^{6}\right)
+
+where :math:`h_{12}` is an approximation to the distance of closest approach between
+the two ellipsoids:
+
+.. math::
+   h_{12}=|\mathbf{r}_{12}|-\sigma_{12}(\mathbf{A}_1, \mathbf{A}_2, \hat{\mathbf{r}}_{12})
+
+.. math::
+   \sigma_{12}(\mathbf{A}_1, \mathbf{A}_2, \hat{\mathbf{r}}_{12})=\left[ \frac{1}{2} \hat{\mathbf{r}}_{12}^T \mathbf{G}_{12}^{-1} \hat{\mathbf{r}}_{12} \right]^{-1/2}
+
+.. math::
+   \mathbf{G}_{12}=\mathbf{A}_1^T \mathbf{S}_1^2 \mathbf{A}_1 + \mathbf{A}_2^T \mathbf{S}_2^2 \mathbf{A}_2
+
+The second term adjusts the energy based on the relative orientations of the two ellipsoids:
+
+.. math::
+   \eta_{12}(\mathbf{A}_1, \mathbf{A}_2)=\left[ \frac{2 s_1 s_2}{\text{det}(\mathbf{G}_{12})} \right]^{1/2}
+
+.. math::
+   s_i=(a_i b_i + c_i^2)\sqrt{a_i b_i}
+
+The third term applies the user-defined scale factors :math:`e_a`, :math:`e_b`,
+and :math:`e_c` that adjust the strength of the interaction along each axis:
+
+.. math::
+   \chi_{12}(\mathbf{A}_1, \mathbf{A}_2, \hat{\mathbf{r}}_{12})=(2 \hat{\mathbf{r}}_{12}^T \mathbf{B}_{12}^{-1} \hat{\mathbf{r}}_{12})^2
+
+.. math::
+   \mathbf{B}_{12}=\mathbf{A}_1^T \mathbf{E}_1 \mathbf{A}_1 + \mathbf{A}_2^T \mathbf{E}_2 \mathbf{A}_2
+
+.. math::
+   \mathbf{E}_i=\begin{bmatrix}
+   e_{ai}^{-1/2} & 0 & 0 \\
+   0 & e_{bi}^{-1/2} & 0 \\
+   0 & 0 & e_{ci}^{-1/2}
+   \end{bmatrix}
+
+When using a cutoff, you can optionally use a switching function to make the energy go smoothly to 0
+at the cutoff distance.  When :math:`r_\mathit{switch} < r < r_\mathit{cutoff}`\ , the energy is multiplied by
+
+.. math::
+   S=1-{6x}^{5}+15{x}^{4}-10{x}^{3}
+
+where :math:`x = (r-r_\mathit{switch})/(r_\mathit{cutoff}-r_\mathit{switch})`. This function decreases smoothly from 1 at
+:math:`r = r_\mathit{switch}` to 0 at :math:`r = r_\mathit{cutoff}`, and has continuous first and
+second derivatives at both ends.
+
+
 AndersenThermostat
 ******************
 

olla/include/openmm/kernels.h

@@ -47,6 +47,7 @@
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/CustomManyParticleForce.h"
 #include "openmm/CustomTorsionForce.h"
+#include "openmm/GayBerneForce.h"
 #include "openmm/GBSAOBCForce.h"
 #include "openmm/HarmonicAngleForce.h"
 #include "openmm/HarmonicBondForce.h"
@@ -897,6 +898,41 @@ public:
     virtual void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force) = 0;
 };
 
+/**
+ * This kernel is invoked by GayBerneForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcGayBerneForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcGayBerneForce";
+    }
+    CalcGayBerneForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the GayBerneForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const GayBerneForce& force) = 0;
+    /**
+     * 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
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the GayBerneForce to copy the parameters from
+     */
+    virtual void copyParametersToContext(ContextImpl& context, const GayBerneForce& force) = 0;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

openmmapi/include/OpenMM.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) 2009-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -49,6 +49,7 @@
 #include "openmm/CustomManyParticleForce.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/Force.h"
+#include "openmm/GayBerneForce.h"
 #include "openmm/GBSAOBCForce.h"
 #include "openmm/HarmonicAngleForce.h"
 #include "openmm/HarmonicBondForce.h"

openmmapi/include/openmm/GayBerneForce.h

+#ifndef OPENMM_GAYBERNEFORCE_H_
+#define OPENMM_GAYBERNEFORCE_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.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "Context.h"
+#include "Force.h"
+#include <map>
+#include <utility>
+#include <vector>
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This class implements the Gay-Berne anisotropic potential.  This is similar to a Lennard-Jones potential,
+ * but it represents the particles as ellipsoids rather than point particles.  In addition to the standard
+ * sigma and epsilon parameters, each particle has three widths sx, sy, and sz that give the diameter of the
+ * ellipsoid along each axis.  It also has three scale factors ex, ey, and ez 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.  The scale factors act as multipliers
+ * for epsilon along each axis, so the strength of the interaction along the ellipsoid's x axis is multiplied by
+ * ex, and likewise for the other axes.  If two particles each have all their widths set to sigma and all their
+ * scale factors set to 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.
+ * The vector to the first particle sets the direction of the x axis.  The vector to the second particle
+ * (after subtracting out any x component) sets the direction of the y axis.  If the ellipsoid is axially
+ * symmetric (sy=sz and ey=ez), you can omit the second particle and define only an x axis direction.
+ * If the ellipsoid is a sphere (all three widths and all three scale factors are equal), both particles
+ * can be omitted.
+ *
+ * To determine the values of sigma and epsilon for an interaction, this class uses Lorentz-Berthelot
+ * combining rules: it takes the arithmetic mean of the sigmas and the geometric mean of the epsilons for
+ * the two interacting particles.  You also can specify "exceptions", particular pairs of particles for
+ * which different values should be used.
+ *
+ * To use this class, create a GayBerneForce object, then call addParticle() once for each particle in the
+ * System to define its parameters.  The number of particles for which you define parameters must be exactly
+ * equal to the number of particles in the System, or else an exception will be thrown when you try to
+ * create a Context.  After a particle has been added, you can modify its force field parameters by calling
+ * setParticleParameters().  This will have no effect on Contexts that already exist unless you call
+ * updateParametersInContext().
+ *
+ * When using a cutoff, by default interactions are sharply truncated at the cutoff distance.  Optionally
+ * you can instead use a switching function to make the interaction smoothly go to zero over a finite
+ * distance range.  To enable this, call setUseSwitchingFunction().  You must also call setSwitchingDistance()
+ * to specify the distance at which the interaction should begin to decrease.  The switching distance must be
+ * less than the cutoff distance.
+ */
+
+class OPENMM_EXPORT GayBerneForce : public Force {
+public:
+    /**
+     * This is an enumeration of the different methods that may be used for handling long range nonbonded forces.
+     */
+    enum NonbondedMethod {
+        /**
+         * No cutoff is applied to nonbonded interactions.  The full set of N^2 interactions is computed exactly.
+         * This necessarily means that periodic boundary conditions cannot be used.  This is the default.
+         */
+        NoCutoff = 0,
+        /**
+         * Interactions beyond the cutoff distance are ignored.
+         */
+        CutoffNonPeriodic = 1,
+        /**
+         * Periodic boundary conditions are used, so that each particle interacts only with the nearest periodic copy of
+         * each other particle.  Interactions beyond the cutoff distance are ignored.
+         */
+        CutoffPeriodic = 2
+    };
+    /**
+     * Create a GayBerneForce.
+     */
+    GayBerneForce();
+    /**
+     * Get the number of particles for which force field parameters have been defined.
+     */
+    int getNumParticles() const {
+        return particles.size();
+    }
+    /**
+     * Get the number of special interactions that should be calculated differently from other interactions.
+     */
+    int getNumExceptions() const {
+        return exceptions.size();
+    }
+    /**
+     * Get the method used for handling long range interactions.
+     */
+    NonbondedMethod getNonbondedMethod() const;
+    /**
+     * Set the method used for handling long range interactions.
+     */
+    void setNonbondedMethod(NonbondedMethod method);
+    /**
+     * Get the cutoff distance (in nm) being used for interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     *
+     * @return the cutoff distance, measured in nm
+     */
+    double getCutoffDistance() const;
+    /**
+     * Set the cutoff distance (in nm) being used for interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     *
+     * @param distance    the cutoff distance, measured in nm
+     */
+    void setCutoffDistance(double distance);
+    /**
+     * Get whether a switching function is applied to the interaction.  If the nonbonded method is set
+     * to NoCutoff, this option is ignored.
+     */
+    bool getUseSwitchingFunction() const;
+    /**
+     * Set whether a switching function is applied to the interaction.  If the nonbonded method is set
+     * to NoCutoff, this option is ignored.
+     */
+    void setUseSwitchingFunction(bool use);
+    /**
+     * Get the distance at which the switching function begins to reduce the interaction.  This must be
+     * less than the cutoff distance.
+     */
+    double getSwitchingDistance() const;
+    /**
+     * Set the distance at which the switching function begins to reduce the interaction.  This must be
+     * less than the cutoff distance.
+     */
+    void setSwitchingDistance(double distance);
+    /**
+     * Add the parameters for a particle.  This should be called once for each particle in the System.
+     * When it is called for the i'th time, it specifies the parameters for the i'th particle.
+     *
+     * @param sigma     the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
+     * @param epsilon   the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
+     * @param xparticle the index of the particle whose position defines the ellipsoid's x axis, or -1 if the ellipsoid is a sphere
+     * @param yparticle the index of the particle whose position defines the ellipsoid's y axis, or -1 if the ellipsoid is axially symmetric
+     * @param sx        the diameter of the ellipsoid along its x axis
+     * @param sy        the diameter of the ellipsoid along its y axis
+     * @param sz        the diameter of the ellipsoid along its z axis
+     * @param ex        the factor by which epsilon is scaled along the ellipsoid's x axis
+     * @param ey        the factor by which epsilon is scaled along the ellipsoid's y axis
+     * @param ez        the factor by which epsilon is scaled along the ellipsoid's z axis
+     * @return the index of the particle that was added
+     */
+    int addParticle(double sigma, double epsilon, int xparticle, int yparticle, double sx, double sy, double sz, double ex, double ey, double ez);
+    /**
+     * Get the parameters for a particle.
+     *
+     * @param index          the index of the particle for which to get parameters
+     * @param[out] sigma     the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
+     * @param[out] epsilon   the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
+     * @param[out] xparticle the index of the particle whose position defines the ellipsoid's x axis, or -1 if the ellipsoid is a sphere
+     * @param[out] yparticle the index of the particle whose position defines the ellipsoid's y axis, or -1 if the ellipsoid is axially symmetric
+     * @param[out] sx        the diameter of the ellipsoid along its x axis
+     * @param[out] sy        the diameter of the ellipsoid along its y axis
+     * @param[out] sz        the diameter of the ellipsoid along its z axis
+     * @param[out] ex        the factor by which epsilon is scaled along the ellipsoid's x axis
+     * @param[out] ey        the factor by which epsilon is scaled along the ellipsoid's y axis
+     * @param[out] ez        the factor by which epsilon is scaled along the ellipsoid's z axis
+     */
+    void getParticleParameters(int index, double& sigma, double& epsilon, int& xparticle, int& yparticle, double& sx, double& sy, double& sz, double& ex, double& ey, double& ez) const;
+    /**
+     * Set the parameters for a particle.
+     *
+     * @param index     the index of the particle for which to set parameters
+     * @param sigma     the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
+     * @param epsilon   the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
+     * @param xparticle the index of the particle whose position defines the ellipsoid's x axis, or -1 if the ellipsoid is a sphere
+     * @param yparticle the index of the particle whose position defines the ellipsoid's y axis, or -1 if the ellipsoid is axially symmetric
+     * @param sx        the diameter of the ellipsoid along its x axis
+     * @param sy        the diameter of the ellipsoid along its y axis
+     * @param sz        the diameter of the ellipsoid along its z axis
+     * @param ex        the factor by which epsilon is scaled along the ellipsoid's x axis
+     * @param ey        the factor by which epsilon is scaled along the ellipsoid's y axis
+     * @param ez        the factor by which epsilon is scaled along the ellipsoid's z axis
+     */
+    void setParticleParameters(int index, double sigma, double epsilon, int xparticle, int yparticle, double sx, double sy, double sz, double ex, double ey, double ez);
+    /**
+     * Add an interaction to the list of exceptions that should be calculated differently from other interactions.  If
+     * epsilon is equal to 0, this will cause the interaction to be completely omitted from force and energy calculations.
+     *
+     * @param particle1  the index of the first particle involved in the interaction
+     * @param particle2  the index of the second particle involved in the interaction
+     * @param sigma      the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
+     * @param epsilon    the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
+     * @param replace    determines the behavior if there is already an exception for the same two particles.  If true, the existing one is replaced.  If false,
+     *                   an exception is thrown.
+     * @return the index of the exception that was added
+     */
+    int addException(int particle1, int particle2, double sigma, double epsilon, bool replace = false);
+    /**
+     * Get the force field parameters for an interaction that should be calculated differently from others.
+     *
+     * @param index           the index of the interaction for which to get parameters
+     * @param[out] particle1  the index of the first particle involved in the interaction
+     * @param[out] particle2  the index of the second particle involved in the interaction
+     * @param[out] sigma      the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
+     * @param[out] epsilon    the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
+     */
+    void getExceptionParameters(int index, int& particle1, int& particle2, double& sigma, double& epsilon) const;
+    /**
+     * Set the force field parameters for an interaction that should be calculated differently from others.  If
+     * epsilon is equal to 0, this will cause the interaction to be completely omitted from force and energy calculations.
+     *
+     * @param index      the index of the interaction for which to get parameters
+     * @param particle1  the index of the first particle involved in the interaction
+     * @param particle2  the index of the second particle involved in the interaction
+     * @param sigma      the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
+     * @param epsilon    the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
+     */
+    void setExceptionParameters(int index, int particle1, int particle2, double sigma, double epsilon);
+    /**
+     * Update the particle and exception parameters in a Context to match those stored in this Force object.  This method
+     * provides an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
+     * Simply call setParticleParameters() and setExceptionParameters() to modify this object's parameters, then call
+     * updateParametersInContext() to copy them over to the Context.
+     *
+     * 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.  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);
+    /**
+     * Returns whether or not this force makes use of periodic boundary
+     * conditions.
+     *
+     * @returns true if force uses PBC and false otherwise
+     */
+    bool usesPeriodicBoundaryConditions() const {
+        return nonbondedMethod == GayBerneForce::CutoffPeriodic;
+    }
+protected:
+    ForceImpl* createImpl() const;
+private:
+    class ParticleInfo;
+    class ExceptionInfo;
+    NonbondedMethod nonbondedMethod;
+    double cutoffDistance, switchingDistance;
+    bool useSwitchingFunction;
+    std::vector<ParticleInfo> particles;
+    std::vector<ExceptionInfo> exceptions;
+    std::map<std::pair<int, int>, int> exceptionMap;
+};
+
+/**
+ * This is an internal class used to record information about a particle.
+ * @private
+ */
+class GayBerneForce::ParticleInfo {
+public:
+    int xparticle, yparticle;
+    double sigma, epsilon, sx, sy, sz, ex, ey, ez;
+    ParticleInfo() {
+        xparticle = yparticle = -1;
+        sigma = epsilon = sx = sy = sz = ex = ey = ez = 0.0;
+    }
+    ParticleInfo(double sigma, double epsilon, int xparticle, int yparticle, double sx, double sy, double sz, double ex, double ey, double ez) :
+        sigma(sigma), epsilon(epsilon), xparticle(xparticle), yparticle(yparticle), sx(sx), sy(sy), sz(sz), ex(ex), ey(ey), ez(ez) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about an exception.
+ * @private
+ */
+class GayBerneForce::ExceptionInfo {
+public:
+    int particle1, particle2;
+    double sigma, epsilon;
+    ExceptionInfo() {
+        particle1 = particle2 = -1;
+        sigma = epsilon = 0.0;
+    }
+    ExceptionInfo(int particle1, int particle2, double sigma, double epsilon) :
+        particle1(particle1), particle2(particle2), sigma(sigma), epsilon(epsilon) {
+    }
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_GAYBERNEFORCE_H_*/

openmmapi/include/openmm/internal/GayBerneForceImpl.h

+#ifndef OPENMM_GAYBERNEFORCEIMPL_H_
+#define OPENMM_GAYBERNEFORCEIMPL_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-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 "ForceImpl.h"
+#include "openmm/GayBerneForce.h"
+#include "openmm/Kernel.h"
+#include <utility>
+#include <set>
+#include <string>
+
+namespace OpenMM {
+
+class System;
+
+/**
+ * This is the internal implementation of GayBerneForce.
+ */
+
+class OPENMM_EXPORT GayBerneForceImpl : public ForceImpl {
+public:
+    GayBerneForceImpl(const GayBerneForce& owner);
+    ~GayBerneForceImpl();
+    void initialize(ContextImpl& context);
+    const GayBerneForce& getOwner() const {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context) {
+        // This force field doesn't update the state directly.
+    }
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
+    std::map<std::string, double> getDefaultParameters() {
+        return std::map<std::string, double>(); // This force field doesn't define any parameters.
+    }
+    std::vector<std::string> getKernelNames();
+    void updateParametersInContext(ContextImpl& context);
+private:
+    const GayBerneForce& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_GAYBERNEFORCEIMPL_H_*/

openmmapi/src/GayBerneForce.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-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 "openmm/Force.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/GayBerneForce.h"
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/internal/GayBerneForceImpl.h"
+#include <cmath>
+#include <map>
+#include <sstream>
+#include <utility>
+
+using namespace OpenMM;
+using std::map;
+using std::pair;
+using std::set;
+using std::string;
+using std::stringstream;
+using std::vector;
+
+GayBerneForce::GayBerneForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), switchingDistance(-1.0), useSwitchingFunction(false) {
+}
+
+GayBerneForce::NonbondedMethod GayBerneForce::getNonbondedMethod() const {
+    return nonbondedMethod;
+}
+
+void GayBerneForce::setNonbondedMethod(NonbondedMethod method) {
+    nonbondedMethod = method;
+}
+
+double GayBerneForce::getCutoffDistance() const {
+    return cutoffDistance;
+}
+
+void GayBerneForce::setCutoffDistance(double distance) {
+    cutoffDistance = distance;
+}
+
+bool GayBerneForce::getUseSwitchingFunction() const {
+    return useSwitchingFunction;
+}
+
+void GayBerneForce::setUseSwitchingFunction(bool use) {
+    useSwitchingFunction = use;
+}
+
+double GayBerneForce::getSwitchingDistance() const {
+    return switchingDistance;
+}
+
+void GayBerneForce::setSwitchingDistance(double distance) {
+    switchingDistance = distance;
+}
+
+int GayBerneForce::addParticle(double sigma, double epsilon, int xparticle, int yparticle, double sx, double sy, double sz, double ex, double ey, double ez) {
+    if (yparticle == -1 && (sy != sz || ey != ez))
+        throw OpenMMException("GayBerneForce: yparticle is -1 for a particle that is not axially symmetric");
+    if (xparticle == -1 && (sx != sz || 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, sx, sy, sz, ex, ey, ez));
+    return particles.size()-1;
+}
+
+void GayBerneForce::getParticleParameters(int index, double& sigma, double& epsilon, int& xparticle, int& yparticle, double& sx, double& sy, double& sz, double& ex, double& ey, double& ez) const {
+    ASSERT_VALID_INDEX(index, particles);
+    sigma = particles[index].sigma;
+    epsilon = particles[index].epsilon;
+    xparticle = particles[index].xparticle;
+    yparticle = particles[index].yparticle;
+    sx = particles[index].sx;
+    sy = particles[index].sy;
+    sz = particles[index].sz;
+    ex = particles[index].ex;
+    ey = particles[index].ey;
+    ez = particles[index].ez;
+}
+
+void GayBerneForce::setParticleParameters(int index, double sigma, double epsilon, int xparticle, int yparticle, double sx, double sy, double sz, double ex, double ey, double ez) {
+    ASSERT_VALID_INDEX(index, particles);
+    if (yparticle == -1 && (sy != sz || ey != ez))
+        throw OpenMMException("GayBerneForce: yparticle is -1 for a particle that is not axially symmetric");
+    if (xparticle == -1 && (sx != sz || 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;
+    particles[index].yparticle = yparticle;
+    particles[index].sx = sx;
+    particles[index].sy = sy;
+    particles[index].sz = sz;
+    particles[index].ex = ex;
+    particles[index].ey = ey;
+    particles[index].ez = ez;
+}
+
+int GayBerneForce::addException(int particle1, int particle2, double sigma, double epsilon, bool replace) {
+    map<pair<int, int>, int>::iterator iter = exceptionMap.find(pair<int, int>(particle1, particle2));
+    int newIndex;
+    if (iter == exceptionMap.end())
+        iter = exceptionMap.find(pair<int, int>(particle2, particle1));
+    if (iter != exceptionMap.end()) {
+        if (!replace) {
+            stringstream msg;
+            msg << "GayBerneForce: There is already an exception for particles ";
+            msg << particle1;
+            msg << " and ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        exceptions[iter->second] = ExceptionInfo(particle1, particle2, sigma, epsilon);
+        newIndex = iter->second;
+        exceptionMap.erase(iter->first);
+    }
+    else {
+        exceptions.push_back(ExceptionInfo(particle1, particle2, sigma, epsilon));
+        newIndex = exceptions.size()-1;
+    }
+    exceptionMap[pair<int, int>(particle1, particle2)] = newIndex;
+    return newIndex;
+}
+
+void GayBerneForce::getExceptionParameters(int index, int& particle1, int& particle2, double& sigma, double& epsilon) const {
+    ASSERT_VALID_INDEX(index, exceptions);
+    particle1 = exceptions[index].particle1;
+    particle2 = exceptions[index].particle2;
+    sigma = exceptions[index].sigma;
+    epsilon = exceptions[index].epsilon;
+}
+
+void GayBerneForce::setExceptionParameters(int index, int particle1, int particle2, double sigma, double epsilon) {
+    ASSERT_VALID_INDEX(index, exceptions);
+    exceptions[index].particle1 = particle1;
+    exceptions[index].particle2 = particle2;
+    exceptions[index].sigma = sigma;
+    exceptions[index].epsilon = epsilon;
+}
+
+ForceImpl* GayBerneForce::createImpl() const {
+    return new GayBerneForceImpl(*this);
+}
+
+void GayBerneForce::updateParametersInContext(Context& context) {
+    dynamic_cast<GayBerneForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
+}

openmmapi/src/GayBerneForceImpl.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-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 "openmm/OpenMMException.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/GayBerneForceImpl.h"
+#include "openmm/kernels.h"
+#include <set>
+#include <sstream>
+
+using namespace OpenMM;
+using namespace std;
+
+GayBerneForceImpl::GayBerneForceImpl(const GayBerneForce& owner) : owner(owner) {
+}
+
+GayBerneForceImpl::~GayBerneForceImpl() {
+}
+
+void GayBerneForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcGayBerneForceKernel::Name(), context);
+
+    // Check for errors in the specification of exceptions.
+
+    const System& system = context.getSystem();
+    if (owner.getNumParticles() != system.getNumParticles())
+        throw OpenMMException("GayBerneForce must have exactly as many particles as the System it belongs to.");
+    if (owner.getUseSwitchingFunction()) {
+        if (owner.getSwitchingDistance() < 0 || owner.getSwitchingDistance() >= owner.getCutoffDistance())
+            throw OpenMMException("GayBerneForce: Switching distance must satisfy 0 <= r_switch < r_cutoff");
+    }
+    for (int i = 0; i < owner.getNumParticles(); i++) {
+        int xparticle, yparticle;
+        double sigma, epsilon, rx, ry, rz, ex, ey, ez;
+        owner.getParticleParameters(i, sigma, epsilon, xparticle, yparticle, rx, ry, rz, ex, ey, ez);
+        if (xparticle < -1 || xparticle >= owner.getNumParticles()) {
+            stringstream msg;
+            msg << "GayBerneForce: Illegal particle index for xparticle: ";
+            msg << xparticle;
+            throw OpenMMException(msg.str());
+        }
+        if (yparticle < -1 || yparticle >= owner.getNumParticles()) {
+            stringstream msg;
+            msg << "GayBerneForce: Illegal particle index for a yparticle: ";
+            msg << yparticle;
+            throw OpenMMException(msg.str());
+        }
+    }
+    vector<set<int> > exceptions(owner.getNumParticles());
+    for (int i = 0; i < owner.getNumExceptions(); i++) {
+        int particle1, particle2;
+        double sigma, epsilon;
+        owner.getExceptionParameters(i, particle1, particle2, sigma, epsilon);
+        if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
+            stringstream msg;
+            msg << "GayBerneForce: Illegal particle index for an exception: ";
+            msg << particle1;
+            throw OpenMMException(msg.str());
+        }
+        if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
+            stringstream msg;
+            msg << "GayBerneForce: Illegal particle index for an exception: ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        if (exceptions[particle1].count(particle2) > 0 || exceptions[particle2].count(particle1) > 0) {
+            stringstream msg;
+            msg << "GayBerneForce: Multiple exceptions are specified for particles ";
+            msg << particle1;
+            msg << " and ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        exceptions[particle1].insert(particle2);
+        exceptions[particle2].insert(particle1);
+    }
+    if (owner.getNonbondedMethod() == GayBerneForce::CutoffPeriodic) {
+        Vec3 boxVectors[3];
+        system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+        double cutoff = owner.getCutoffDistance();
+        if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
+            throw OpenMMException("GayBerneForce: The cutoff distance cannot be greater than half the periodic box size.");
+    }
+    kernel.getAs<CalcGayBerneForceKernel>().initialize(context.getSystem(), owner);
+}
+
+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);
+    return 0.0;
+}
+
+std::vector<std::string> GayBerneForceImpl::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(CalcGayBerneForceKernel::Name());
+    return names;
+}
+
+void GayBerneForceImpl::updateParametersInContext(ContextImpl& context) {
+    kernel.getAs<CalcGayBerneForceKernel>().copyParametersToContext(context, owner);
+}

platforms/cpu/include/CpuGayBerneForce.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 OPENMM_CPU_GAYBERNEFORCE_H__
+#define OPENMM_CPU_GAYBERNEFORCE_H__
+
+#include "openmm/GayBerneForce.h"
+#include "openmm/internal/ThreadPool.h"
+#include "CpuNeighborList.h"
+#include "CpuPlatform.h"
+#include "RealVec.h"
+#include <set>
+#include <utility>
+
+namespace OpenMM {
+
+class CpuGayBerneForce {
+public:
+    struct Matrix;
+    class ComputeTask;
+
+    /**
+     * Constructor.
+     */
+    CpuGayBerneForce(const GayBerneForce& force);
+
+    /**
+     * Compute the interaction.
+     *
+     * @param positions     the positions of the atoms
+     * @param forces        forces will be added to this vector
+     * @param threadForce   individual threads can add their forces to this vector
+     * @param boxVectors    the periodic box vectors
+     * @param data          the platform data for the current context
+     * @return the energy of the interaction
+     */
+    RealOpenMM calculateForce(const std::vector<RealVec>& positions, std::vector<RealVec>& forces, std::vector<AlignedArray<float> >& threadForce, RealVec* boxVectors, CpuPlatform::PlatformData& data);
+
+    /**
+     * This routine contains the code executed by each thread.
+     */
+    void threadComputeForce(ThreadPool& threads, int threadIndex, CpuNeighborList* neighborList);
+    
+    /**
+     * Get the exclusions being used by the force.
+     */
+    const std::vector<std::set<int> >& getExclusions() const;
+
+private:
+    struct ParticleInfo;
+    struct ExceptionInfo;
+    std::vector<ParticleInfo> particles;
+    std::vector<ExceptionInfo> exceptions;
+    std::set<std::pair<int, int> > exclusions;
+    std::vector<std::set<int> > particleExclusions;
+    GayBerneForce::NonbondedMethod nonbondedMethod;
+    RealOpenMM cutoffDistance, switchingDistance;
+    bool useSwitchingFunction;
+    std::vector<RealOpenMM> s;
+    std::vector<Matrix> A, B, G;
+    std::vector<double> threadEnergy;
+    std::vector<std::vector<RealVec> > threadTorque;
+    // The following variables are used to make information accessible to the individual threads.
+    RealVec const* positions;
+    std::vector<AlignedArray<float> >* threadForce;
+    RealVec* boxVectors;
+    void* atomicCounter;
+
+    void computeEllipsoidFrames(const std::vector<RealVec>& positions);
+    
+    void applyTorques(const std::vector<RealVec>& positions, std::vector<RealVec>& forces);
+
+    RealOpenMM computeOneInteraction(int particle1, int particle2, RealOpenMM sigma, RealOpenMM epsilon, const RealVec* positions,
+            float* forces, std::vector<RealVec>& torques, const RealVec* boxVectors);
+};
+
+struct CpuGayBerneForce::ParticleInfo {
+    int xparticle, yparticle;
+    RealOpenMM sigmaOver2, sqrtEpsilon, rx, ry, rz, ex, ey, ez;
+    bool isPointParticle;
+};
+
+struct CpuGayBerneForce::ExceptionInfo {
+    int particle1, particle2;
+    RealOpenMM sigma, epsilon;
+};
+
+struct CpuGayBerneForce::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] = 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;
+    }
+};
+
+static RealVec operator*(const RealVec& r, CpuGayBerneForce::Matrix& m) {
+    return RealVec(m.v[0][0]*r[0] + m.v[1][0]*r[1] + m.v[2][0]*r[2],
+                   m.v[0][1]*r[0] + m.v[1][1]*r[1] + m.v[2][1]*r[2],
+                   m.v[0][2]*r[0] + m.v[1][2]*r[1] + m.v[2][2]*r[2]);
+}
+
+} // namespace OpenMM
+
+#endif // OPENMM_CPU_GAYBERNEFORCE_H__

platforms/cpu/include/CpuKernels.h

@@ -36,6 +36,7 @@
 #include "CpuCustomGBForce.h"
 #include "CpuCustomManyParticleForce.h"
 #include "CpuCustomNonbondedForce.h"
+#include "CpuGayBerneForce.h"
 #include "CpuGBSAOBCForce.h"
 #include "CpuLangevinDynamics.h"
 #include "CpuNeighborList.h"
@@ -446,6 +447,42 @@ private:
     NonbondedMethod nonbondedMethod;
 };
 
+/**
+ * This kernel is invoked by GayBerneForce to calculate the forces acting on the system.
+ */
+class CpuCalcGayBerneForceKernel : public CalcGayBerneForceKernel {
+public:
+    CpuCalcGayBerneForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) : CalcGayBerneForceKernel(name, platform),
+            data(data), ixn(NULL) {
+    }
+    ~CpuCalcGayBerneForceKernel();
+    /**
+     * 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:
+    CpuPlatform::PlatformData& data;
+    CpuGayBerneForce* ixn;
+};
+
 /**
  * This kernel is invoked by LangevinIntegrator to take one time step.
  */

platforms/cpu/include/CpuPlatform.h

@@ -82,7 +82,7 @@ class CpuPlatform::PlatformData {
 public:
     PlatformData(int numParticles, int numThreads);
     ~PlatformData();
-    void requestNeighborList(double cutoffDistance, double padding, bool useExclusions, std::vector<std::set<int> >& exclusionList);
+    void requestNeighborList(double cutoffDistance, double padding, bool useExclusions, const std::vector<std::set<int> >& exclusionList);
     AlignedArray<float> posq;
     std::vector<AlignedArray<float> > threadForce;
     ThreadPool threads;

platforms/cpu/src/CpuGayBerneForce.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.                                     *
+ * -------------------------------------------------------------------------- */
+
+#ifdef _MSC_VER
+    // Prevent Windows from defining macros that interfere with other code.
+    #define NOMINMAX
+#endif
+#include "CpuGayBerneForce.h"
+#include "ReferenceForce.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/GayBerneForce.h"
+#include "openmm/internal/gmx_atomic.h"
+#include <algorithm>
+#include <cmath>
+
+using namespace OpenMM;
+using namespace std;
+
+class CpuGayBerneForce::ComputeTask : public ThreadPool::Task {
+public:
+    ComputeTask(CpuGayBerneForce& owner, CpuNeighborList* neighborList) : owner(owner), neighborList(neighborList) {
+    }
+    void execute(ThreadPool& threads, int threadIndex) {
+        owner.threadComputeForce(threads, threadIndex, neighborList);
+    }
+    CpuGayBerneForce& owner;
+    CpuNeighborList* neighborList;
+};
+
+CpuGayBerneForce::CpuGayBerneForce(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];
+        double sigma, epsilon, sx, sy, sz, ex, ey, ez;
+        force.getParticleParameters(i, sigma, epsilon, p.xparticle, p.yparticle, sx, sy, sz, ex, ey, ez);
+        p.sigmaOver2 = 0.5*sigma;
+        p.sqrtEpsilon = sqrt(epsilon);
+        p.rx = 0.5*sx;
+        p.ry = 0.5*sy;
+        p.rz = 0.5*sz;
+        p.ex = ex;
+        p.ey = ey;
+        p.ez = ez;
+        p.isPointParticle = (sx == sigma && sy == sigma && sz == sigma && ex == 1.0 && ey == 1.0 && ez == 1.0);
+    }
+    int numExceptions = force.getNumExceptions();
+    exceptions.resize(numExceptions);
+    particleExclusions.resize(numParticles);
+    for (int i = 0; i < numExceptions; i++) {
+        ExceptionInfo& e = exceptions[i];
+        double sigma, epsilon;
+        force.getExceptionParameters(i, e.particle1, e.particle2, sigma, epsilon);
+        e.sigma = sigma;
+        e.epsilon = epsilon;
+        exclusions.insert(make_pair(min(e.particle1, e.particle2), max(e.particle1, e.particle2)));
+        particleExclusions[e.particle1].insert(e.particle2);
+        particleExclusions[e.particle2].insert(e.particle1);
+    }
+    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)*sqrtf(p.rx*p.ry);
+    }
+}
+
+const vector<set<int> >& CpuGayBerneForce::getExclusions() const {
+    return particleExclusions;
+}
+
+RealOpenMM CpuGayBerneForce::calculateForce(const vector<RealVec>& positions, std::vector<RealVec>& forces, std::vector<AlignedArray<float> >& threadForce, RealVec* boxVectors, CpuPlatform::PlatformData& data) {
+    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);
+
+    // Record the parameters for the threads.
+    
+    ThreadPool& threads = data.threads;
+    int numThreads = threads.getNumThreads();
+    this->positions = &positions[0];
+    this->threadForce = &threadForce;
+    this->boxVectors = boxVectors;
+    threadEnergy.resize(numThreads);
+    threadTorque.resize(numThreads);
+    gmx_atomic_t counter;
+    gmx_atomic_set(&counter, 0);
+    this->atomicCounter = &counter;
+    
+    // Signal the threads to compute the pairwise interactions.
+    
+    ComputeTask task(*this, data.neighborList);
+    threads.execute(task);
+    threads.waitForThreads();
+    
+    // Signal the threads to compute exceptions.
+    
+    gmx_atomic_set(&counter, 0);
+    threads.resumeThreads();
+    threads.waitForThreads();
+    
+    // Combine the energies from all the threads.
+    
+    double energy = 0;
+    for (int i = 0; i < numThreads; i++)
+        energy += threadEnergy[i];
+    
+    // Apply torques.
+    
+    applyTorques(positions, forces);
+    return energy;
+}
+
+void CpuGayBerneForce::threadComputeForce(ThreadPool& threads, int threadIndex, CpuNeighborList* neighborList) {
+    int numParticles = particles.size();
+    int numThreads = threads.getNumThreads();
+    threadEnergy[threadIndex] = 0;
+    float* forces = &(*threadForce)[threadIndex][0];
+    vector<RealVec>& torques = threadTorque[threadIndex];
+    torques.resize(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        torques[i] = RealVec();
+    double energy = 0.0;
+
+    // Compute this thread's subset of interactions.
+    
+    if (neighborList == NULL) {
+        while (true) {
+            int i = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (i >= numParticles)
+                break;
+            if (particles[i].sqrtEpsilon == 0.0f)
+                continue;
+            for (int j = 0; j < i; j++) {
+                if (particles[j].sqrtEpsilon == 0.0f)
+                    continue;
+                if (particleExclusions[i].find(j) != particleExclusions[i].end())
+                    continue; // This interaction will be handled by an exception.
+                RealOpenMM sigma = particles[i].sigmaOver2+particles[j].sigmaOver2;
+                RealOpenMM epsilon = particles[i].sqrtEpsilon*particles[j].sqrtEpsilon;
+                energy += computeOneInteraction(i, j, sigma, epsilon, positions, forces, torques, boxVectors);
+            }
+        }
+    }
+    else {
+        while (true) {
+            int blockIndex = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (blockIndex >= neighborList->getNumBlocks())
+                break;
+            const int blockSize = neighborList->getBlockSize();
+            const int* blockAtom = &neighborList->getSortedAtoms()[blockSize*blockIndex];
+            const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
+            const vector<char>& exclusions = neighborList->getBlockExclusions(blockIndex);
+            for (int i = 0; i < (int) neighbors.size(); i++) {
+                int first = neighbors[i];
+                if (particles[first].sqrtEpsilon == 0.0f)
+                    continue;
+                for (int k = 0; k < blockSize; k++) {
+                    if ((exclusions[i] & (1<<k)) == 0) {
+                        int second = blockAtom[k];
+                        if (particles[second].sqrtEpsilon == 0.0f)
+                            continue;
+                        RealOpenMM sigma = particles[first].sigmaOver2+particles[second].sigmaOver2;
+                        RealOpenMM epsilon = particles[first].sqrtEpsilon*particles[second].sqrtEpsilon;
+                        energy += computeOneInteraction(first, second, sigma, epsilon, positions, forces, torques, boxVectors);
+                    }
+                }
+            }
+        }
+    }
+
+    // Compute exceptions.
+
+    threads.syncThreads();
+    int numExceptions = exceptions.size();
+    const int groupSize = max(1, numExceptions/(10*numThreads));
+    while (true) {
+        int start = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), groupSize);
+        if (start >= numExceptions)
+            break;
+        int end = min(start+groupSize, numExceptions);
+        for (int i = start; i < end; i++) {
+            ExceptionInfo& e = exceptions[i];
+            energy += computeOneInteraction(e.particle1, e.particle2, e.sigma, e.epsilon, positions, forces, torques, boxVectors);
+        }
+    }
+    threadEnergy[threadIndex] = energy;
+}
+
+void CpuGayBerneForce::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(1/sqrt(p.ex), 1/sqrt(p.ey), 1/sqrt(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];
+                }
+            }
+    }
+}
+
+void CpuGayBerneForce::applyTorques(const vector<RealVec>& positions, vector<RealVec>& forces) {
+    int numParticles = particles.size();
+    int numThreads = threadTorque.size();
+    for (int particle = 0; particle < numParticles; particle++) {
+        ParticleInfo& p = particles[particle];
+        RealVec pos = positions[particle];
+        if (p.xparticle != -1) {
+            // Add up the torques from the individual threads.
+            
+            RealVec torque;
+            for (int i = 0; i < numThreads; i++)
+                torque += threadTorque[i][particle];
+            
+            // Apply a force to the x particle.
+            
+            RealVec dx = positions[p.xparticle]-pos;
+            double dx2 = dx.dot(dx);
+            RealVec f = torque.cross(dx)/dx2;
+            forces[p.xparticle] += f;
+            forces[particle] -= f;
+            if (p.yparticle != -1) {
+                // Apply a force to the y particle.  This is based on the component of the torque
+                // that was not already applied to the x particle.
+                
+                RealVec dy = positions[p.yparticle]-pos;
+                double dy2 = dy.dot(dy);
+                RealVec torque2 = dx*(torque.dot(dx)/dx2);
+                f = torque2.cross(dy)/dy2;
+                forces[p.yparticle] += f;
+                forces[particle] -= f;
+            }
+        }
+    }
+}
+
+RealOpenMM CpuGayBerneForce::computeOneInteraction(int particle1, int particle2, RealOpenMM sigma, RealOpenMM epsilon, const RealVec* positions,
+        float* forces, vector<RealVec>& torques, 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 r = deltaR[ReferenceForce::RIndex];
+    if (nonbondedMethod != GayBerneForce::NoCutoff && r >= cutoffDistance)
+        return 0;
+    RealOpenMM rInv = 1/r;
+    RealVec dr(deltaR[ReferenceForce::XIndex], deltaR[ReferenceForce::YIndex], deltaR[ReferenceForce::ZIndex]);
+    RealVec drUnit = dr*rInv;
+    
+    // Compute the switching function.
+
+    RealOpenMM switchValue = 1, switchDeriv = 0;
+    if (useSwitchingFunction && r > switchingDistance) {
+        RealOpenMM t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
+        switchValue = 1+t*t*t*(-10+t*(15-t*6));
+        switchDeriv = t*t*(-30+t*(60-t*30))/(cutoffDistance-switchingDistance);
+    }
+
+    // Interactions between two point particles can be computed more easily.
+    
+    if (particles[particle1].isPointParticle && particles[particle2].isPointParticle) {
+        RealOpenMM sig = sigma*rInv;
+        RealOpenMM sig2 = sig*sig;
+        RealOpenMM sig6 = sig2*sig2*sig2;
+        RealOpenMM energy = 4*epsilon*(sig6-1)*sig6;
+        RealVec force = drUnit*(switchValue*4*epsilon*(12*sig6 - 6)*sig6*rInv - energy*switchDeriv);
+        forces[4*particle1] += force[0];
+        forces[4*particle1+1] += force[1];
+        forces[4*particle1+2] += force[2];
+        forces[4*particle2] -= force[0];
+        forces[4*particle2+1] -= force[1];
+        forces[4*particle2+2] -= force[2];
+        return energy*switchValue;
+    }
+
+    // Compute vectors and matrices we'll be needing.
+
+    Matrix B12 = B[particle1]+B[particle2];
+    Matrix G12 = G[particle1]+G[particle2];
+    Matrix B12inv = B12.inverse();
+    Matrix G12inv = G12.inverse();
+    RealOpenMM detG12 = G12.determinant();
+
+    // Estimate the distance between the ellipsoids and compute the first terms needed for the energy.
+
+    RealOpenMM sigma12 = 1/SQRT(0.5*drUnit.dot(G12inv*drUnit));
+    RealOpenMM h12 = r - sigma12;
+    RealOpenMM rho = sigma/(h12+sigma);
+    RealOpenMM rho2 = rho*rho;
+    RealOpenMM rho6 = rho2*rho2*rho2;
+    RealOpenMM u = 4*epsilon*(rho6*rho6-rho6);
+    RealOpenMM eta = SQRT(2*s[particle1]*s[particle2]/detG12);
+    RealOpenMM chi = 2*drUnit.dot(B12inv*drUnit);
+    chi *= chi;
+    RealOpenMM energy = u*eta*chi;
+    
+    // Compute the terms needed for the force.
+
+    RealVec kappa = G12inv*dr;
+    RealVec iota = B12inv*dr;
+    RealOpenMM rInv2 = rInv*rInv;
+    RealOpenMM dUSLJdr = 24*epsilon*(2*rho6-1)*rho6*rho/sigma;
+    RealOpenMM temp = 0.5*sigma12*sigma12*sigma12*rInv2;
+    RealVec dudr = (drUnit + (kappa-drUnit*kappa.dot(drUnit))*temp)*dUSLJdr;
+    RealVec dchidr = (iota-drUnit*iota.dot(drUnit))*(-8*rInv2*SQRT(chi));
+    RealVec force = (dchidr*u + dudr*chi)*(eta*switchValue) - drUnit*(energy*switchDeriv);
+    forces[4*particle1] += force[0];
+    forces[4*particle1+1] += force[1];
+    forces[4*particle1+2] += force[2];
+    forces[4*particle2] -= force[0];
+    forces[4*particle2+1] -= force[1];
+    forces[4*particle2+2] -= force[2];
+
+    // Compute the terms needed for the torque.
+
+    for (int j = 0; j < 2; j++) {
+        int particle = (j == 0 ? particle1 : particle2);
+        ParticleInfo& p = particles[particle];
+        if (p.isPointParticle)
+            continue;
+        RealVec dudq = (kappa*G[particle]).cross(kappa*(temp*dUSLJdr));
+        RealVec dchidq = (iota*B[particle]).cross(iota)*(-4*rInv2);
+        RealOpenMM (&g12)[3][3] = G12.v;
+        RealOpenMM (&a)[3][3] = A[particle].v;
+        RealVec scale = RealVec(p.rx*p.rx, p.ry*p.ry, p.rz*p.rz)*(-0.5*eta/detG12);
+        Matrix D;
+        RealOpenMM (&d)[3][3] = D.v;
+        d[0][0] = scale[0]*(2*a[0][0]*(g12[1][1]*g12[2][2] - g12[1][2]*g12[2][1]) +
+                              a[0][2]*(g12[1][2]*g12[0][1] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])) +
+                              a[0][1]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])));
+        d[0][1] = scale[0]*(  a[0][0]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])) +
+                            2*a[0][1]*(g12[0][0]*g12[2][2] - g12[2][0]*g12[0][2]) +
+                              a[0][2]*(g12[1][0]*g12[0][2] + g12[2][0]*g12[0][1] - g12[0][0]*(g12[1][2] + g12[2][1])));
+        d[0][2] = scale[0]*(  a[0][0]*(g12[0][1]*g12[1][2] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])) +
+                              a[0][1]*(g12[1][0]*g12[0][2] + g12[2][0]*g12[0][1] - g12[0][0]*(g12[1][2] + g12[2][1])) +
+                            2*a[0][2]*(g12[1][1]*g12[0][0] - g12[1][0]*g12[0][1]));
+        d[1][0] = scale[1]*(2*a[1][0]*(g12[1][1]*g12[2][2] - g12[1][2]*g12[2][1]) +
+                              a[1][1]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])) +
+                              a[1][2]*(g12[1][2]*g12[0][1] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])));
+        d[1][1] = scale[1]*(  a[1][0]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])) +
+                            2*a[1][1]*(g12[2][2]*g12[0][0] - g12[2][0]*g12[0][2]) +
+                              a[1][2]*(g12[1][0]*g12[0][2] + g12[0][1]*g12[2][0] - g12[0][0]*(g12[1][2] + g12[2][1])));
+        d[1][2] = scale[1]*(  a[1][0]*(g12[0][1]*g12[1][2] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])) +
+                              a[1][1]*(g12[1][0]*g12[0][2] + g12[0][1]*g12[2][0] - g12[0][0]*(g12[1][2] + g12[2][1])) +
+                            2*a[1][2]*(g12[1][1]*g12[0][0] - g12[1][0]*g12[0][1]));
+        d[2][0] = scale[2]*(2*a[2][0]*(g12[1][1]*g12[2][2] - g12[2][1]*g12[1][2]) +
+                              a[2][1]*(g12[0][2]*g12[2][1] + g12[1][2]*g12[2][0] - g12[2][2]*(g12[0][1] + g12[1][0])) +
+                              a[2][2]*(g12[0][1]*g12[1][2] + g12[2][1]*g12[1][0] - g12[1][1]*(g12[0][2] + g12[2][0])));
+        d[2][1] = scale[2]*(  a[2][0]*(g12[0][2]*g12[2][1] + g12[1][2]*g12[2][0] - g12[2][2]*(g12[0][1] + g12[1][0])) +
+                            2*a[2][1]*(g12[0][0]*g12[2][2] - g12[0][2]*g12[2][0]) +
+                              a[2][2]*(g12[1][0]*g12[0][2] + g12[0][1]*g12[2][0] - g12[0][0]*(g12[1][2] + g12[2][1])));
+        d[2][2] = scale[2]*(  a[2][0]*(g12[0][1]*g12[1][2] + g12[2][1]*g12[1][0] - g12[1][1]*(g12[0][2] + g12[2][0])) +
+                              a[2][1]*(g12[1][0]*g12[0][2] + g12[2][0]*g12[0][1] - g12[0][0]*(g12[1][2] + g12[2][1])) +
+                            2*a[2][2]*(g12[1][1]*g12[0][0] - g12[1][0]*g12[0][1]));
+        RealVec detadq;
+        for (int i = 0; i < 3; i++)
+            detadq += RealVec(a[i][0], a[i][1], a[i][2]).cross(RealVec(d[i][0], d[i][1], d[i][2]));
+        RealVec torque = (dchidq*(u*eta) + detadq*(u*chi) + dudq*(eta*chi))*switchValue;
+        torques[particle] -= torque;
+    }
+    return switchValue*energy;
+}

platforms/cpu/src/CpuKernelFactory.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) 2013 Stanford University and the Authors.           *
+ * Portions copyright (c) 2013-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -57,6 +57,8 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform&
         return new CpuCalcGBSAOBCForceKernel(name, platform, data);
     if (name == CalcCustomGBForceKernel::Name())
         return new CpuCalcCustomGBForceKernel(name, platform, data);
+    if (name == CalcGayBerneForceKernel::Name())
+        return new CpuCalcGayBerneForceKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())
         return new CpuIntegrateLangevinStepKernel(name, platform, data);
     throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str());

platforms/cpu/src/CpuKernels.cpp

@@ -1236,6 +1236,30 @@ void CpuCalcCustomManyParticleForceKernel::copyParametersToContext(ContextImpl&
     }
 }
 
+CpuCalcGayBerneForceKernel::~CpuCalcGayBerneForceKernel() {
+    if (ixn != NULL)
+        delete ixn;
+}
+
+void CpuCalcGayBerneForceKernel::initialize(const System& system, const GayBerneForce& force) {
+    ixn = new CpuGayBerneForce(force);
+    data.isPeriodic = (force.getNonbondedMethod() == GayBerneForce::CutoffPeriodic);
+    if (force.getNonbondedMethod() != GayBerneForce::NoCutoff) {
+        double cutoff = force.getCutoffDistance();
+        data.requestNeighborList(cutoff, 0.1*cutoff, true, ixn->getExclusions());
+    }
+}
+
+double CpuCalcGayBerneForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    return ixn->calculateForce(extractPositions(context), extractForces(context), data.threadForce, extractBoxVectors(context), data);
+}
+
+void CpuCalcGayBerneForceKernel::copyParametersToContext(ContextImpl& context, const GayBerneForce& force) {
+    delete ixn;
+    ixn = NULL;
+    ixn = new CpuGayBerneForce(force);
+}
+
 CpuIntegrateLangevinStepKernel::~CpuIntegrateLangevinStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/cpu/src/CpuPlatform.cpp

@@ -71,6 +71,7 @@ CpuPlatform::CpuPlatform() {
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
+    registerKernelFactory(CalcGayBerneForceKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     platformProperties.push_back(CpuThreads());
     int threads = getNumProcessors();
@@ -155,7 +156,7 @@ CpuPlatform::PlatformData::~PlatformData() {
 
 bool isVec8Supported();
 
-void CpuPlatform::PlatformData::requestNeighborList(double cutoffDistance, double padding, bool useExclusions, vector<set<int> >& exclusionList) {
+void CpuPlatform::PlatformData::requestNeighborList(double cutoffDistance, double padding, bool useExclusions, const vector<set<int> >& exclusionList) {
     if (neighborList == NULL)
         neighborList = new CpuNeighborList(isVec8Supported() ? 8 : 4);
     if (cutoffDistance > cutoff)

platforms/cpu/tests/TestCpuGayBerneForce.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 "CpuTests.h"
+#include "TestGayBerneForce.h"
+
+void runPlatformTests() {
+}

platforms/cuda/include/CudaKernels.h

@@ -1104,6 +1104,74 @@ private:
     CUevent event;
 };
 
+/**
+ * This kernel is invoked by GayBerneForce to calculate the forces acting on the system.
+ */
+class CudaCalcGayBerneForceKernel : public CalcGayBerneForceKernel {
+public:
+    CudaCalcGayBerneForceKernel(std::string name, const Platform& platform, CudaContext& cu) : CalcGayBerneForceKernel(name, platform), cu(cu),
+            hasInitializedKernels(false), sortedParticles(NULL), axisParticleIndices(NULL), sigParams(NULL), epsParams(NULL), scale(NULL), exceptionParticles(NULL),
+            exceptionParams(NULL), aMatrix(NULL),
+            bMatrix(NULL), gMatrix(NULL), exclusions(NULL), exclusionStartIndex(NULL), blockCenter(NULL), blockBoundingBox(NULL), neighbors(NULL),
+            neighborIndex(NULL), neighborBlockCount(NULL), sortedPos(NULL), torque(NULL) {
+    }
+    ~CudaCalcGayBerneForceKernel();
+    /**
+     * 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:
+    class ReorderListener;
+    void sortAtoms();
+    CudaContext& cu;
+    bool hasInitializedKernels;
+    int numRealParticles, numExceptions, maxNeighborBlocks;
+    GayBerneForce::NonbondedMethod nonbondedMethod;
+    CudaArray* sortedParticles;
+    CudaArray* axisParticleIndices;
+    CudaArray* sigParams;
+    CudaArray* epsParams;
+    CudaArray* scale;
+    CudaArray* exceptionParticles;
+    CudaArray* exceptionParams;
+    CudaArray* aMatrix;
+    CudaArray* bMatrix;
+    CudaArray* gMatrix;
+    CudaArray* exclusions;
+    CudaArray* exclusionStartIndex;
+    CudaArray* blockCenter;
+    CudaArray* blockBoundingBox;
+    CudaArray* neighbors;
+    CudaArray* neighborIndex;
+    CudaArray* neighborBlockCount;
+    CudaArray* sortedPos;
+    CudaArray* torque;
+    std::vector<bool> isRealParticle;
+    std::vector<std::pair<int, int> > exceptionAtoms;
+    std::vector<std::pair<int, int> > excludedPairs;
+    std::vector<void*> framesArgs, blockBoundsArgs, neighborsArgs, forceArgs, torqueArgs;
+    CUfunction framesKernel, blockBoundsKernel, neighborsKernel, forceKernel, torqueKernel;
+    CUevent event;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/cuda/src/CudaKernelFactory.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-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -110,6 +110,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcCustomCompoundBondForceKernel(name, platform, cu, context.getSystem());
     if (name == CalcCustomManyParticleForceKernel::Name())
         return new CudaCalcCustomManyParticleForceKernel(name, platform, cu, context.getSystem());
+    if (name == CalcGayBerneForceKernel::Name())
+        return new CudaCalcGayBerneForceKernel(name, platform, cu);
     if (name == IntegrateVerletStepKernel::Name())
         return new CudaIntegrateVerletStepKernel(name, platform, cu);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/cuda/src/CudaKernels.cpp

@@ -5988,6 +5988,432 @@ void CudaCalcCustomManyParticleForceKernel::copyParametersToContext(ContextImpl&
     cu.invalidateMolecules();
 }
 
+class CudaGayBerneForceInfo : public CudaForceInfo {
+public:
+    CudaGayBerneForceInfo(const GayBerneForce& force) : force(force) {
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        int xparticle1, yparticle1;
+        double sigma1, epsilon1, sx1, sy1, sz1, ex1, ey1, ez1;
+        int xparticle2, yparticle2;
+        double sigma2, epsilon2, sx2, sy2, sz2, ex2, ey2, ez2;
+        force.getParticleParameters(particle1, sigma1, epsilon1, xparticle1, yparticle1, sx1, sy1, sz1, ex1, ey1, ez1);
+        force.getParticleParameters(particle2, sigma2, epsilon2, xparticle2, yparticle2, sx2, sy2, sz2, ex2, ey2, ez2);
+        return (sigma1 == sigma2 && epsilon1 == epsilon2 && sx1 == sx2 && sy1 == sy2 && sz1 == sz2 && ex1 == ex2 && ey1 == ey2 && ez1 == ez2);
+    }
+    int getNumParticleGroups() {
+        return force.getNumExceptions()+force.getNumParticles();
+    }
+    void getParticlesInGroup(int index, vector<int>& particles) {
+        if (index < force.getNumExceptions()) {
+            int particle1, particle2;
+            double sigma, epsilon;
+            force.getExceptionParameters(index, particle1, particle2, sigma, epsilon);
+            particles.resize(2);
+            particles[0] = particle1;
+            particles[1] = particle2;
+        }
+        else {
+            int particle = index-force.getNumExceptions();
+            int xparticle, yparticle;
+            double sigma, epsilon, sx, sy, sz, ex, ey, ez;
+            force.getParticleParameters(particle, sigma, epsilon, xparticle, yparticle, sx, sy, sz, ex, ey, ez);
+            particles.clear();
+            particles.push_back(particle);
+            if (xparticle > -1)
+                particles.push_back(xparticle);
+            if (yparticle > -1)
+                particles.push_back(yparticle);
+        }
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        if (group1 < force.getNumExceptions() && group2 < force.getNumExceptions()) {
+            int particle1, particle2;
+            double sigma1, sigma2, epsilon1, epsilon2;
+            force.getExceptionParameters(group1, particle1, particle2, sigma1, epsilon1);
+            force.getExceptionParameters(group2, particle1, particle2, sigma2, epsilon2);
+            return (sigma1 == sigma2 && epsilon1 == epsilon2);
+        }
+        return true;
+    }
+private:
+    const GayBerneForce& force;
+};
+
+class CudaCalcGayBerneForceKernel::ReorderListener : public CudaContext::ReorderListener {
+public:
+    ReorderListener(CudaCalcGayBerneForceKernel& owner) : owner(owner) {
+    }
+    void execute() {
+        owner.sortAtoms();
+    }
+private:
+    CudaCalcGayBerneForceKernel& owner;
+};
+
+CudaCalcGayBerneForceKernel::~CudaCalcGayBerneForceKernel() {
+    if (sortedParticles != NULL)
+        delete sortedParticles;
+    if (axisParticleIndices != NULL)
+        delete axisParticleIndices;
+    if (sigParams != NULL)
+        delete sigParams;
+    if (epsParams != NULL)
+        delete epsParams;
+    if (scale != NULL)
+        delete scale;
+    if (exceptionParticles != NULL)
+        delete exceptionParticles;
+    if (exceptionParams != NULL)
+        delete exceptionParams;
+    if (aMatrix != NULL)
+        delete aMatrix;
+    if (bMatrix != NULL)
+        delete bMatrix;
+    if (gMatrix != NULL)
+        delete gMatrix;
+    if (exclusions != NULL)
+        delete exclusions;
+    if (exclusionStartIndex != NULL)
+        delete exclusionStartIndex;
+    if (blockCenter != NULL)
+        delete blockCenter;
+    if (blockBoundingBox != NULL)
+        delete blockBoundingBox;
+    if (neighbors != NULL)
+        delete neighbors;
+    if (neighborIndex != NULL)
+        delete neighborIndex;
+    if (neighborBlockCount != NULL)
+        delete neighborBlockCount;
+    if (sortedPos != NULL)
+        delete sortedPos;
+    if (torque != NULL)
+        delete torque;
+}
+
+void CudaCalcGayBerneForceKernel::initialize(const System& system, const GayBerneForce& force) {
+    // Initialize interactions.
+
+    int numParticles = force.getNumParticles();
+    sigParams = CudaArray::create<float4>(cu, cu.getPaddedNumAtoms(), "sigParams");
+    epsParams = CudaArray::create<float2>(cu, cu.getPaddedNumAtoms(), "epsParams");
+    scale = CudaArray::create<float4>(cu, cu.getPaddedNumAtoms(), "scale");
+    axisParticleIndices = CudaArray::create<int2>(cu, cu.getPaddedNumAtoms(), "axisParticleIndices");
+    sortedParticles = CudaArray::create<int>(cu, cu.getPaddedNumAtoms(), "sortedParticles");
+    aMatrix = CudaArray::create<float>(cu, 9*cu.getPaddedNumAtoms(), "aMatrix");
+    bMatrix = CudaArray::create<float>(cu, 9*cu.getPaddedNumAtoms(), "bMatrix");
+    gMatrix = CudaArray::create<float>(cu, 9*cu.getPaddedNumAtoms(), "gMatrix");
+    vector<float4> sigParamsVector(cu.getPaddedNumAtoms(), make_float4(0, 0, 0, 0));
+    vector<float2> epsParamsVector(cu.getPaddedNumAtoms(), make_float2(0, 0));
+    vector<float4> scaleVector(cu.getPaddedNumAtoms(), make_float4(0, 0, 0, 0));
+    vector<int2> axisParticleVector(cu.getPaddedNumAtoms(), make_int2(0, 0));
+    isRealParticle.resize(cu.getPaddedNumAtoms());
+    for (int i = 0; i < numParticles; i++) {
+        int xparticle, yparticle;
+        double sigma, epsilon, sx, sy, sz, ex, ey, ez;
+        force.getParticleParameters(i, sigma, epsilon, xparticle, yparticle, sx, sy, sz, ex, ey, ez);
+        axisParticleVector[i] = make_int2(xparticle, yparticle);
+        sigParamsVector[i] = make_float4((float) (0.5*sigma), (float) (0.25*sx*sx), (float) (0.25*sy*sy), (float) (0.25*sz*sz));
+        epsParamsVector[i] = make_float2((float) sqrt(epsilon), (float) (0.125*(sx*sy + sz*sz)*sqrt(sx*sy)));
+        scaleVector[i] = make_float4((float) (1/sqrt(ex)), (float) (1/sqrt(ey)), (float) (1/sqrt(ez)), 0);
+        isRealParticle[i] = (epsilon != 0.0);
+    }
+    sigParams->upload(sigParamsVector);
+    epsParams->upload(epsParamsVector);
+    scale->upload(scaleVector);
+    axisParticleIndices->upload(axisParticleVector);
+    
+    // Record exceptions and exclusions.
+
+    vector<float2> exceptionParamsVec;
+    for (int i = 0; i < force.getNumExceptions(); i++) {
+        int particle1, particle2;
+        double sigma, epsilon;
+        force.getExceptionParameters(i, particle1, particle2, sigma, epsilon);
+        if (epsilon != 0.0) {
+            exceptionParamsVec.push_back(make_float2((float) sigma, (float) epsilon));
+            exceptionAtoms.push_back(make_pair(particle1, particle2));
+            isRealParticle[particle1] = true;
+            isRealParticle[particle2] = true;
+        }
+        if (isRealParticle[particle1] && isRealParticle[particle2])
+            excludedPairs.push_back(pair<int, int>(particle1, particle2));
+    }
+    numRealParticles = 0;
+    for (int i = 0; i < isRealParticle.size(); i++)
+        if (isRealParticle[i])
+            numRealParticles++;
+    numExceptions = exceptionParamsVec.size();
+    exclusions = CudaArray::create<int>(cu, max(1, (int) excludedPairs.size()), "exclusions");
+    exclusionStartIndex = CudaArray::create<int>(cu, numRealParticles+1, "exclusionStartIndex");
+    exceptionParticles = CudaArray::create<int4>(cu, max(1, numExceptions), "exceptionParticles");
+    exceptionParams = CudaArray::create<float2>(cu, max(1, numExceptions), "exceptionParams");
+    if (numExceptions > 0)
+        exceptionParams->upload(exceptionParamsVec);
+    
+    // Create data structures used for the neighbor list.
+
+    int numAtomBlocks = (numRealParticles+31)/32;
+    int elementSize = (cu.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
+    blockCenter = new CudaArray(cu, numAtomBlocks, 4*elementSize, "blockCenter");
+    blockBoundingBox = new CudaArray(cu, numAtomBlocks, 4*elementSize, "blockBoundingBox");
+    sortedPos = new CudaArray(cu, numRealParticles, 4*elementSize, "sortedPos");
+    maxNeighborBlocks = numRealParticles*2;
+    neighbors = CudaArray::create<int>(cu, maxNeighborBlocks*32, "neighbors");
+    neighborIndex = CudaArray::create<int>(cu, maxNeighborBlocks, "neighborIndex");
+    neighborBlockCount = CudaArray::create<int>(cu, 1, "neighborBlockCount");
+    if (force.getNonbondedMethod() != GayBerneForce::NoCutoff)
+        CHECK_RESULT(cuEventCreate(&event, CU_EVENT_DISABLE_TIMING), "Error creating event for CustomManyParticleForce");
+
+    // Create array for accumulating torques.
+    
+    torque = CudaArray::create<long long>(cu, 3*cu.getPaddedNumAtoms(), "torque");
+    cu.addAutoclearBuffer(*torque);
+
+    // Create the kernels.
+    
+    nonbondedMethod = force.getNonbondedMethod();
+    bool useCutoff = (nonbondedMethod != GayBerneForce::NoCutoff);
+    bool usePeriodic = (nonbondedMethod == GayBerneForce::CutoffPeriodic);
+    map<string, string> defines;
+    defines["USE_SWITCH"] = (useCutoff && force.getUseSwitchingFunction() ? "1" : "0");
+    double cutoff = force.getCutoffDistance();
+    defines["CUTOFF_SQUARED"] = cu.doubleToString(cutoff*cutoff);
+    if (useCutoff) {
+        defines["USE_CUTOFF"] = 1;
+        if (usePeriodic)
+            defines["USE_PERIODIC"] = "1";
+        
+        // Compute the switching coefficients.
+        
+        if (force.getUseSwitchingFunction()) {
+            defines["SWITCH_CUTOFF"] = cu.doubleToString(force.getSwitchingDistance());
+            defines["SWITCH_C3"] = cu.doubleToString(10/pow(force.getSwitchingDistance()-cutoff, 3.0));
+            defines["SWITCH_C4"] = cu.doubleToString(15/pow(force.getSwitchingDistance()-cutoff, 4.0));
+            defines["SWITCH_C5"] = cu.doubleToString(6/pow(force.getSwitchingDistance()-cutoff, 5.0));
+        }
+    }
+    defines["PADDED_NUM_ATOMS"] = cu.intToString(cu.getPaddedNumAtoms());
+    CUmodule module = cu.createModule(CudaKernelSources::vectorOps+CudaKernelSources::gayBerne, defines);
+    framesKernel = cu.getKernel(module, "computeEllipsoidFrames");
+    blockBoundsKernel = cu.getKernel(module, "findBlockBounds");
+    neighborsKernel = cu.getKernel(module, "findNeighbors");
+    forceKernel = cu.getKernel(module, "computeForce");
+    torqueKernel = cu.getKernel(module, "applyTorques");
+    cu.addForce(new CudaGayBerneForceInfo(force));
+    cu.addReorderListener(new ReorderListener(*this));
+}
+
+double CudaCalcGayBerneForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    if (!hasInitializedKernels) {
+        hasInitializedKernels = true;
+        sortAtoms();
+        framesArgs.push_back(&numRealParticles);
+        framesArgs.push_back(&cu.getPosq().getDevicePointer());
+        framesArgs.push_back(&axisParticleIndices->getDevicePointer());
+        framesArgs.push_back(&sigParams->getDevicePointer());
+        framesArgs.push_back(&scale->getDevicePointer());
+        framesArgs.push_back(&aMatrix->getDevicePointer());
+        framesArgs.push_back(&bMatrix->getDevicePointer());
+        framesArgs.push_back(&gMatrix->getDevicePointer());
+        framesArgs.push_back(&sortedParticles->getDevicePointer());
+        blockBoundsArgs.push_back(&numRealParticles);
+        blockBoundsArgs.push_back(cu.getPeriodicBoxSizePointer());
+        blockBoundsArgs.push_back(cu.getInvPeriodicBoxSizePointer());
+        blockBoundsArgs.push_back(cu.getPeriodicBoxVecXPointer());
+        blockBoundsArgs.push_back(cu.getPeriodicBoxVecYPointer());
+        blockBoundsArgs.push_back(cu.getPeriodicBoxVecZPointer());
+        blockBoundsArgs.push_back(&sortedParticles->getDevicePointer());
+        blockBoundsArgs.push_back(&cu.getPosq().getDevicePointer());
+        blockBoundsArgs.push_back(&sortedPos->getDevicePointer());
+        blockBoundsArgs.push_back(&blockCenter->getDevicePointer());
+        blockBoundsArgs.push_back(&blockBoundingBox->getDevicePointer());
+        blockBoundsArgs.push_back(&neighborBlockCount->getDevicePointer());
+        neighborsArgs.push_back(&numRealParticles);
+        neighborsArgs.push_back(&maxNeighborBlocks);
+        neighborsArgs.push_back(cu.getPeriodicBoxSizePointer());
+        neighborsArgs.push_back(cu.getInvPeriodicBoxSizePointer());
+        neighborsArgs.push_back(cu.getPeriodicBoxVecXPointer());
+        neighborsArgs.push_back(cu.getPeriodicBoxVecYPointer());
+        neighborsArgs.push_back(cu.getPeriodicBoxVecZPointer());
+        neighborsArgs.push_back(&sortedPos->getDevicePointer());
+        neighborsArgs.push_back(&blockCenter->getDevicePointer());
+        neighborsArgs.push_back(&blockBoundingBox->getDevicePointer());
+        neighborsArgs.push_back(&neighbors->getDevicePointer());
+        neighborsArgs.push_back(&neighborIndex->getDevicePointer());
+        neighborsArgs.push_back(&neighborBlockCount->getDevicePointer());
+        neighborsArgs.push_back(&exclusions->getDevicePointer());
+        neighborsArgs.push_back(&exclusionStartIndex->getDevicePointer());
+        forceArgs.push_back(&cu.getForce().getDevicePointer());
+        forceArgs.push_back(&torque->getDevicePointer());
+        forceArgs.push_back(&numRealParticles);
+        forceArgs.push_back(&numExceptions);
+        forceArgs.push_back(&cu.getEnergyBuffer().getDevicePointer());
+        forceArgs.push_back(&sortedPos->getDevicePointer());
+        forceArgs.push_back(&sigParams->getDevicePointer());
+        forceArgs.push_back(&epsParams->getDevicePointer());
+        forceArgs.push_back(&sortedParticles->getDevicePointer());
+        forceArgs.push_back(&aMatrix->getDevicePointer());
+        forceArgs.push_back(&bMatrix->getDevicePointer());
+        forceArgs.push_back(&gMatrix->getDevicePointer());
+        forceArgs.push_back(&exclusions->getDevicePointer());
+        forceArgs.push_back(&exclusionStartIndex->getDevicePointer());
+        forceArgs.push_back(&exceptionParticles->getDevicePointer());
+        forceArgs.push_back(&exceptionParams->getDevicePointer());
+        if (nonbondedMethod != GayBerneForce::NoCutoff) {
+            forceArgs.push_back(&maxNeighborBlocks);
+            forceArgs.push_back(&neighbors->getDevicePointer());
+            forceArgs.push_back(&neighborIndex->getDevicePointer());
+            forceArgs.push_back(&neighborBlockCount->getDevicePointer());
+            forceArgs.push_back(cu.getPeriodicBoxSizePointer());
+            forceArgs.push_back(cu.getInvPeriodicBoxSizePointer());
+            forceArgs.push_back(cu.getPeriodicBoxVecXPointer());
+            forceArgs.push_back(cu.getPeriodicBoxVecYPointer());
+            forceArgs.push_back(cu.getPeriodicBoxVecZPointer());
+        }
+        torqueArgs.push_back(&cu.getForce().getDevicePointer());
+        torqueArgs.push_back(&torque->getDevicePointer());
+        torqueArgs.push_back(&numRealParticles);
+        torqueArgs.push_back(&cu.getPosq().getDevicePointer());
+        torqueArgs.push_back(&axisParticleIndices->getDevicePointer());
+        torqueArgs.push_back(&sortedParticles->getDevicePointer());
+    }
+    cu.executeKernel(framesKernel, &framesArgs[0], numRealParticles);
+    cu.executeKernel(blockBoundsKernel, &blockBoundsArgs[0], (numRealParticles+31)/32);
+    if (nonbondedMethod == GayBerneForce::NoCutoff) {
+        cu.executeKernel(forceKernel, &forceArgs[0], cu.getNonbondedUtilities().getNumForceThreadBlocks()*cu.getNonbondedUtilities().getForceThreadBlockSize());
+    }
+    else {
+        while (true) {
+            cu.executeKernel(neighborsKernel, &neighborsArgs[0], numRealParticles);
+            int* count = (int*) cu.getPinnedBuffer();
+            neighborBlockCount->download(count, false);
+            CHECK_RESULT(cuEventRecord(event, 0), "Error recording event for GayBerneForce");
+            cu.executeKernel(forceKernel, &forceArgs[0], cu.getNonbondedUtilities().getNumForceThreadBlocks()*cu.getNonbondedUtilities().getForceThreadBlockSize());
+            CHECK_RESULT(cuEventSynchronize(event), "Error synchronizing on event for GayBerneForce");
+            if (*count <= maxNeighborBlocks)
+                break;
+            
+            // There wasn't enough room for the neighbor list, so we need to recreate it.
+
+            delete neighbors;
+            neighbors = NULL;
+            delete neighborIndex;
+            neighborIndex = NULL;
+            maxNeighborBlocks = (int) ceil((*count)*1.1);
+            neighbors = CudaArray::create<int>(cu, maxNeighborBlocks*32, "neighbors");
+            neighborIndex = CudaArray::create<int>(cu, maxNeighborBlocks, "neighborIndex");
+            neighborsArgs[10] = &neighbors->getDevicePointer();
+            neighborsArgs[11] = &neighborIndex->getDevicePointer();
+            forceArgs[17] = &neighbors->getDevicePointer();
+            forceArgs[18] = &neighborIndex->getDevicePointer();
+        }
+    }
+    cu.executeKernel(torqueKernel, &torqueArgs[0], numRealParticles);
+    return 0.0;
+}
+
+void CudaCalcGayBerneForceKernel::copyParametersToContext(ContextImpl& context, const GayBerneForce& force) {
+    // Make sure the new parameters are acceptable.
+    
+    if (force.getNumParticles() != cu.getNumAtoms())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+    vector<int> exceptions;
+    for (int i = 0; i < force.getNumExceptions(); i++) {
+        int particle1, particle2;
+        double sigma, epsilon;
+        force.getExceptionParameters(i, particle1, particle2, sigma, epsilon);
+        if (exceptionAtoms.size() > exceptions.size() && make_pair(particle1, particle2) == exceptionAtoms[exceptions.size()])
+            exceptions.push_back(i);
+        else if (epsilon != 0.0)
+            throw OpenMMException("updateParametersInContext: The set of non-excluded exceptions has changed");
+    }
+    int numExceptions = exceptionAtoms.size();
+    
+    // Record the per-particle parameters.
+    
+    vector<float4> sigParamsVector(cu.getPaddedNumAtoms(), make_float4(0, 0, 0, 0));
+    vector<float2> epsParamsVector(cu.getPaddedNumAtoms(), make_float2(0, 0));
+    vector<float4> scaleVector(cu.getPaddedNumAtoms(), make_float4(0, 0, 0, 0));
+    for (int i = 0; i < force.getNumParticles(); i++) {
+        int xparticle, yparticle;
+        double sigma, epsilon, sx, sy, sz, ex, ey, ez;
+        force.getParticleParameters(i, sigma, epsilon, xparticle, yparticle, sx, sy, sz, ex, ey, ez);
+        sigParamsVector[i] = make_float4((float) (0.5*sigma), (float) (0.25*sx*sx), (float) (0.25*sy*sy), (float) (0.25*sz*sz));
+        epsParamsVector[i] = make_float2((float) sqrt(epsilon), (float) (0.125*(sx*sy + sz*sz)*sqrt(sx*sy)));
+        scaleVector[i] = make_float4((float) (1/sqrt(ex)), (float) (1/sqrt(ey)), (float) (1/sqrt(ez)), 0);
+        if (epsilon != 0.0 && !isRealParticle[i])
+            throw OpenMMException("updateParametersInContext: The set of ignored particles (ones with epsilon=0) has changed");
+    }
+    sigParams->upload(sigParamsVector);
+    epsParams->upload(epsParamsVector);
+    scale->upload(scaleVector);
+    
+    // Record the exceptions.
+    
+    if (numExceptions > 0) {
+        vector<float2> exceptionParamsVec(numExceptions);
+        for (int i = 0; i < numExceptions; i++) {
+            int atom1, atom2;
+            double sigma, epsilon;
+            force.getExceptionParameters(exceptions[i], atom1, atom2, sigma, epsilon);
+            exceptionParamsVec[i] = make_float2((float) sigma, (float) epsilon);
+        }
+        exceptionParams->upload(exceptionParamsVec);
+    }
+    cu.invalidateMolecules();
+    sortAtoms();
+}
+
+void CudaCalcGayBerneForceKernel::sortAtoms() {
+    // Sort the list of atoms by type to avoid thread divergence.  This is executed every time
+    // the atoms are reordered.
+    
+    int nextIndex = 0;
+    vector<int> particles(cu.getPaddedNumAtoms(), 0);
+    const vector<int>& order = cu.getAtomIndex();
+    vector<int> inverseOrder(order.size(), -1);
+    for (int i = 0; i < cu.getNumAtoms(); i++) {
+        int atom = order[i];
+        if (isRealParticle[atom]) {
+            inverseOrder[atom] = nextIndex;
+            particles[nextIndex++] = atom;
+        }
+    }
+    sortedParticles->upload(particles);
+    
+    // Update the list of exception particles.
+    
+    int numExceptions = exceptionAtoms.size();
+    if (numExceptions > 0) {
+        vector<int4> exceptionParticlesVec(numExceptions);
+        for (int i = 0; i < numExceptions; i++)
+            exceptionParticlesVec[i] = make_int4(exceptionAtoms[i].first, exceptionAtoms[i].second, inverseOrder[exceptionAtoms[i].first], inverseOrder[exceptionAtoms[i].second]);
+        exceptionParticles->upload(exceptionParticlesVec);
+    }
+    
+    // Rebuild the list of exclusions.
+    
+    vector<vector<int> > excludedAtoms(numRealParticles);
+    for (int i = 0; i < excludedPairs.size(); i++) {
+        int first = inverseOrder[min(excludedPairs[i].first, excludedPairs[i].second)];
+        int second = inverseOrder[max(excludedPairs[i].first, excludedPairs[i].second)];
+        excludedAtoms[first].push_back(second);
+    }
+    int index = 0;
+    vector<int> exclusionVec(exclusions->getSize());
+    vector<int> startIndexVec(exclusionStartIndex->getSize());
+    for (int i = 0; i < numRealParticles; i++) {
+        startIndexVec[i] = index;
+        for (int j = 0; j < excludedAtoms[i].size(); j++)
+            exclusionVec[index++] = excludedAtoms[i][j];
+    }
+    startIndexVec[numRealParticles] = index;
+    exclusions->upload(exclusionVec);
+    exclusionStartIndex->upload(startIndexVec);
+}
+
 CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
 }
 

platforms/cuda/src/CudaPlatform.cpp

@@ -92,6 +92,7 @@ CudaPlatform::CudaPlatform() {
     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);