Added W-H combining rule

plugins/amoeba/openmmapi/include/openmm/AmoebaVdwForce.h

@@ -136,14 +136,14 @@ public:
     /**
      * Set epsilon combining rule
      *
-     * @param epsilonCombiningRule   epsilon combining rule:   'ARITHMETIC', 'GEOMETRIC'. 'HARMONIC', 'HHG'
+     * @param epsilonCombiningRule   epsilon combining rule:   'ARITHMETIC', 'GEOMETRIC'. 'HARMONIC', 'W-H', 'HHG'
      */
     void setEpsilonCombiningRule(const std::string& epsilonCombiningRule);
 
     /**
      * Get epsilon combining rule
      *
-     * @return epsilonCombiningRule   epsilon combining rule:  'ARITHMETIC', 'GEOMETRIC'. 'HARMONIC', 'HHG'
+     * @return epsilonCombiningRule   epsilon combining rule:  'ARITHMETIC', 'GEOMETRIC'. 'HARMONIC', 'W-H', 'HHG'
      */
     const std::string& getEpsilonCombiningRule(void) const;
 

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp

@@ -2419,8 +2419,10 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
         replacements["EPSILON_COMBINING_RULE"] = "2";
     else if (epsilonCombiningRule == "HARMONIC")
         replacements["EPSILON_COMBINING_RULE"] = "3";
-    else if (epsilonCombiningRule == "HHG")
+    else if (epsilonCombiningRule == "W-H")
         replacements["EPSILON_COMBINING_RULE"] = "4";
+    else if (epsilonCombiningRule == "HHG")
+        replacements["EPSILON_COMBINING_RULE"] = "5";
     else
         throw OpenMMException("Illegal combining rule for sigma: "+sigmaCombiningRule);
     double cutoff = force.getCutoffDistance();

plugins/amoeba/platforms/cuda/src/kernels/amoebaVdwForce2.cu

@@ -22,6 +22,13 @@
 #elif EPSILON_COMBINING_RULE == 3
     real epssum = sigmaEpsilon1.y+sigmaEpsilon2.y;
     real epsilon = (epssum == 0.0f ? (real) 0 : 2*(sigmaEpsilon1.y*sigmaEpsilon2.y)/(sigmaEpsilon1.y+sigmaEpsilon2.y));
+#elif EPSILON_COMBINING_RULE == 4
+    real sigma1_3 = sigmaEpsilon1.x*sigmaEpsilon1.x*sigmaEpsilon1.x;
+    real sigma2_3 = sigmaEpsilon2.x*sigmaEpsilon2.x*sigmaEpsilon2.x;
+    real sigma1_6 = sigma1_3*sigma1_3; 
+    real sigma2_6 = sigma2_3*sigma2_3;
+    real eps_s = SQRT(sigmaEpsilon1.y*sigmaEpsilon2.y);
+    real epsilon = (eps_s == 0.0f ? (real) 0 : 2*eps_s*sigma1_3*sigma2_3/(sigma1_6 + sigma2_6));
 #else
     real epsilon_s = SQRT(sigmaEpsilon1.y) + SQRT(sigmaEpsilon2.y);
     real epsilon = (epsilon_s == 0.0f ? (real) 0 : 4*sigmaEpsilon1.y*sigmaEpsilon2.y/(epsilon_s*epsilon_s));

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

@@ -368,6 +368,36 @@ void testVdwAmmoniaCubicMeanHhg() {
     compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
 }
 
