AmoebaVdwForce works with triclinic boxes

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaVdwForce.cpp

@@ -39,6 +39,7 @@
 #include "openmm/System.h"
 #include "openmm/AmoebaVdwForce.h"
 #include "openmm/LangevinIntegrator.h"
+#include "sfmt/SFMT.h"
 #include <iostream>
 #include <vector>
 #include <stdlib.h>
@@ -50,6 +51,7 @@
 
 
 using namespace OpenMM;
+using namespace std;
 
 extern "C" void registerAmoebaCudaKernelFactories();
 
@@ -1981,6 +1983,62 @@ void testVdwWater( int includeVdwDispersionCorrection, FILE* log ) {
     }
 }
 
+void testTriclinic() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    Vec3 a(3.1, 0, 0);
+    Vec3 b(0.4, 3.5, 0);
+    Vec3 c(-0.1, -0.5, 4.0);
+    system.setDefaultPeriodicBoxVectors(a, b, c);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    AmoebaVdwForce* vdw = new AmoebaVdwForce();
+    vdw->setUseDispersionCorrection(false);
+    vdw->addParticle(0, 0.5, 1.0, 0.0);
+    vdw->addParticle(1, 0.5, 1.0, 0.0);
+    vdw->setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
+    const double cutoff = 1.5;
+    vdw->setCutoff(cutoff);
+    system.addForce(vdw);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+    vector<Vec3> positions(2);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int iteration = 0; iteration < 50; iteration++) {
+        // Generate random positions for the two particles.
+
+        positions[0] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        positions[1] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        context.setPositions(positions);
+
+        // Loop over all possible periodic copies and find the nearest one.
+
+        Vec3 delta;
+        double distance2 = 100.0;
+        for (int i = -1; i < 2; i++)
+            for (int j = -1; j < 2; j++)
+                for (int k = -1; k < 2; k++) {
+                    Vec3 d = positions[1]-positions[0]+a*i+b*j+c*k;
+                    if (d.dot(d) < distance2) {
+                        delta = d;
+                        distance2 = d.dot(d);
+                    }
+                }
+        double distance = sqrt(distance2);
+
+        // See if the energy is correct.
+
+        State state = context.getState(State::Energy);
+        if (distance >= cutoff) {
+            ASSERT_EQUAL(0.0, state.getPotentialEnergy());
+        }
+        else if (distance < 0.9*cutoff) {
+            const double energy = pow(1.07/(distance+0.07), 7.0)*(1.12/(pow(distance, 7.0)+0.12)-2);
+            ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+        }
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         std::cout << "TestCudaAmoebaVdwForce running test..." << std::endl;
@@ -2029,6 +2087,10 @@ int main(int argc, char* argv[]) {
  
         includeVdwDispersionCorrection     = 1;
         testVdwWater( includeVdwDispersionCorrection, log );
+        
+        // test triclinic boxes
+
+        testTriclinic();
 
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp

@@ -976,14 +976,14 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
         computeNeighborListVoxelHash( *neighborList, numParticles, posData, allExclusions, extractBoxVectors(context), usePBC, cutoff, 0.0);
         if( usePBC ){
             vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::CutoffPeriodic);
-            RealVec& box = extractBoxSize(context);
+            RealVec* boxVectors = extractBoxVectors(context);
             double minAllowedSize = 1.999999*cutoff;
-            if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize){
+            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 cutoff.");
             }
-            vdwForce.setPeriodicBox(box);
+            vdwForce.setPeriodicBox(boxVectors);
             energy  = vdwForce.calculateForceAndEnergy( numParticles, posData, indexIVs, sigmas, epsilons, reductions, *neighborList, forceData);
-            energy += dispersionCoefficient/(box[0]*box[1]*box[2]);
+            energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
         } else {
             vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::CutoffNonPeriodic);
         }

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp

@@ -24,18 +24,18 @@
 
 #include "AmoebaReferenceForce.h"
 #include "AmoebaReferenceVdwForce.h"
