Continuing reference implementation of Gay-Berne force

openmmapi/include/openmm/GayBerneForce.h

@@ -47,7 +47,9 @@ namespace OpenMM {
  * sigma and epsilon parameters, each particle has three radii r_x, r_y, and r_z that give the size of the
  * ellipsoid along each axis.  It also has three scale factors e_x, e_y, and e_z that scale the strength
  * of the interaction along each axis.  You can think of this force as a Lennard-Jones interaction computed
- * based on the distance between the nearest points on two ellipsoids.  If two particles each have all their
+ * based on the distance between the nearest points on two ellipsoids.  The interpretation of the scale
+ * factors is slightly complicated, but roughly speaking, the strength of the interaction along the ellipsoid's
+ * x axis is divided by e_x^2, and likewise for the other axes.  If two particles each have all their
  * radii set to sigma/2 and all their scale factors set 1, the interaction simplifies to a standard Lennard-Jones
  * force between point particles.
  *

platforms/reference/include/ReferenceGayBerneForce.h

@@ -77,7 +77,6 @@ private:
 struct ReferenceGayBerneForce::ParticleInfo {
     int xparticle, yparticle;
     RealOpenMM sigma, epsilon, rx, ry, rz, ex, ey, ez;
-    bool radiiAreZero, scalesAreZero;
 };
 
 struct ReferenceGayBerneForce::ExceptionInfo {

platforms/reference/src/SimTKReference/ReferenceGayBerneForce.cpp

@@ -45,8 +45,6 @@ ReferenceGayBerneForce::ReferenceGayBerneForce(const GayBerneForce& force) {
     for (int i = 0; i < numParticles; i++) {
         ParticleInfo& p = particles[i];
         force.getParticleParameters(i, p.sigma, p.epsilon, p.xparticle, p.yparticle, p.rx, p.ry, p.rz, p.ex, p.ey, p.ez);
-        p.radiiAreZero = (p.rx == 0 && p.ry == 0 && p.rz == 0);
-        p.scalesAreZero = (p.ex == 0 && p.ey == 0 && p.ez == 0);
     }
     int numExceptions = force.getNumExceptions();
     exceptions.resize(numExceptions);
@@ -90,14 +88,19 @@ RealOpenMM ReferenceGayBerneForce::calculateForce(const vector<RealVec>& positio
 
     RealOpenMM energy = 0;
     int numParticles = particles.size();
-    for (int i = 1; i < numParticles; i++)
+    for (int i = 1; i < numParticles; i++) {
+        if (particles[i].epsilon == 0.0)
+            continue;
         for (int j = 0; j < i; j++) {
+            if (particles[j].epsilon == 0.0)
+                continue;
             if (exclusions.find(make_pair(j, i)) != exclusions.end())
                 continue; // This interaction will be handled by an exception.
             RealOpenMM sigma = 0.5*(particles[i].sigma+particles[j].sigma);
             RealOpenMM epsilon = SQRT(particles[i].epsilon*particles[j].epsilon);
             energy += computeOneInteraction(i, j, sigma, epsilon, positions, forces, boxVectors);
         }
+    }
 
     // Compute exceptions.
 
@@ -173,38 +176,42 @@ RealOpenMM ReferenceGayBerneForce::computeOneInteraction(int particle1, int part
         ReferenceForce::getDeltaRPeriodic(positions[particle2], positions[particle1], boxVectors, deltaR);
     else
         ReferenceForce::getDeltaR(positions[particle2], positions[particle1], deltaR);
-    RealOpenMM dist = deltaR[ReferenceForce::RIndex];
-    if (nonbondedMethod != GayBerneForce::NoCutoff && dist >= cutoffDistance)
+    RealOpenMM r = deltaR[ReferenceForce::RIndex];
+    if (nonbondedMethod != GayBerneForce::NoCutoff && r >= cutoffDistance)
         return 0;
 
     // Compute vectors and matrices we'll be needing.
 
+    RealOpenMM rInv = 1/r;
     RealVec dr(deltaR[ReferenceForce::XIndex], deltaR[ReferenceForce::YIndex], deltaR[ReferenceForce::ZIndex]);
-    RealVec drUnit = dr/dist;
+    RealVec drUnit = dr*rInv;
     Matrix B12 = B[particle1]+B[particle2];
     Matrix G12 = G[particle1]+G[particle2];
     Matrix B12inv = B12.inverse();
     Matrix G12inv = G12.inverse();
 
-    // Estimate the distance between the ellipsoids and compute the first term in the energy.
+    // Estimate the distance between the ellipsoids and compute the first terms needed for the energy.
 
     ParticleInfo& p1 = particles[particle1];
     ParticleInfo& p2 = particles[particle2];
-    RealOpenMM h12 = dist;
-    if (!p1.radiiAreZero || !p2.radiiAreZero)
-        h12 -= 1/SQRT(0.5*drUnit.dot(G12inv*drUnit));
+    RealOpenMM 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);
-
-    // Compute the second term in the energy.
-
     RealOpenMM eta = SQRT(2*s[particle1]*s[particle2]/G12.determinant());
-
-    // Compute the third term in the energy.
-
     RealOpenMM chi = 2*drUnit.dot(B12inv*drUnit);
     chi *= chi;
+
+    // Compute the terms needed for the force.
+
+    RealVec kappa = G12inv*dr;
+    RealVec iota = B12inv*dr;
+    RealVec dudr = (drUnit + (kappa-drUnit*kappa.dot(drUnit))*(0.5*sigma12*sigma12*sigma12*rInv*rInv))*(24*epsilon*(2*rho6-1)*rho6*rho/sigma);
+    RealVec dchidr = (iota-drUnit*iota.dot(drUnit))*(-8*rInv*rInv*SQRT(chi));
+    RealVec force = (dchidr*u + dudr*chi)*eta;
+    forces[particle1] += force;
+    forces[particle2] -= force;
     return u*eta*chi;
 }

tests/TestGayBerneForce.h

@@ -77,6 +77,62 @@ void testPointParticles() {
     State state1 = context.getState(State::Forces | State::Energy, false, 1);
     State state2 = context.getState(State::Forces | State::Energy, false, 2);
     ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
+}
+
+void testEnergyScales() {
+    // Create two Lennard-Jones particles for which the energy scale factors vary.
+
+    const double sigma = 0.5;
+    const double epsilon = 1.5;
+    System system;
+    for (int i = 0; i < 6; i++)
+        system.addParticle(1.0);
+    GayBerneForce* gb = new GayBerneForce();
+    system.addForce(gb);
+    gb->addParticle(sigma, epsilon, 1, 2, sigma/2, sigma/2, sigma/2, 1.1, 1.5, 1.8);
+    gb->addParticle(1, 0, -1, -1, 1, 1, 1, 1, 1, 1);
+    gb->addParticle(1, 0, -1, -1, 1, 1, 1, 1, 1, 1);
+    gb->addParticle(sigma, epsilon, 4, 5, sigma/2, sigma/2, sigma/2, 1.2, 1.6, 1.7);
+    gb->addParticle(1, 0, -1, -1, 1, 1, 1, 1, 1, 1);
+    gb->addParticle(1, 0, -1, -1, 1, 1, 1, 1, 1, 1);
+    vector<Vec3> positions(6);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    positions[2] = Vec3(0, 1, 0);
+    positions[3] = Vec3(1, 0, 0);
+    positions[4] = Vec3(2, 0, 0);
+    positions[5] = Vec3(1, 1, 0);
+    VerletIntegrator integ(0.001);
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+
+    // Depending on their relative orientations, the interaction should be equivalent
+    // to LJ with different values of epsilon.
+
+    double expectedEnergy = 4*epsilon*(pow(sigma, 12.0)-pow(sigma, 6.0));
+    double expectedForce = 4*epsilon*(12*pow(sigma, 12.0)-6*pow(sigma, 6.0));
+    double expectedScale = pow(2.0/(1.1*1.1+1.2*1.2), 2.0);
+    State state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(expectedEnergy*expectedScale, state.getPotentialEnergy(), 1e-5);
+    ASSERT_EQUAL_VEC(Vec3(expectedForce*expectedScale, 0, 0), state.getForces()[3], 1e-5);
+    positions[3] = Vec3(0, 1, 0);
+    positions[4] = Vec3(1, 1, 0);
+    positions[5] = Vec3(0, 2, 0);
+    context.setPositions(positions);
+    expectedScale = pow(2.0/(1.5*1.5+1.6*1.6), 2.0);
+    state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(expectedEnergy*expectedScale, state.getPotentialEnergy(), 1e-5);
+    ASSERT_EQUAL_VEC(Vec3(0, expectedForce*expectedScale, 0), state.getForces()[3], 1e-5);
+    positions[3] = Vec3(0, 1, 0);
+    positions[4] = Vec3(1, 1, 0);
+    positions[5] = Vec3(0, 1, 1);
+    context.setPositions(positions);
+    expectedScale = pow(2.0/(1.5*1.5+1.7*1.7), 2.0);
+    state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(expectedEnergy*expectedScale, state.getPotentialEnergy(), 1e-5);
+    ASSERT_EQUAL_VEC(Vec3(0, expectedForce*expectedScale, 0), state.getForces()[3], 1e-5);
 }
 
 void runPlatformTests();
@@ -85,6 +141,7 @@ int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
         testPointParticles();
+        testEnergyScales();
         runPlatformTests();
     }
     catch(const exception& e) {