+// test VDW w/ sigmaRule=CubicMean and epsilonRule=W-H
+
+void testVdwAmmoniaCubicMeanWH() {
+
+    std::string testName      = "testVdwAmmoniaCubicMeanWH";
+
+    int numberOfParticles     = 8;
+    double boxDimension       = -1.0;
+    double cutoff             = 9000000.0;
+    std::vector<Vec3> forces;
+    double energy;
+
+    setupAndGetForcesEnergyVdwAmmonia("CUBIC-MEAN", "W-H", cutoff, boxDimension, forces, energy);
+    std::vector<Vec3> expectedForces(numberOfParticles);
+
+    double expectedEnergy     =  3.771794e+00;
+
+    expectedForces[0]         = Vec3( 2.3979839e+02,  -1.1829842e-02,  -5.3258772e+00);
+    expectedForces[1]         = Vec3(-1.9942459e+00,   4.3142144e-01,  -1.7290171e-01);
+    expectedForces[2]         = Vec3(-1.9935442e+00,  -4.2937965e-01,  -1.7369876e-01);
+    expectedForces[3]         = Vec3(-8.9050582e-01,  -4.8659920e-04,   2.2190848e-01);
+    expectedForces[4]         = Vec3(-5.2306326e+01,   1.4895040e-03,  -3.3588483e-01);
+    expectedForces[5]         = Vec3( 1.4153288e+00,  -2.7130186e-01,  -1.5480591e-01);
+    expectedForces[6]         = Vec3(-1.8544507e+02,   8.4027272e-03,   6.0950274e+00);
+    expectedForces[7]         = Vec3( 1.4159723e+00,   2.7168386e-01,  -1.5376786e-01);
+
+    double tolerance          = 1.0e-04;
+    compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
+}
+
 // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
 
 void testVdwAmmoniaArithmeticArithmetic() {
@@ -2025,6 +2055,10 @@ int main(int argc, char* argv[]) {
 
         testVdwAmmoniaCubicMeanHhg();
 
+	// test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG
+
+        testVdwAmmoniaCubicMeanWH();
+
         // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
 
         testVdwAmmoniaArithmeticArithmetic();

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

@@ -129,6 +129,8 @@ void AmoebaReferenceVdwForce::setEpsilonCombiningRule(const std::string& epsilon
          _combineEpsilons = &AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule;
     } else if (_epsilonCombiningRule == "HARMONIC") {
          _combineEpsilons = &AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule;
+    } else if (_epsilonCombiningRule == "W-H") {
+         _combineEpsilons = &AmoebaReferenceVdwForce::whEpsilonCombiningRule;
     } else if (_epsilonCombiningRule == "HHG") {
          _combineEpsilons = &AmoebaReferenceVdwForce::hhgEpsilonCombiningRule;
     } else {
@@ -151,6 +153,14 @@ double AmoebaReferenceVdwForce::geometricEpsilonCombiningRule(double epsilonI, d
 double AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ) const {
     return (epsilonI != 0.0 && epsilonJ != 0.0) ? 2.0*(epsilonI*epsilonJ)/(epsilonI + epsilonJ) : 0.0;
 }
+double AmoebaReferenceVdwForce::whEpsilonCombiningRule(double sigmaI, double sigmaJ, double epsilonI, double epsilonJ) const {
+    double  sigmaI3 = sigmaI * sigmaI * sigmaI;
+    double  sigmaJ3 = sigmaJ * sigmaJ * sigmaJ;
+    double  sigmaI6 = sigmaI3 * sigmaI3;
+    double  sigmaJ6 = sigmaJ3 * sigmaJ3;
+    double  eps_s = sqrt(epsilonI*epsilonJ);
+    return (epsilonI != 0.0 && epsilonJ != 0.0) ? 2.0*eps_s*sigmaI3*sigmaJ3/(sigmaI6+sigmaJ6) : 0.0;
+}
 
 double AmoebaReferenceVdwForce::hhgEpsilonCombiningRule(double epsilonI, double epsilonJ) const {
     double denominator = sqrt(epsilonI) + sqrt(epsilonJ);

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

@@ -252,6 +252,7 @@ private:
     double arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
     double  geometricEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
     double  harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
+    double  whEpsilonCombiningRule(double sigmaI, double sigmaJ, double epsilonI, double epsilonJ) const;
     double  hhgEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
 
     /**---------------------------------------------------------------------------------------

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

@@ -379,6 +379,35 @@ void testVdwAmmoniaCubicMeanHhg() {
     compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
 }
 
+// test VDW w/ sigmaRule=CubicMean and epsilonRule=W-H
+void testVdwAmmoniaCubicMeanWH() {
+
+    std::string testName      = "testVdwAmmoniaCubicMeanWH";
+
+    int numberOfParticles     = 8;
+    double boxDimension       = -1.0;
+    double cutoff             = 9000000.0;
+    std::vector<Vec3> forces;
+    double energy;
+
+    setupAndGetForcesEnergyVdwAmmonia("CUBIC-MEAN", "W-H", cutoff, boxDimension, forces, energy);
+    std::vector<Vec3> expectedForces(numberOfParticles);
+
+    double expectedEnergy     =  3.771794e+00;
+
+    expectedForces[0]         = Vec3( 2.3979839e+02,  -1.1829842e-02,  -5.3258772e+00);
+    expectedForces[1]         = Vec3(-1.9942459e+00,   4.3142144e-01,  -1.7290171e-01);
+    expectedForces[2]         = Vec3(-1.9935442e+00,  -4.2937965e-01,  -1.7369876e-01);
+    expectedForces[3]         = Vec3(-8.9050582e-01,  -4.8659920e-04,   2.2190848e-01);
+    expectedForces[4]         = Vec3(-5.2306326e+01,   1.4895040e-03,  -3.3588483e-01);
+    expectedForces[5]         = Vec3( 1.4153288e+00,  -2.7130186e-01,  -1.5480591e-01);
+    expectedForces[6]         = Vec3(-1.8544507e+02,   8.4027272e-03,   6.0950274e+00);
+    expectedForces[7]         = Vec3( 1.4159723e+00,   2.7168386e-01,  -1.5376786e-01);
+
+    double tolerance          = 1.0e-04;
+    compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
+}
+
 // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
 
 void testVdwAmmoniaArithmeticArithmetic() {
@@ -2037,6 +2066,10 @@ int main(int numberOfArguments, char* argv[]) {
 
         testVdwAmmoniaCubicMeanHhg();
 
+	// test VDW w/ sigmaRule=CubicMean and epsilonRule=W-H
+
+        testVdwAmmoniaCubicMeanWH();
+
         // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
 
         testVdwAmmoniaArithmeticArithmetic();

wrappers/python/simtk/openmm/app/forcefield.py

@@ -4431,7 +4431,7 @@ class AmoebaVdwGenerator(object):
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
 
         sigmaMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'CUBIC-MEAN':1}
-        epsilonMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'HARMONIC':1, 'HHG':1}
+        epsilonMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'HARMONIC':1, 'W-H':1, 'HHG':1}
 
         force = mm.AmoebaVdwForce()
         sys.addForce(force)