+#include "ReferenceForce.h"
 #include <algorithm>
 #include <cctype>
 
 using std::vector;
-using OpenMM::RealVec;
+using namespace OpenMM;
 
 AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( ) : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9) {
 
     setTaperCoefficients( _cutoff );
     setSigmaCombiningRule( "ARITHMETIC" );
     setEpsilonCombiningRule( "GEOMETRIC" );
-    _periodicBoxDimensions = RealVec( 0.0, 0.0, 0.0 );
 }
 
 
@@ -44,7 +44,6 @@ AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( const std::string& sigmaCombin
     setTaperCoefficients( _cutoff );
     setSigmaCombiningRule( sigmaCombiningRule );
     setEpsilonCombiningRule( epsilonCombiningRule );
-    _periodicBoxDimensions = RealVec( 0.0, 0.0, 0.0 );
 }
 
 AmoebaReferenceVdwForce::NonbondedMethod AmoebaReferenceVdwForce::getNonbondedMethod( void ) const {
@@ -77,12 +76,10 @@ double AmoebaReferenceVdwForce::getCutoff( void ) const {
     return _cutoff;
 }
 
-void AmoebaReferenceVdwForce::setPeriodicBox( const RealVec& box ){
-    _periodicBoxDimensions = box;
-}
-
-RealVec AmoebaReferenceVdwForce::getPeriodicBox( void ) const {
-    return _periodicBoxDimensions;
+void AmoebaReferenceVdwForce::setPeriodicBox(OpenMM::RealVec* vectors) {
+    _periodicBoxVectors[0] = vectors[0];
+    _periodicBoxVectors[1] = vectors[1];
+    _periodicBoxVectors[2] = vectors[2];
 }
 
 void AmoebaReferenceVdwForce::setSigmaCombiningRule( const std::string& sigmaCombiningRule ){
@@ -181,7 +178,7 @@ void AmoebaReferenceVdwForce::addReducedForce( unsigned int particleI, unsigned
     forces[particleIV][2] += sign*force[2]*(one - reduction);
 }
 
-RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma, RealOpenMM combindedEpsilon,
+RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combinedSigma, RealOpenMM combinedEpsilon,
                                                       const Vec3& particleIPosition,
                                                       const Vec3& particleJPosition,
                                                       Vec3& force ) const {
@@ -197,36 +194,34 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma,
 
    // ---------------------------------------------------------------------------------------
 
-    RealOpenMM xr           = particleIPosition[0] - particleJPosition[0];
-    RealOpenMM yr           = particleIPosition[1] - particleJPosition[1];
-    RealOpenMM zr           = particleIPosition[2] - particleJPosition[2];
+    // get deltaR, R2, and R between 2 atoms
 
-    if( _nonbondedMethod == CutoffPeriodic ){
-               xr          -= static_cast<RealOpenMM>((floor(xr/_periodicBoxDimensions[0]+0.5)*_periodicBoxDimensions[0]));
-               yr          -= static_cast<RealOpenMM>((floor(yr/_periodicBoxDimensions[1]+0.5)*_periodicBoxDimensions[1]));
-               zr          -= static_cast<RealOpenMM>((floor(zr/_periodicBoxDimensions[2]+0.5)*_periodicBoxDimensions[2]));
-    }
+    RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+    if (_nonbondedMethod == CutoffPeriodic)
+        ReferenceForce::getDeltaRPeriodic(particleJPosition, particleIPosition, _periodicBoxVectors, deltaR);
+    else
+        ReferenceForce::getDeltaR(particleJPosition, particleIPosition, deltaR);
 
-    RealOpenMM r_ij_2       = xr*xr + yr*yr + zr*zr;
-    RealOpenMM r_ij         = SQRT(r_ij_2);
-    RealOpenMM sigma_7      = combindedSigma*combindedSigma*combindedSigma;
-               sigma_7      = sigma_7*sigma_7*combindedSigma;
+    RealOpenMM r_ij_2       = deltaR[ReferenceForce::R2Index];
+    RealOpenMM r_ij         = deltaR[ReferenceForce::RIndex];
+    RealOpenMM sigma_7      = combinedSigma*combinedSigma*combinedSigma;
+               sigma_7      = sigma_7*sigma_7*combinedSigma;
 
     RealOpenMM r_ij_6       = r_ij_2*r_ij_2*r_ij_2;
     RealOpenMM r_ij_7       = r_ij_6*r_ij;
 
     RealOpenMM rho          = r_ij_7 + ghal*sigma_7;
 
-    RealOpenMM tau          = (dhal + one)/(r_ij + dhal*combindedSigma);
+    RealOpenMM tau          = (dhal + one)/(r_ij + dhal*combinedSigma);
     RealOpenMM tau_7        = tau*tau*tau;
                tau_7        = tau_7*tau_7*tau;
 
     RealOpenMM dtau         = tau/(dhal + one);
 
     RealOpenMM ratio        = (sigma_7/rho);
-    RealOpenMM gtau         = combindedEpsilon*tau_7*r_ij_6*(ghal+one)*ratio*ratio;
+    RealOpenMM gtau         = combinedEpsilon*tau_7*r_ij_6*(ghal+one)*ratio*ratio;
 
-    RealOpenMM energy       = combindedEpsilon*tau_7*sigma_7*( (ghal+one)*sigma_7/rho - two);
+    RealOpenMM energy       = combinedEpsilon*tau_7*sigma_7*( (ghal+one)*sigma_7/rho - two);
 
     RealOpenMM dEdR         = -seven*(dtau*energy + gtau);
 
@@ -242,9 +237,9 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma,
 
     dEdR                   /= r_ij;
 
-    force[0]                = dEdR*xr;
-    force[1]                = dEdR*yr;
-    force[2]                = dEdR*zr;
+    force[0]                = dEdR*deltaR[0];
+    force[1]                = dEdR*deltaR[1];
+    force[2]                = dEdR*deltaR[2];
 
     return energy;
 
@@ -315,11 +310,11 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
         for( unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++ ){
             if( exclusions[jj] == 0 ){
 
-                RealOpenMM combindedSigma   = (this->*_combineSigmas)(sigmaI, sigmas[jj] );
-                RealOpenMM combindedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj] );
+                RealOpenMM combinedSigma   = (this->*_combineSigmas)(sigmaI, sigmas[jj] );
+                RealOpenMM combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj] );
 
                 Vec3 force;
-                energy                     += calculatePairIxn( combindedSigma, combindedEpsilon,
+                energy                     += calculatePairIxn( combinedSigma, combinedEpsilon,
                                                                 reducedPositions[ii], reducedPositions[jj],
                                                                 force );
                 
@@ -384,11 +379,11 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
         int siteI                   = pair.first;
         int siteJ                   = pair.second;
 
-        RealOpenMM combindedSigma   = (this->*_combineSigmas)( sigmas[siteI], sigmas[siteJ] );
-        RealOpenMM combindedEpsilon = (this->*_combineEpsilons)( epsilons[siteI], epsilons[siteJ] );
+        RealOpenMM combinedSigma   = (this->*_combineSigmas)( sigmas[siteI], sigmas[siteJ] );
+        RealOpenMM combinedEpsilon = (this->*_combineEpsilons)( epsilons[siteI], epsilons[siteJ] );
 
         Vec3 force;
-        energy                     += calculatePairIxn( combindedSigma, combindedEpsilon,
+        energy                     += calculatePairIxn( combinedSigma, combinedEpsilon,
                                                         reducedPositions[siteI], reducedPositions[siteJ], force );
                 
         if( indexIVs[siteI] == siteI ){

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.h

@@ -31,7 +31,7 @@
 #include <string>
 #include <vector>
 
-using namespace OpenMM;
+namespace OpenMM {
 
 class AmoebaReferenceVdwForce;
 typedef  RealOpenMM (AmoebaReferenceVdwForce::*CombiningFunction)( RealOpenMM x, RealOpenMM y) const;
@@ -174,21 +174,11 @@ public:
     
        Set box dimensions
     
-       @param box box dimensions
+       @param vectors    the vectors defining the periodic box
     
        --------------------------------------------------------------------------------------- */
     
-    void setPeriodicBox( const RealVec& box );
-
-    /**---------------------------------------------------------------------------------------
-    
-       Get box dimensions
-    
-       @return box dimensions
-    
-       --------------------------------------------------------------------------------------- */
-    
-    RealVec getPeriodicBox( void ) const;
+    void setPeriodicBox(OpenMM::RealVec* vectors);
 
     /**---------------------------------------------------------------------------------------
     
@@ -253,7 +243,7 @@ private:
     double _taperCutoffFactor;
     double _taperCutoff;
     RealOpenMM _taperCoefficients[3];
-    RealVec _periodicBoxDimensions;
+    RealVec _periodicBoxVectors[3];
     CombiningFunction _combineSigmas;
     RealOpenMM arithmeticSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
     RealOpenMM  geometricSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
@@ -333,6 +323,7 @@ private:
 
 };
 
+}
 // ---------------------------------------------------------------------------------------
 
 #endif // _AmoebaReferenceVdwForce___

plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaVdwForce.cpp

@@ -39,6 +39,7 @@
 #include "openmm/System.h"
 #include "openmm/AmoebaVdwForce.h"
 #include "openmm/LangevinIntegrator.h"
+#include "sfmt/SFMT.h"
 #include <iostream>
 #include <vector>
 #include <stdlib.h>
@@ -50,6 +51,7 @@
 
 
 using namespace OpenMM;
+using namespace std;
 
 extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
 
@@ -1981,6 +1983,62 @@ void testVdwWater( int includeVdwDispersionCorrection, FILE* log ) {
     }
 }
 
+void testTriclinic() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    Vec3 a(3.1, 0, 0);
+    Vec3 b(0.4, 3.5, 0);
+    Vec3 c(-0.1, -0.5, 4.0);
+    system.setDefaultPeriodicBoxVectors(a, b, c);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    AmoebaVdwForce* vdw = new AmoebaVdwForce();
+    vdw->setUseDispersionCorrection(false);
+    vdw->addParticle(0, 0.5, 1.0, 0.0);
+    vdw->addParticle(1, 0.5, 1.0, 0.0);
+    vdw->setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
+    const double cutoff = 1.5;
+    vdw->setCutoff(cutoff);
+    system.addForce(vdw);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    vector<Vec3> positions(2);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int iteration = 0; iteration < 50; iteration++) {
+        // Generate random positions for the two particles.
+
+        positions[0] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        positions[1] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        context.setPositions(positions);
+
+        // Loop over all possible periodic copies and find the nearest one.
+
+        Vec3 delta;
+        double distance2 = 100.0;
+        for (int i = -1; i < 2; i++)
+            for (int j = -1; j < 2; j++)
+                for (int k = -1; k < 2; k++) {
+                    Vec3 d = positions[1]-positions[0]+a*i+b*j+c*k;
+                    if (d.dot(d) < distance2) {
+                        delta = d;
+                        distance2 = d.dot(d);
+                    }
+                }
+        double distance = sqrt(distance2);
+
+        // See if the energy is correct.
+
+        State state = context.getState(State::Energy);
+        if (distance >= cutoff) {
+            ASSERT_EQUAL(0.0, state.getPotentialEnergy());
+        }
+        else if (distance < 0.9*cutoff) {
+            const double energy = pow(1.07/(distance+0.07), 7.0)*(1.12/(pow(distance, 7.0)+0.12)-2);
+            ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+        }
+    }
+}
+
 int main( int numberOfArguments, char* argv[] ) {
 
     try {
@@ -2029,6 +2087,10 @@ int main( int numberOfArguments, char* argv[] ) {
  
         includeVdwDispersionCorrection     = 1;
         testVdwWater( includeVdwDispersionCorrection, log );
+        
+        // test triclinic boxes
+        
+        testTriclinic();
 
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;