Cleanup to reference implementation of AmoebaVdwForce

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.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-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors:                                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -1000,9 +1000,8 @@ ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(const std::
 }
 
 ReferenceCalcAmoebaVdwForceKernel::~ReferenceCalcAmoebaVdwForceKernel() {
-    if (neighborList) {
+    if (neighborList)
         delete neighborList;
-    } 
 }
 
 void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) {
@@ -1010,43 +1009,11 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
     // per-particle parameters
 
     numParticles = system.getNumParticles();
-
-    indexIVs.resize(numParticles);
-    allExclusions.resize(numParticles);
-    sigmas.resize(numParticles);
-    epsilons.resize(numParticles);
-    reductions.resize(numParticles);
-    isAlchemical.resize(numParticles);
-
-    for (int ii = 0; ii < numParticles; ii++) {
-
-        int indexIV, type;
-        double sigma, epsilon, reduction;
-        bool alchemical;
-        std::vector<int> exclusions;
-
-        force.getParticleParameters(ii, indexIV, sigma, epsilon, reduction, alchemical, type);
-        force.getParticleExclusions(ii, exclusions);
-        for (unsigned int jj = 0; jj < exclusions.size(); jj++) {
-           allExclusions[ii].insert(exclusions[jj]);
-        }
-
-        indexIVs[ii]      = indexIV;
-        sigmas[ii]        = sigma;
-        epsilons[ii]      = epsilon;
-        reductions[ii]    = reduction;
-        isAlchemical[ii]  = alchemical;
-    }   
-    sigmaCombiningRule     = force.getSigmaCombiningRule();
-    epsilonCombiningRule   = force.getEpsilonCombiningRule();
     useCutoff              = (force.getNonbondedMethod() != AmoebaVdwForce::NoCutoff);
     usePBC                 = (force.getNonbondedMethod() == AmoebaVdwForce::CutoffPeriodic);
     cutoff                 = force.getCutoffDistance();
     neighborList           = useCutoff ? new NeighborList() : NULL;
     dispersionCoefficient  = force.getUseDispersionCorrection() ?  AmoebaVdwForceImpl::calcDispersionCorrection(system, force) : 0.0;
-    alchemicalMethod       = force.getAlchemicalMethod();
-    n                      = force.getSoftcorePower();
-    alpha                  = force.getSoftcoreAlpha();
     vdwForce.initialize(force);
 }
 
@@ -1057,41 +1024,25 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
     double energy;
     double lambda = context.getParameter(AmoebaVdwForce::Lambda());
     if (useCutoff) {
-        computeNeighborListVoxelHash(*neighborList, numParticles, posData, allExclusions, extractBoxVectors(context), usePBC, cutoff, 0.0);
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, vdwForce.getExclusions(), extractBoxVectors(context), usePBC, cutoff, 0.0);
         if (usePBC) {
             Vec3* boxVectors = extractBoxVectors(context);
             double minAllowedSize = 1.999999*cutoff;
-            if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][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(boxVectors);
-            energy  = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, indexIVs, sigmas, epsilons, 
-                                                       reductions, isAlchemical, *neighborList, forceData);
+            energy  = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, *neighborList, forceData);
             energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
         }
     }
     else
-        energy = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, indexIVs, sigmas, epsilons, 
-                                                  reductions, isAlchemical, allExclusions, forceData);
+        energy = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, forceData);
     return static_cast<double>(energy);
 }
 
 void ReferenceCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaVdwForce& force) {
     if (numParticles != force.getNumParticles())
         throw OpenMMException("updateParametersInContext: The number of particles has changed");
-
-    // Record the values.
-    for (int i = 0; i < numParticles; ++i) {
-        int indexIV, type;
-        double sigma, epsilon, reduction;
-        bool alchemical;
-        force.getParticleParameters(i, indexIV, sigma, epsilon, reduction, alchemical, type);
-        indexIVs[i] = indexIV;
-        sigmas[i] = sigma;
-        epsilons[i] = epsilon;
-        reductions[i]= reduction;
-        isAlchemical[i]= alchemical;
-    }
     vdwForce.initialize(force);
 }
 

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.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) 2008-2018 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors:                                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -503,17 +503,6 @@ private:
     double cutoff;
     double dispersionCoefficient;
     AmoebaReferenceVdwForce vdwForce;
