Merge pull request #2378 from leucinw/master

Added the W-H combining rule for AMOEBA+ VdW

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

@@ -180,14 +180,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/openmmapi/src/AmoebaVdwForceImpl.cpp

@@ -192,7 +192,8 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
             // ARITHMETIC = 1
             // GEOMETRIC  = 2
             // HARMONIC   = 3
-            // HHG        = 4
+            // W-H        = 4
+            // HHG        = 5
             if (epsilonCombiningRule == "ARITHMETIC") {
               epsilon = 0.5f * (iEpsilon + jEpsilon);
             } else if (epsilonCombiningRule == "GEOMETRIC") {
@@ -203,6 +204,13 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
               } else {
                 epsilon = 0.0;
               }
+            } else if (epsilonCombiningRule == "W-H") {
+              double iSigma3 = iSigma * iSigma * iSigma;
+              double jSigma3 = jSigma * jSigma * jSigma;
+              double iSigma6 = iSigma3 * iSigma3;
+              double jSigma6 = jSigma3 * jSigma3;
+              double eps_s = std::sqrt(iEpsilon*jEpsilon);
+              epsilon = (eps_s == 0.0 ? 0.0 : 2.0f*eps_s*iSigma3*jSigma3/(iSigma6 + jSigma6));
             } else {
               double epsilonS = std::sqrt(iEpsilon) + std::sqrt(jEpsilon);
               if (epsilonS != 0.0) {

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

@@ -2442,8 +2442,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);
 

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

@@ -396,11 +396,12 @@ void testVdwAmmoniaCubicMeanHhg() {
     compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
 }
 
-// test alchemical VDW 
 
-void testVdwAlchemical(int power, double alpha, double lambda, AmoebaVdwForce::AlchemicalMethod method) {
+// test VDW w/ sigmaRule=CubicMean and epsilonRule=W-H
 
-    std::string testName      = "testVdwAlchemical";
+void testVdwAmmoniaCubicMeanWH() {
+
+    std::string testName      = "testVdwAmmoniaCubicMeanWH";
 
     int numberOfParticles     = 8;
     double boxDimension       = -1.0;
@@ -408,6 +409,32 @@ void testVdwAlchemical(int power, double alpha, double lambda, AmoebaVdwForce::A
     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 alchemical VDW 
+
+void testVdwAlchemical(int power, double alpha, double lambda, AmoebaVdwForce::AlchemicalMethod method) {
+
+    std::string testName      = "testVdwAlchemical";
+
     setupAndGetForcesEnergyVdwAmmonia2("CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy,
                                       method, power, alpha, lambda);
     std::vector<Vec3> expectedForces(numberOfParticles);
@@ -2090,6 +2117,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

@@ -155,6 +155,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 {
@@ -166,19 +168,27 @@ std::string AmoebaReferenceVdwForce::getEpsilonCombiningRule() const {
     return _epsilonCombiningRule;
 }
 
-double AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ) const {
+double AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const {
     return 0.5*(epsilonI + epsilonJ);
 }
 
-double AmoebaReferenceVdwForce::geometricEpsilonCombiningRule(double epsilonI, double epsilonJ) const {
+double AmoebaReferenceVdwForce::geometricEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const {
     return sqrt(epsilonI*epsilonJ);
 }
 
-double AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ) const {
+double AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const {
     return (epsilonI != 0.0 && epsilonJ != 0.0) ? 2.0*(epsilonI*epsilonJ)/(epsilonI + epsilonJ) : 0.0;
 }
+double AmoebaReferenceVdwForce::whEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) 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 AmoebaReferenceVdwForce::hhgEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const {
     double denominator = sqrt(epsilonI) + sqrt(epsilonJ);
     return (epsilonI != 0.0 && epsilonJ != 0.0) ? 4.0*(epsilonI*epsilonJ)/(denominator*denominator) : 0.0;
 }
@@ -310,7 +320,10 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
             if (exclusions[jj] == 0) {
 
                 double combinedSigma   = (this->*_combineSigmas)(sigmaI, sigmas[jj]);
-                double combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj]);
+
+                double combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj], sigmaI, sigmas[jj]);
+
+              
                 double softcore = 0.0;
 
                 if (this->_alchemicalMethod == Decouple && (isAlchemicalI != isAlchemical[jj])) {
@@ -321,6 +334,7 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
                    softcore = this->_alpha * pow(1.0 - lambda, 2);
                 }
 
+
                 Vec3 force;
                 energy += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
                                            reducedPositions[ii], reducedPositions[jj], force);
@@ -379,7 +393,9 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
         int siteJ                   = pair.second;
 
         double combinedSigma   = (this->*_combineSigmas)(sigmas[siteI], sigmas[siteJ]);
-        double combinedEpsilon = (this->*_combineEpsilons)(epsilons[siteI], epsilons[siteJ]);
+
+        double combinedEpsilon = (this->*_combineEpsilons)(epsilons[siteI], epsilons[siteJ], sigmas[siteI], sigmas[siteJ]);
+
         double softcore        = 0.0;
         int isAlchemicalI      = isAlchemical[siteI];
         int isAlchemicalJ      = isAlchemical[siteJ];
@@ -392,6 +408,7 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
            softcore = this->_alpha * pow(1.0 - lambda, 2);
         }
 
+
         Vec3 force;
         energy += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
                                    reducedPositions[siteI], reducedPositions[siteJ], force);

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

@@ -34,6 +34,7 @@ namespace OpenMM {
 
 class AmoebaReferenceVdwForce;
 typedef double (AmoebaReferenceVdwForce::*CombiningFunction)(double x, double y) const;
+typedef double (AmoebaReferenceVdwForce::*CombiningFunctionEpsilon)(double x, double y, double z, double w) const;
 
 // ---------------------------------------------------------------------------------------
 
@@ -333,13 +334,16 @@ private:
     Vec3 _periodicBoxVectors[3];
     CombiningFunction _combineSigmas;
     double arithmeticSigmaCombiningRule(double sigmaI, double sigmaJ) const;
-    double geometricSigmaCombiningRule(double sigmaI, double sigmaJ) const;
-    double cubicMeanSigmaCombiningRule(double sigmaI, double sigmaJ) const;
-    CombiningFunction _combineEpsilons;
-    double arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
-    double geometricEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
-    double harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
-    double hhgEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
+    double  geometricSigmaCombiningRule(double sigmaI, double sigmaJ) const;
+    double  cubicMeanSigmaCombiningRule(double sigmaI, double sigmaJ) const;
+
+    CombiningFunctionEpsilon _combineEpsilons;
+    double arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const;
+    double  geometricEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const;
+    double  harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const;
+    double  whEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const;
+    double  hhgEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const;
+
 
     /**---------------------------------------------------------------------------------------
     

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

@@ -405,6 +405,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 alchemical VDW 
 
 void testVdwAlchemical(int power, double alpha, double lambda, AmoebaVdwForce::AlchemicalMethod method) {
@@ -441,7 +471,6 @@ void testVdwAlchemical(int power, double alpha, double lambda, AmoebaVdwForce::A
     compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
 }
 
-
 // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
 
 void testVdwAmmoniaArithmeticArithmetic() {
@@ -2100,6 +2129,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)