-    AmoebaVdwForce::AlchemicalMethod alchemicalMethod;
-    int n;
-    double alpha;
-    std::vector<int> indexIVs;
-    std::vector< std::set<int> > allExclusions;
-    std::vector<double> sigmas;
-    std::vector<double> epsilons;
-    std::vector<double> reductions;
-    std::vector<bool> isAlchemical;
-    std::string sigmaCombiningRule;
-    std::string epsilonCombiningRule;
     const System& system;
     NeighborList* neighborList;
 };

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

@@ -34,10 +34,7 @@ using std::vector;
 using namespace OpenMM;
 
 AmoebaReferenceVdwForce::AmoebaReferenceVdwForce() : _nonbondedMethod(AmoebaVdwForce::NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9), _n(5), _alpha(0.7), _alchemicalMethod(AmoebaVdwForce::None) {
-
     setTaperCoefficients(_cutoff);
-    setSigmaCombiningRule("ARITHMETIC");
-    setEpsilonCombiningRule("GEOMETRIC");
 }
 
 void AmoebaReferenceVdwForce::initialize(const AmoebaVdwForce& force) {
@@ -49,32 +46,25 @@ void AmoebaReferenceVdwForce::initialize(const AmoebaVdwForce& force) {
     if (force.getNonbondedMethod() != AmoebaVdwForce::NoCutoff)
         setCutoff(force.getCutoffDistance());
     AmoebaVdwForceImpl::createParameterMatrix(force, particleType, sigmaMatrix, epsilonMatrix);
+    int numParticles = force.getNumParticles();
+    indexIVs.resize(numParticles);
+    reductions.resize(numParticles);
+    isAlchemical.resize(numParticles);
+    allExclusions.clear();
+    allExclusions.resize(numParticles);
+    for (int i = 0; i < numParticles; i++) {
+        int type;
+        double sigma, epsilon;
+        bool alchemical;
+        vector<int> exclusions;
+        force.getParticleParameters(i, indexIVs[i], sigma, epsilon, reductions[i], alchemical, type);
+        isAlchemical[i] = alchemical;
+        force.getParticleExclusions(i, exclusions);
+        for (unsigned int j = 0; j < exclusions.size(); j++)
+           allExclusions[i].insert(exclusions[j]);
+    }
 }
 
-AmoebaReferenceVdwForce::AmoebaReferenceVdwForce(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule) : _nonbondedMethod(AmoebaVdwForce::NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9),  _n(5), _alpha(0.7), _alchemicalMethod(AmoebaVdwForce::None) {
-
-    setTaperCoefficients(_cutoff);
-    setSigmaCombiningRule(sigmaCombiningRule);
-    setEpsilonCombiningRule(epsilonCombiningRule);
-}
-
-void AmoebaReferenceVdwForce::setSoftcorePower(int n) {
-    _n = n;
-}
-
-int AmoebaReferenceVdwForce::getSoftcorePower() const {
-    return _n;
-}
-
-void AmoebaReferenceVdwForce::setSoftcoreAlpha(double alpha) {
-    _alpha = alpha;
-}
-
-double AmoebaReferenceVdwForce::getSoftcoreAlpha() const {
-    return _alpha;
-}
-
-
 void AmoebaReferenceVdwForce::setTaperCoefficients(double cutoff) {
     _taperCutoff = cutoff*_taperCutoffFactor;
     if (_taperCutoff != cutoff) {
@@ -93,98 +83,14 @@ void AmoebaReferenceVdwForce::setCutoff(double cutoff) {
     setTaperCoefficients(_cutoff);
 }
 
-double AmoebaReferenceVdwForce::getCutoff() const {
-    return _cutoff;
-}
-
 void AmoebaReferenceVdwForce::setPeriodicBox(OpenMM::Vec3* vectors) {
     _periodicBoxVectors[0] = vectors[0];
     _periodicBoxVectors[1] = vectors[1];
     _periodicBoxVectors[2] = vectors[2];
 }
 
-void AmoebaReferenceVdwForce::setSigmaCombiningRule(const std::string& sigmaCombiningRule) {
-
-    _sigmaCombiningRule = sigmaCombiningRule;
-
-    // convert to upper case and set combining function
-
-    std::transform(_sigmaCombiningRule.begin(), _sigmaCombiningRule.end(), _sigmaCombiningRule.begin(),  (int(*)(int)) std::toupper);
-    if (_sigmaCombiningRule == "GEOMETRIC") {
-        _combineSigmas = &AmoebaReferenceVdwForce::geometricSigmaCombiningRule;
-    } else if (_sigmaCombiningRule == "CUBIC-MEAN") {
-        _combineSigmas = &AmoebaReferenceVdwForce::cubicMeanSigmaCombiningRule;
-    } else {
-        _combineSigmas = &AmoebaReferenceVdwForce::arithmeticSigmaCombiningRule;
-    }
-}
-
-std::string AmoebaReferenceVdwForce::getSigmaCombiningRule() const {
-    return _sigmaCombiningRule;
-}
-
-double AmoebaReferenceVdwForce::arithmeticSigmaCombiningRule(double sigmaI, double sigmaJ) const {
-    return (sigmaI + sigmaJ);
-}
-
-double AmoebaReferenceVdwForce::geometricSigmaCombiningRule(double sigmaI, double sigmaJ) const {
-    return 2.0*sqrt(sigmaI*sigmaJ);
-}
-
-double AmoebaReferenceVdwForce::cubicMeanSigmaCombiningRule(double sigmaI, double sigmaJ) const {
-    double sigmaI2 = sigmaI*sigmaI;
-    double sigmaJ2 = sigmaJ*sigmaJ;
-
-    return sigmaI != 0.0 && sigmaJ != 0.0 ? 2.0*(sigmaI2*sigmaI + sigmaJ2*sigmaJ)/(sigmaI2 + sigmaJ2) : 0.0;
-}
-
-void AmoebaReferenceVdwForce::setEpsilonCombiningRule(const std::string& epsilonCombiningRule) {
-
-    _epsilonCombiningRule = epsilonCombiningRule;
-    std::transform(_epsilonCombiningRule.begin(), _epsilonCombiningRule.end(), _epsilonCombiningRule.begin(),  (int(*)(int)) std::toupper);
-
-    // convert to upper case and set combining function
-
-    if (_epsilonCombiningRule == "ARITHMETIC") {
-         _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 {
-         _combineEpsilons = &AmoebaReferenceVdwForce::geometricEpsilonCombiningRule;
-    }
-}
-
-std::string AmoebaReferenceVdwForce::getEpsilonCombiningRule() const {
-    return _epsilonCombiningRule;
-}
-
-double AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const {
-    return 0.5*(epsilonI + epsilonJ);
-}
-
-double AmoebaReferenceVdwForce::geometricEpsilonCombiningRule(double epsilonI, double epsilonJ, double sigmaI, double sigmaJ) const {
-    return sqrt(epsilonI*epsilonJ);
-}
-
-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, 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;
+std::vector<std::set<int> >& AmoebaReferenceVdwForce::getExclusions() {
+    return allExclusions;
 }
 
 void AmoebaReferenceVdwForce::addReducedForce(unsigned int particleI, unsigned int particleIV,
@@ -217,8 +123,6 @@ double AmoebaReferenceVdwForce::calculatePairIxn(double combinedSigma, double co
         ReferenceForce::getDeltaRPeriodic(particleJPosition, particleIPosition, _periodicBoxVectors, deltaR);
     else
         ReferenceForce::getDeltaR(particleJPosition, particleIPosition, deltaR);
-
-    double r_ij_2       = deltaR[ReferenceForce::R2Index];
     double r_ij         = deltaR[ReferenceForce::RIndex];
 
     double energy, dEdR;
@@ -282,20 +186,14 @@ void AmoebaReferenceVdwForce::setReducedPositions(int numParticles,
             reducedPositions[ii]   = Vec3(reductions[ii]*(particlePositions[ii][0] - particlePositions[reductionIndex][0]) + particlePositions[reductionIndex][0], 
                                           reductions[ii]*(particlePositions[ii][1] - particlePositions[reductionIndex][1]) + particlePositions[reductionIndex][1], 
                                           reductions[ii]*(particlePositions[ii][2] - particlePositions[reductionIndex][2]) + particlePositions[reductionIndex][2]); 
-        } else {
-            reducedPositions[ii]   = Vec3(particlePositions[ii][0], particlePositions[ii][1], particlePositions[ii][2]); 
         }
+        else
+            reducedPositions[ii]   = Vec3(particlePositions[ii][0], particlePositions[ii][1], particlePositions[ii][2]); 
     }
 }
 
 double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double lambda,
                                                         const vector<Vec3>& particlePositions,
-                                                        const std::vector<int>& indexIVs, 
-                                                        const std::vector<double>& sigmas,
-                                                        const std::vector<double>& epsilons,
-                                                        const std::vector<double>& reductions,
-                                                        const std::vector<bool>& isAlchemical,
-                                                        const std::vector< std::set<int> >& allExclusions,
                                                         vector<Vec3>& forces) const {
 
     // set reduced coordinates
@@ -365,11 +263,6 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
 
 double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double lambda,
                                                         const vector<Vec3>& particlePositions,
-                                                        const std::vector<int>& indexIVs, 
-                                                        const std::vector<double>& sigmas,
-                                                        const std::vector<double>& epsilons,
-                                                        const std::vector<double>& reductions,
-                                                        const std::vector<bool>& isAlchemical,
                                                         const NeighborList& neighborList,
                                                         vector<Vec3>& forces) const {
 

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

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -28,6 +28,7 @@
 #include "openmm/Vec3.h"
 #include "openmm/AmoebaVdwForce.h"
 #include "ReferenceNeighborList.h"
+#include <set>
 #include <string>
 #include <vector>
 
@@ -54,36 +55,8 @@ public:
     void initialize(const AmoebaVdwForce& force);
  
     /**---------------------------------------------------------------------------------------
-       
-       Constructor
-       
-       --------------------------------------------------------------------------------------- */
- 
-    AmoebaReferenceVdwForce(const std::string& sigmaCombiningRule,
-                            const std::string& epsilonCombiningRule);
- 
-    /**---------------------------------------------------------------------------------------
-       
-       Destructor
-       
-       --------------------------------------------------------------------------------------- */
- 
-    ~AmoebaReferenceVdwForce() {};
- 
-
-    /**---------------------------------------------------------------------------------------
-    
-       Get cutoff
-    
-       @return cutoff
-    
-       --------------------------------------------------------------------------------------- */
     
-    double getCutoff() const;
-
-    /**---------------------------------------------------------------------------------------
-    
-       Set cutof
+       Set cutoff
     
        @param cutoff
     
@@ -92,95 +65,22 @@ public:
     void setCutoff(double cutoff);
 
     /**---------------------------------------------------------------------------------------
-   
-       Get softcore power
-   
-       @return n
-   
-       --------------------------------------------------------------------------------------- */
-   
-    int getSoftcorePower() const;
-
-    /**---------------------------------------------------------------------------------------
-   
-       Set softcore power
-   
-       @param n
-
-       --------------------------------------------------------------------------------------- */
-
-    void setSoftcorePower(int n);
-
-   /**---------------------------------------------------------------------------------------
-
-       Get softcore alpha
-  
-       @return alpha
-  
-       --------------------------------------------------------------------------------------- */
-  
-    double getSoftcoreAlpha() const;
-
-    /**---------------------------------------------------------------------------------------
-  
-       Set softcore alpha 
-  
-       @param alpha
-
-       --------------------------------------------------------------------------------------- */
-
-    void setSoftcoreAlpha(double alpha);
-
-
-    /**---------------------------------------------------------------------------------------
-    
-       Set sigma combining rule
-    
-       @param sigmaCombiningRule      rule: GEOMETRIC, CUBIC-MEAN, ARITHMETIC (default)
-    
-       --------------------------------------------------------------------------------------- */
-    
-    void setSigmaCombiningRule(const std::string& sigmaCombiningRule);
-
-    /**---------------------------------------------------------------------------------------
-    
-       Get sigma combining rule
-    
-       @return sigmaCombiningRule
-    
-       --------------------------------------------------------------------------------------- */
-    
-    std::string getSigmaCombiningRule() const;
-
-    /**---------------------------------------------------------------------------------------
-    
-       Set epsilon combining rule
-    
-       @param epsilonCombiningRule      rule: GEOMETRIC, CUBIC-MEAN, ARITHMETIC (default)
     
-       --------------------------------------------------------------------------------------- */
-    
-    void setEpsilonCombiningRule(const std::string& epsilonCombiningRule);
-
-    /**---------------------------------------------------------------------------------------
-    
-       Get epsilon combining rule
+       Set box dimensions
     
-       @return epsilonCombiningRule
+       @param vectors    the vectors defining the periodic box
     
        --------------------------------------------------------------------------------------- */
     
-    std::string getEpsilonCombiningRule() const;
-
+    void setPeriodicBox(OpenMM::Vec3* vectors);
+ 
     /**---------------------------------------------------------------------------------------
     
-       Set box dimensions
-    
-       @param vectors    the vectors defining the periodic box
+       Get the set of exclusions for each particle.
     
        --------------------------------------------------------------------------------------- */
     
-    void setPeriodicBox(OpenMM::Vec3* vectors);
+    std::vector<std::set<int> >& getExclusions();
 
     /**---------------------------------------------------------------------------------------
     
@@ -189,12 +89,6 @@ public:
        @param numParticles            number of particles
        @param lambda                  lambda value
        @param particlePositions       Cartesian coordinates of particles
-       @param indexIVs                position index for associated reducing particle
-       @param sigmas                  particle sigmas 
-       @param epsilons                particle epsilons
-       @param reductions              particle reduction factors
-       @param isAlchemical            particle alchemical flag
-       @param vdwExclusions           particle exclusions
        @param forces                  add forces to this vector
     
        @return energy
@@ -202,11 +96,6 @@ public:
        --------------------------------------------------------------------------------------- */
     
     double calculateForceAndEnergy(int numParticles, double lambda, const std::vector<OpenMM::Vec3>& particlePositions,
-                                   const std::vector<int>& indexIVs, 
-                                   const std::vector<double>& sigmas, const std::vector<double>& epsilons,
-                                   const std::vector<double>& reductions,
-                                   const std::vector<bool>& isAlchemical,
-                                   const std::vector< std::set<int> >& vdwExclusions,
                                    std::vector<OpenMM::Vec3>& forces) const;
          
     /**---------------------------------------------------------------------------------------
@@ -216,11 +105,6 @@ public:
        @param numParticles            number of particles
        @param lambda                  lambda value
        @param particlePositions       Cartesian coordinates of particles
-       @param indexIVs                position index for associated reducing particle
-       @param sigmas                  particle sigmas 
-       @param epsilons                particle epsilons
-       @param reductions              particle reduction factors
-       @param isAlchemical            particle alchemical flag
        @param neighborList            neighbor list
        @param forces                  add forces to this vector
     
@@ -229,12 +113,7 @@ public:
        --------------------------------------------------------------------------------------- */
     
     double calculateForceAndEnergy(int numParticles, double lambda, const std::vector<OpenMM::Vec3>& particlePositions, 
-                                   const std::vector<int>& indexIVs, 
-                                   const std::vector<double>& sigmas, const std::vector<double>& epsilons,
-                                   const std::vector<double>& reductions,
-                                   const std::vector<bool>& isAlchemical,
-                                   const NeighborList& neighborList,
-                                   std::vector<OpenMM::Vec3>& forces) const;
+                                   const NeighborList& neighborList, std::vector<OpenMM::Vec3>& forces) const;
          
 private:
     // taper coefficient indices
@@ -242,8 +121,6 @@ private:
     static const int C4=1;
     static const int C5=2;
 
-    std::string _sigmaCombiningRule;
-    std::string _epsilonCombiningRule;
     AmoebaVdwForce::NonbondedMethod _nonbondedMethod;
     AmoebaVdwForce::AlchemicalMethod _alchemicalMethod;
     AmoebaVdwForce::PotentialFunction potentialFunction;
@@ -256,19 +133,11 @@ private:
     std::vector<int> particleType;
     std::vector<std::vector<double> > sigmaMatrix;
     std::vector<std::vector<double> > epsilonMatrix;
+    std::vector<int> indexIVs;
+    std::vector<double> reductions;
+    std::vector<bool> isAlchemical;
+    std::vector<std::set<int> > allExclusions;
     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;
-
-    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;
-
 
     /**---------------------------------------------------------------------------------------