Merge pull request #2765 from peastman/amoeba

Added new features to AmoebaVdwForce

platforms/cuda/src/kernels/nonbonded.cu

@@ -18,6 +18,10 @@ static __inline__ __device__ float real_shfl(float var, int srcLane) {
     return SHFL(var, srcLane);
 }
 
+static __inline__ __device__ float real_shfl(int var, int srcLane) {
+    return SHFL(var, srcLane);
+}
+
 static __inline__ __device__ double real_shfl(double var, int srcLane) {
     int hi, lo;
     asm volatile("mov.b64 { %0, %1 }, %2;" : "=r"(lo), "=r"(hi) : "d"(var));

plugins/amoeba/openmmapi/include/openmm/AmoebaVdwForce.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-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs, Peter Eastman                                    *
  * Contributors:                                                              *
  *                                                                            *
@@ -34,16 +34,30 @@
 
 #include "openmm/Force.h"
 #include "internal/windowsExportAmoeba.h"
+#include <string>
 #include <vector>
 
 namespace OpenMM {
 
 /**
- * This class implements a buffered 14-7 potential used to model van der Waals forces.
+ * This class models van der Waals forces in the AMOEBA force field.  It can use
+ * either buffered 14-7 potential or a Lennard-Jones 12-6 potential.
  *
- * To use it, create an AmoebaVdwForce object then call addParticle() once for each particle.  After
- * a particle has been added, you can modify its force field parameters by calling setParticleParameters().
- * This will have no effect on Contexts that already exist unless you call updateParametersInContext().
+ * This class can operate in two different modes.  In one mode, force field parameters
+ * are defined for each particle.  When two particles interact, a combining rule is
+ * used to calculate the interaction parameters based on the parameters for the two
+ * particles.  To use the class in this mode, call the version of addParticle() that
+ * takes sigma and epsilon values.  It should be called once for each particle in the
+ * System.
+ * 
+ * In the other mode, each particle has a type index, and parameters are specified for
+ * each type rather than each individual particle.  By default this mode also uses a
+ * combining rule, but you can override it by defining alternate parameters to use for
+ * specific pairs of particle types.  To use the class in this mode, call the version of
+ * addParticle() that takes a type index.  It should be called once for each particle
+ * in the System.  You also must call addParticleType() once for each type.  If you
+ * wish to override the combining for particular pairs of types, do so by calling
+ * addTypePair().
  * 
  * A unique feature of this class is that the interaction site for a particle does not need to be
  * exactly at the particle's location.  Instead, it can be placed a fraction of the distance from that
@@ -91,6 +105,20 @@ public:
         CutoffPeriodic = 1,
     };
 
+    /**
+     * This is an enumeration of the different potential functions that can be used.
+     */
+    enum PotentialFunction {
+        /**
+         * Use a buffered 14-7 potential.  This is the default.
+         */
+        Buffered147 = 0,
+        /**
+         * Use a Lennard-Jones 12-6 potential.
+         */
+        LennardJones = 1,
+    };
+
     /**
      * This is an enumeration of the different alchemical methods used when applying softcore interactions.
      */
@@ -121,6 +149,20 @@ public:
         return parameters.size();
     }
 
+    /**
+     * Get the number of particle types.
+     */
+    int getNumParticleTypes() const {
+        return types.size();
+    }
+
+    /**
+     * Get the number of type pairs.
+     */
+    int getNumTypePairs() const {
+        return pairs.size();
+    }
+
     /**
      * Set the force field parameters for a vdw particle.
      *
@@ -131,9 +173,10 @@ public:
      * @param reductionFactor the fraction of the distance along the line from the parent particle to this particle
      *                        at which the interaction site should be placed
      * @param isAlchemical    if true, this vdW particle is undergoing an alchemical change.
+     * @param typeIndex       the index of the particle type for this particle
      */
     void setParticleParameters(int particleIndex, int parentIndex, double sigma, double epsilon, 
-                               double reductionFactor, bool isAlchemical = false);
+                               double reductionFactor, bool isAlchemical=false, int typeIndex=-1);
 
     /**
      * Get the force field parameters for a vdw particle.
@@ -145,13 +188,14 @@ public:
      * @param[out] reductionFactor the fraction of the distance along the line from the parent particle to this particle
      *                             at which the interaction site should be placed
      * @param[out] isAlchemical    if true, this vdW particle is undergoing an alchemical change.
+     * @param[out] typeIndex       the index of the particle type for this particle
      */
     void getParticleParameters(int particleIndex, int& parentIndex, double& sigma, double& epsilon, 
-                               double& reductionFactor, bool& isAlchemical) const;
-
+                               double& reductionFactor, bool& isAlchemical, int& typeIndex) const;
 
     /**
-     * Add the force field parameters for a vdw particle.
+     * Add the force field parameters for a vdw particle.  This version is used when parameters
+     * are defined for each particle.
      *
      * @param parentIndex     the index of the parent particle
      * @param sigma           vdw sigma
@@ -163,6 +207,82 @@ public:
      */
     int addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical = false);
 
+    /**
+     * Add the force field parameters for a vdw particle. This version is used when parameters
+     * are defined by particle type.
+     *
+     * @param parentIndex     the index of the parent particle
+     * @param typeIndex       the index of the particle type for this particle
+     * @param reductionFactor the fraction of the distance along the line from the parent particle to this particle
+     *                        at which the interaction site should be placed
+     * @param isAlchemical    if true, this vdW particle is undergoing an alchemical change.
+     * @return index of added particle
+     */
+    int addParticle(int parentIndex, int typeIndex, double reductionFactor, bool isAlchemical = false);
+
+    /**
+     * Add a particle type.
+     * 
+     * @param sigma     the sigma value for particles of this type
+     * @param epsilon   the epsilon value for particles of this type
+     * @return the index of the particle type that was just added.
+     */
+    int addParticleType(double sigma, double epsilon);
+
+    /**
+     * Get the force field parameters for a particle type.
+     * 
+     * @param typeIndex      the index of the particle type
+     * @param[out] sigma     the sigma value for particles of this type
+     * @param[out] epsilon   the epsilon value for particles of this type
+     */
+    void getParticleTypeParameters(int typeIndex, double& sigma, double& epsilon) const;
+
+    /**
+     * Set the force field parameters for a particle type.
+     * 
+     * @param typeIndex the index of the particle type
+     * @param sigma     the sigma value for particles of this type
+     * @param epsilon   the epsilon value for particles of this type
+     */
+    void setParticleTypeParameters(int typeIndex, double sigma, double epsilon);
+
+    /**
+     * Add a type pair.  This overrides the standard combining rule for interactions
+     * between particles of two particular types.
+     * 
+     * @param type1     the index of the first particle type
+     * @param type2     the index of the second particle type
+     * @param sigma     the sigma value for interactions between particles of these two types
+     * @param epsilon   the epsilon  value for interactions between particles of these two types
+     * @return the index of the type pair that was just added.
+     */
+    int addTypePair(int type1, int type2, double sigma, double epsilon);
+
+    /**
+     * Get the force field parameters for a type pair.  This overrides the standard
+     * combining rule for interactions between particles of two particular types.
+     * 
+     * @param pairIndex      the index of the type pair
+     * @param[out] type1     the index of the first particle type
+     * @param[out] type2     the index of the second particle type
+     * @param[out] sigma     the sigma value for interactions between particles of these two types
+     * @param[out] epsilon   the epsilon  value for interactions between particles of these two types
+     */
+    void getTypePairParameters(int pairIndex, int& type1, int& type2, double& sigma, double& epsilon) const;
+
+    /**
+     * Set the force field parameters for a type pair.  This overrides the standard
+     * combining rule for interactions between particles of two particular types.
+     * 
+     * @param pairIndex the index of the type pair
+     * @param type1     the index of the first particle type
+     * @param type2     the index of the second particle type
+     * @param sigma     the sigma value for interactions between particles of these two types
+     * @param epsilon   the epsilon  value for interactions between particles of these two types
+     */
+    void setTypePairParameters(int pairIndex, int type1, int type2, double sigma, double epsilon);
+    
     /**
      * Set sigma combining rule
      *
@@ -211,6 +331,13 @@ public:
         useDispersionCorrection = useCorrection;
     }
     
+    /**
+     * Get whether parameters were specified by particle or by particle type.
+     */
+    bool getUseParticleTypes() const {
+        return useTypes;
+    }
+
     /**
      * Set exclusions for specified particle
      *
@@ -268,6 +395,20 @@ public:
      */
     void setNonbondedMethod(NonbondedMethod method);
 
+    /**
+     * Get the potential function to use.
+     */
+    PotentialFunction getPotentialFunction() const {
+        return potentialFunction;
+    }
+
+    /**
+     * Set the potential function to use.
+     */
+    void setPotentialFunction(PotentialFunction potential) {
+        potentialFunction = potential;
+    }
+
     /**
      * Set the softcore power on lambda (default = 5).
      */
@@ -322,9 +463,12 @@ protected:
 private:
 
     class VdwInfo;
+    class ParticleTypeInfo;
+    class TypePairInfo;
     NonbondedMethod nonbondedMethod;
+    PotentialFunction potentialFunction;
     double cutoff;
-    bool useDispersionCorrection;
+    bool useDispersionCorrection, useTypes;
     AlchemicalMethod alchemicalMethod;
     int n;
     double alpha;
@@ -334,7 +478,8 @@ private:
 
     std::vector< std::vector<int> > exclusions;
     std::vector<VdwInfo> parameters;
-    std::vector< std::vector< std::vector<double> > > sigEpsTable;
+    std::vector<ParticleTypeInfo> types;
+    std::vector<TypePairInfo> pairs;
 };
 
 /**
@@ -343,7 +488,7 @@ private:
  */
 class AmoebaVdwForce::VdwInfo {
 public:
-    int parentIndex;
+    int parentIndex, typeIndex;
     double reductionFactor, sigma, epsilon, cutoff;
     bool isAlchemical;
     VdwInfo() {
@@ -353,8 +498,36 @@ public:
         epsilon              = 0.0;
         isAlchemical         = false;
     }
-    VdwInfo(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical) :
-        parentIndex(parentIndex), reductionFactor(reductionFactor), sigma(sigma), epsilon(epsilon), isAlchemical(isAlchemical)  {
+    VdwInfo(int parentIndex, double sigma, double epsilon, int typeIndex, double reductionFactor, bool isAlchemical) :
+        parentIndex(parentIndex), reductionFactor(reductionFactor), sigma(sigma), epsilon(epsilon), typeIndex(typeIndex), isAlchemical(isAlchemical) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about a particle type.
+ * @private
+ */
+class AmoebaVdwForce::ParticleTypeInfo {
+public:
+    double sigma, epsilon;
+    ParticleTypeInfo() : sigma(1.0), epsilon(0.0) {
+    }
+    ParticleTypeInfo(double sigma, double epsilon) : sigma(sigma), epsilon(epsilon) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about a type pair.
+ * @private
+ */
+class AmoebaVdwForce::TypePairInfo {
+public:
+    int type1, type2;
+    double sigma, epsilon;
+    TypePairInfo() : type1(-1), type2(-1), sigma(1.0), epsilon(0.0) {
+    }
+    TypePairInfo(int type1, int type2, double sigma, double epsilon) :
+        type1(type1), type2(type2), sigma(sigma), epsilon(epsilon) {
     }
 };
 

plugins/amoeba/openmmapi/include/openmm/internal/AmoebaVdwForceImpl.h

@@ -65,6 +65,18 @@ public:
        return parameters;
     }
     std::vector<std::string> getKernelNames();
+    /**
+     * Compute the matrix of sigma and epsilon to use between every pair of particle types.
+     * 
+     * @param force               the force for which to calculate it
+     * @param[out] type           on exit, this contains the type index of every particle
+     * @param[out] sigmaMatrix    on exit, sigma[i][j] contains the value to use for interactions between particles
+     *                            of types i and j
+     * @param[out] epsilonMatrix  on exit, epsilon[i][j] contains the value to use for interactions between particles
+     *                            of types i and j
+     */
+    static void createParameterMatrix(const AmoebaVdwForce& force, std::vector<int>& type,
+        std::vector<std::vector<double> >& sigmaMatrix, std::vector<std::vector<double> >& epsilonMatrix);
     /**
      * Compute the coefficient which, when divided by the periodic box volume, gives the
      * long range dispersion correction to the energy.

plugins/amoeba/openmmapi/src/AmoebaVdwForce.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:                                                              *
  *                                                                            *
@@ -33,35 +33,90 @@
 #include "openmm/OpenMMException.h"
 #include "openmm/AmoebaVdwForce.h"
 #include "openmm/internal/AmoebaVdwForceImpl.h"
+#include "openmm/internal/AssertionUtilities.h"
 
 using namespace OpenMM;
 using std::string;
 using std::vector;
 
-AmoebaVdwForce::AmoebaVdwForce() : nonbondedMethod(NoCutoff), sigmaCombiningRule("CUBIC-MEAN"), epsilonCombiningRule("HHG"), cutoff(1.0e+10), useDispersionCorrection(true), alchemicalMethod(None), n(5), alpha(0.7) {
+AmoebaVdwForce::AmoebaVdwForce() : nonbondedMethod(NoCutoff), potentialFunction(Buffered147),
+        sigmaCombiningRule("CUBIC-MEAN"), epsilonCombiningRule("HHG"), cutoff(1.0e+10), useDispersionCorrection(true),
+        useTypes(false), alchemicalMethod(None), n(5), alpha(0.7) {
 }
 
 int AmoebaVdwForce::addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical) {
-    parameters.push_back(VdwInfo(parentIndex, sigma, epsilon, reductionFactor, isAlchemical));
+    if (useTypes)
+        throw OpenMMException("AmoebaVdwForce: must use the same version of addParticle() for all particles");
+    parameters.push_back(VdwInfo(parentIndex, sigma, epsilon, -1, reductionFactor, isAlchemical));
+    return parameters.size()-1;
+}
+
+int AmoebaVdwForce::addParticle(int parentIndex, int typeIndex, double reductionFactor, bool isAlchemical) {
+    if (parameters.size() > 0 && !useTypes)
+        throw OpenMMException("AmoebaVdwForce: must use the same version of addParticle() for all particles");
+    useTypes = true;
+    parameters.push_back(VdwInfo(parentIndex, 1.0, 0.0, typeIndex, reductionFactor, isAlchemical));
     return parameters.size()-1;
 }
 
 void AmoebaVdwForce::getParticleParameters(int particleIndex, int& parentIndex,
-                                           double& sigma, double& epsilon, double& reductionFactor, bool& isAlchemical) const {
+                                           double& sigma, double& epsilon, double& reductionFactor, bool& isAlchemical, int& typeIndex) const {
+    ASSERT_VALID_INDEX(particleIndex, parameters);
     parentIndex     = parameters[particleIndex].parentIndex;
     sigma           = parameters[particleIndex].sigma;
     epsilon         = parameters[particleIndex].epsilon;
     reductionFactor = parameters[particleIndex].reductionFactor;
     isAlchemical    = parameters[particleIndex].isAlchemical;
+    typeIndex       = parameters[particleIndex].typeIndex;
 }
 
 void AmoebaVdwForce::setParticleParameters(int particleIndex, int parentIndex,
-                                           double sigma, double epsilon, double reductionFactor, bool isAlchemical) {
+                                           double sigma, double epsilon, double reductionFactor, bool isAlchemical, int typeIndex) {
+    ASSERT_VALID_INDEX(particleIndex, parameters);
     parameters[particleIndex].parentIndex     = parentIndex;
     parameters[particleIndex].sigma           = sigma;
     parameters[particleIndex].epsilon         = epsilon;
     parameters[particleIndex].reductionFactor = reductionFactor;
     parameters[particleIndex].isAlchemical    = isAlchemical;
+    parameters[particleIndex].typeIndex       = typeIndex;
+}
+
+int AmoebaVdwForce::addParticleType(double sigma, double epsilon) {
+    types.push_back(ParticleTypeInfo(sigma, epsilon));
+    return types.size()-1;
+}
+
+void AmoebaVdwForce::getParticleTypeParameters(int typeIndex, double& sigma, double& epsilon) const {
+    ASSERT_VALID_INDEX(typeIndex, types);
+    sigma = types[typeIndex].sigma;
+    epsilon = types[typeIndex].epsilon;
+}
+
+void AmoebaVdwForce::setParticleTypeParameters(int typeIndex, double sigma, double epsilon) {
+    ASSERT_VALID_INDEX(typeIndex, types);
+    types[typeIndex].sigma = sigma;
+    types[typeIndex].epsilon = epsilon;
+}
+
+int AmoebaVdwForce::addTypePair(int type1, int type2, double sigma, double epsilon) {
+    pairs.push_back(TypePairInfo(type1, type2, sigma, epsilon));
+    return pairs.size()-1;
+}
+
+void AmoebaVdwForce::getTypePairParameters(int pairIndex, int& type1, int& type2, double& sigma, double& epsilon) const {
+    ASSERT_VALID_INDEX(pairIndex, pairs);
+    type1 = pairs[pairIndex].type1;
+    type2 = pairs[pairIndex].type2;
+    sigma = pairs[pairIndex].sigma;
+    epsilon = pairs[pairIndex].epsilon;
+}
+
+void AmoebaVdwForce::setTypePairParameters(int pairIndex, int type1, int type2, double sigma, double epsilon) {
+    ASSERT_VALID_INDEX(pairIndex, pairs);
+    pairs[pairIndex].type1 = type1;
+    pairs[pairIndex].type2 = type2;
+    pairs[pairIndex].sigma = sigma;
+    pairs[pairIndex].epsilon = epsilon;
 }
 
 void AmoebaVdwForce::setSigmaCombiningRule(const std::string& inputSigmaCombiningRule) {

plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp

@@ -41,10 +41,6 @@
 using namespace OpenMM;
 using namespace std;
 
-using std::pair;
-using std::vector;
-using std::set;
-
 AmoebaVdwForceImpl::AmoebaVdwForceImpl(const AmoebaVdwForce& owner) : owner(owner) {
 }
 
@@ -77,6 +73,139 @@ double AmoebaVdwForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includ
     return 0.0;
 }
 
+void AmoebaVdwForceImpl::createParameterMatrix(const AmoebaVdwForce& force, vector<int>& type,
+        vector<vector<double> >& sigmaMatrix, vector<vector<double> >& epsilonMatrix) {
+    int numParticles = force.getNumParticles();
+    type.resize(numParticles);
+    int numTypes;
+    vector<double> typeSigma, typeEpsilon;
+    if (force.getUseParticleTypes()) {
+        // We get the types directly from the particles.
+
+        double sigma, epsilon, reduction;
+        int parent;
+        bool isAlchemical;
+        for (int i = 0; i < numParticles; i++)
+            force.getParticleParameters(i, parent, sigma, epsilon, reduction, isAlchemical, type[i]);
+        numTypes = force.getNumParticleTypes();
+        typeSigma.resize(numTypes);
+        typeEpsilon.resize(numTypes);
+        for (int i = 0; i < numTypes; i++)
+            force.getParticleTypeParameters(i, typeSigma[i], typeEpsilon[i]);
+    }
+    else {
+        // Identify types by finding every unique sigma/epsilon pair.
+
+        map<pair<double, double>, int> typeForParams;
+        for (int i = 0; i < numParticles; i++) {
+            double sigma, epsilon, reduction;
+            int parent, typeIndex;
+            bool isAlchemical;
+            force.getParticleParameters(i, parent, sigma, epsilon, reduction, isAlchemical, typeIndex);
+            pair<double, double> params = make_pair(sigma, epsilon);
+            map<pair<double, double>, int>::iterator entry = typeForParams.find(params);
+            if (entry == typeForParams.end()) {
+                int index = typeForParams.size();
+                typeForParams[params] = index;
+            }
+            type[i] = typeForParams[params];
+        }
+        numTypes = typeForParams.size();
+        typeSigma.resize(numTypes);
+        typeEpsilon.resize(numTypes);
+        for (auto params : typeForParams) {
+            typeSigma[params.second] = params.first.first;
+            typeEpsilon[params.second] = params.first.second;
+        }
+    }
+    
+    // Build the matrices by applying combining rules.
+
+    sigmaMatrix.clear();
+    epsilonMatrix.clear();
+    sigmaMatrix.resize(numTypes, vector<double>(numTypes));
+    epsilonMatrix.resize(numTypes, vector<double>(numTypes));
+    string sigmaCombiningRule = force.getSigmaCombiningRule();
+    string epsilonCombiningRule = force.getEpsilonCombiningRule();
+    for (int i = 0; i < numTypes; i++) {
+        double iSigma = typeSigma[i];
+        double iEpsilon = typeEpsilon[i];
+        for (int j = 0; j < numTypes; j++) {
+            double jSigma = typeSigma[j];
+            double jEpsilon = typeEpsilon[j];
+            double sigma, epsilon;
+            // ARITHMETIC = 1
+            // GEOMETRIC  = 2
+            // CUBIC-MEAN = 3
+            if (sigmaCombiningRule == "ARITHMETIC")
+              sigma = iSigma+jSigma;
+            else if (sigmaCombiningRule == "GEOMETRIC")
+              sigma = 2*sqrt(iSigma*jSigma);
+            else if (sigmaCombiningRule == "CUBIC-MEAN") {
+              double iSigma2 = iSigma*iSigma;
+              double jSigma2 = jSigma*jSigma;
+              if ((iSigma2+jSigma2) != 0.0)
+                sigma = 2*(iSigma2*iSigma + jSigma2*jSigma) / (iSigma2+jSigma2);
+              else
+                sigma = 0.0;
+            }
+            else
+                throw OpenMMException("AmoebaVdwForce: Unknown value for sigma combining rule: "+sigmaCombiningRule);
+            sigmaMatrix[i][j] = sigma;
+            sigmaMatrix[j][i] = sigma;
+
+            // ARITHMETIC = 1
+            // GEOMETRIC  = 2
+            // HARMONIC   = 3
+            // W-H        = 4
+            // HHG        = 5
+            if (epsilonCombiningRule == "ARITHMETIC")
+              epsilon = 0.5*(iEpsilon+jEpsilon);
+            else if (epsilonCombiningRule == "GEOMETRIC")
+              epsilon = sqrt(iEpsilon*jEpsilon);
+            else if (epsilonCombiningRule == "HARMONIC") {
+              if ((iEpsilon+jEpsilon) != 0.0)
+                epsilon = 2*(iEpsilon*jEpsilon) / (iEpsilon+jEpsilon);
+              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 = sqrt(iEpsilon*jEpsilon);
+              epsilon = (eps_s == 0.0 ? 0.0 : 2*eps_s*iSigma3*jSigma3/(iSigma6+jSigma6));
+            }
+            else if (epsilonCombiningRule == "HHG") {
+              double epsilonS = sqrt(iEpsilon)+sqrt(jEpsilon);
+              if (epsilonS != 0.0)
+                epsilon = 4*(iEpsilon*jEpsilon) / (epsilonS*epsilonS);
+              else
+                epsilon = 0.0;
+            }
+            else
+                throw OpenMMException("AmoebaVdwForce: Unknown value for epsilon combining rule: "+epsilonCombiningRule);
+            epsilonMatrix[i][j] = epsilon;
+            epsilonMatrix[j][i] = epsilon;
+        }
+    }
+    
+    // Record any type pairs that override the combining rules.
+
+    if (force.getUseParticleTypes()) {
+        for (int i = 0; i < force.getNumTypePairs(); i++) {
+            int type1, type2;
+            double sigma, epsilon;
+            force.getTypePairParameters(i, type1, type2, sigma, epsilon);
+            sigmaMatrix[type1][type2] = sigma;
+            sigmaMatrix[type2][type1] = sigma;
+            epsilonMatrix[type1][type2] = epsilon;
+            epsilonMatrix[type2][type1] = epsilon;
+        }
+    }
+}
+
 double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const AmoebaVdwForce& force) {
 
     // Amoeba VdW dispersion correction implemented by LPW
@@ -84,24 +213,17 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
     if (force.getNonbondedMethod() == AmoebaVdwForce::NoCutoff)
         return 0.0;
 
-    // Identify all particle classes (defined by sigma and epsilon and reduction), and count the number of
+    // Identify all particle classes (defined by sigma and epsilon), and count the number of
     // particles in each class.
 
-    map<pair<double, double>, int> classCounts;
-    for (int i = 0; i < force.getNumParticles(); i++) {
-        double sigma, epsilon, reduction;
-        bool isAlchemical;
-        // The variables reduction, ivindex and isAlchemical are not used.
-        int ivindex;
-        // Get the sigma and epsilon parameters, ignoring everything else.
-        force.getParticleParameters(i, ivindex, sigma, epsilon, reduction, isAlchemical);
-        pair<double, double> key = make_pair(sigma, epsilon);
-        map<pair<double, double>, int>::iterator entry = classCounts.find(key);
-        if (entry == classCounts.end())
-            classCounts[key] = 1;
-        else
-            entry->second++;
-    }
+    vector<int> type;
+    vector<vector<double> > sigmaMatrix;
+    vector<vector<double> > epsilonMatrix;
+    createParameterMatrix(force, type, sigmaMatrix, epsilonMatrix);
+    int numTypes = sigmaMatrix.size();
+    vector<int> typeCounts(numTypes, 0);
+    for (int i = 0; i < force.getNumParticles(); i++)
+        typeCounts[type[i]]++;
 
     // Compute the VdW tapering coefficients.  Mostly copied from amoebaCudaGpu.cpp.
     double cutoff = force.getCutoffDistance();
@@ -143,7 +265,7 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
     c4 = 15.0 * (vdwCut + vdwTaperCut) * denom;
     c5 = -6.0 * denom;
 
-    // Loop over all pairs of classes to compute the coefficient.
+    // Loop over all pairs of types to compute the coefficient.
     // Copied over from TINKER - numerical integration.
     double range = 20.0;
     double cut = vdwTaperCut; // This is where tapering BEGINS
@@ -159,67 +281,14 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
 
     double elrc = 0.0; // This number is incremented and passed out at the end
     double e = 0.0;
-    double sigma, epsilon; // The pairwise sigma and epsilon parameters.
-    int i = 0, k = 0; // Loop counters.
 
     // Double loop over different atom types.
-    std::string sigmaCombiningRule = force.getSigmaCombiningRule();
-    std::string epsilonCombiningRule = force.getEpsilonCombiningRule();
-    for (auto& class1 : classCounts) {
-        k = 0;
-        for (auto& class2 : classCounts) { 
-            // AMOEBA combining rules, copied over from the CUDA code.
-            double iSigma = class1.first.first;
-            double jSigma = class2.first.first;
-            double iEpsilon = class1.first.second;
-            double jEpsilon = class2.first.second;
-            // ARITHMETIC = 1
-            // GEOMETRIC  = 2
-            // CUBIC-MEAN = 3
-            if (sigmaCombiningRule == "ARITHMETIC") {
-              sigma = iSigma + jSigma;
-            } else if (sigmaCombiningRule == "GEOMETRIC") {
-              sigma = 2.0f * std::sqrt(iSigma * jSigma);
-            } else {
-              double iSigma2 = iSigma*iSigma;
-              double jSigma2 = jSigma*jSigma;
-              if ((iSigma2 + jSigma2) != 0.0) {
-                sigma = 2.0f * (iSigma2 * iSigma + jSigma2 * jSigma) / (iSigma2 + jSigma2);
-              } else {
-                sigma = 0.0;
-              }
-            }
-            // ARITHMETIC = 1
-            // GEOMETRIC  = 2
-            // HARMONIC   = 3
-            // W-H        = 4
-            // HHG        = 5
-            if (epsilonCombiningRule == "ARITHMETIC") {
-              epsilon = 0.5f * (iEpsilon + jEpsilon);
-            } else if (epsilonCombiningRule == "GEOMETRIC") {
-              epsilon = std::sqrt(iEpsilon * jEpsilon);
-            } else if (epsilonCombiningRule == "HARMONIC") {
-              if ((iEpsilon + jEpsilon) != 0.0) {
-                epsilon = 2.0f * (iEpsilon * jEpsilon) / (iEpsilon + jEpsilon);
-              } 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) {
-                epsilon = 4.0f * (iEpsilon * jEpsilon) / (epsilonS * epsilonS);
-              } else {
-                epsilon = 0.0;
-              }
-            }
-            int count = class1.second * class2.second;
+    
+    for (int i = 0; i < numTypes; i++) {
+        for (int j = 0; j < numTypes; j++) {
+            double sigma = sigmaMatrix[i][j];
+            double epsilon = epsilonMatrix[i][j];
+            int count = typeCounts[i]*typeCounts[j];
             // Below is an exact copy of stuff from the previous block.
             double rv = sigma;
             double termik = 2.0 * M_PI * count; // termik is equivalent to 2 * pi * count.
@@ -233,18 +302,18 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
                 r2 = r*r;
                 double r3 = r2 * r;
                 double r6 = r3 * r3;
+                if (force.getPotentialFunction() == AmoebaVdwForce::LennardJones) {
+                    double p6 = rv6 / r6;
+                    double p12 = p6 * p6;
+                    e = 4 * epsilon * (p12 - p6);
+                }
+                else {
                     double r7 = r6 * r;
-                // The following is for buffered 14-7 only.
-                /*
-                double rho = r/rv;
-                double term1 = pow(((dhal + 1.0) / (dhal + rho)),7);
-                double term2 = ((ghal + 1.0) / (ghal + pow(rho,7))) - 2.0;
-                e = epsilon * term1 * term2;
-                */
-                double rho = r7 + ghal*rv7;
-                double tau = (dhal + 1.0) / (r + dhal * rv);
+                    double rho = r7 + ghal * rv7;
+                    double tau = (dhal+1.0) / (r+dhal*rv);
                     double tau7 = pow(tau, 7);
-                e = epsilon * rv7 * tau7 * ((ghal + 1.0) * rv7 / rho - 2.0);
+                    e = epsilon * rv7 * tau7 * ((ghal+1.0)*rv7/rho-2.0);
+                }
                 double taper = 0.0;
                 if (r < off) {
                     double r4 = r2 * r2;
@@ -255,9 +324,7 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
                 etot = etot + e * rdelta * r2;
             }
             elrc = elrc + termik * etot;
-            k++;
         }
-        i++;
     }
     return elrc;
 }

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

@@ -2341,12 +2341,12 @@ public:
     ForceInfo(const AmoebaVdwForce& force) : force(force) {
     }
     bool areParticlesIdentical(int particle1, int particle2) {
-        int iv1, iv2;
+        int iv1, iv2, type1, type2;
         double sigma1, sigma2, epsilon1, epsilon2, reduction1, reduction2;
         bool isAlchemical1, isAlchemical2;
-        force.getParticleParameters(particle1, iv1, sigma1, epsilon1, reduction1, isAlchemical1);
-        force.getParticleParameters(particle2, iv2, sigma2, epsilon2, reduction2, isAlchemical2);
-        return (sigma1 == sigma2 && epsilon1 == epsilon2 && reduction1 == reduction2 && isAlchemical1 == isAlchemical2);
+        force.getParticleParameters(particle1, iv1, sigma1, epsilon1, reduction1, isAlchemical1, type1);
+        force.getParticleParameters(particle2, iv2, sigma2, epsilon2, reduction2, isAlchemical2, type2);
+        return (sigma1 == sigma2 && epsilon1 == epsilon2 && reduction1 == reduction2 && isAlchemical1 == isAlchemical2 && type1 == type2);
     }
 private:
     const AmoebaVdwForce& force;
@@ -2366,40 +2366,51 @@ CudaCalcAmoebaVdwForceKernel::~CudaCalcAmoebaVdwForceKernel() {
 
 void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) {
     cu.setAsCurrent();
-    sigmaEpsilon.initialize<float2>(cu, cu.getPaddedNumAtoms(), "sigmaEpsilon");
-    bondReductionAtoms.initialize<int>(cu, cu.getPaddedNumAtoms(), "bondReductionAtoms");
-    bondReductionFactors.initialize<float>(cu, cu.getPaddedNumAtoms(), "bondReductionFactors");
-    tempPosq.initialize(cu, cu.getPaddedNumAtoms(), cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4), "tempPosq");
-    tempForces.initialize<long long>(cu, 3*cu.getPaddedNumAtoms(), "tempForces");
+    int paddedNumAtoms = cu.getPaddedNumAtoms();
+    bondReductionAtoms.initialize<int>(cu, paddedNumAtoms, "bondReductionAtoms");
+    bondReductionFactors.initialize<float>(cu, paddedNumAtoms, "bondReductionFactors");
+    tempPosq.initialize(cu, paddedNumAtoms, cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4), "tempPosq");
+    tempForces.initialize<long long>(cu, 3*paddedNumAtoms, "tempForces");
     
     // Record atom parameters.
-    vector<float2> sigmaEpsilonVec(cu.getPaddedNumAtoms(), make_float2(0, 1));
-    vector<float> isAlchemicalVec(cu.getPaddedNumAtoms(), 0);
-    vector<int> bondReductionAtomsVec(cu.getPaddedNumAtoms(), 0);
-    vector<float> bondReductionFactorsVec(cu.getPaddedNumAtoms(), 0);
+    vector<int> atomTypeVec;
+    vector<vector<double> > sigmaMatrix, epsilonMatrix;
+    AmoebaVdwForceImpl::createParameterMatrix(force, atomTypeVec, sigmaMatrix, epsilonMatrix);
+    atomTypeVec.resize(paddedNumAtoms, 0);
+    int numTypes = sigmaMatrix.size();
+    atomType.initialize<int>(cu, paddedNumAtoms, "atomType");
+    sigmaEpsilon.initialize<float2>(cu, numTypes*numTypes, "sigmaEpsilon");
+    vector<float2> sigmaEpsilonVec(sigmaEpsilon.getSize());
+    for (int i = 0; i < numTypes; i++)
+        for (int j = 0; j < numTypes; j++)
+            sigmaEpsilonVec[i*numTypes+j] = make_float2((float) sigmaMatrix[i][j], (float) epsilonMatrix[i][j]);
+    atomType.upload(atomTypeVec);
+    sigmaEpsilon.upload(sigmaEpsilonVec);
+    
+    vector<float> isAlchemicalVec(paddedNumAtoms, 0);
+    vector<int> bondReductionAtomsVec(paddedNumAtoms, 0);
+    vector<float> bondReductionFactorsVec(paddedNumAtoms, 0);
     vector<vector<int> > exclusions(cu.getNumAtoms());
 
     // Handle Alchemical parameters.
     hasAlchemical = force.getAlchemicalMethod() != AmoebaVdwForce::None;
     if (hasAlchemical) {
-       isAlchemical.initialize<float>(cu, cu.getPaddedNumAtoms(), "isAlchemical");
+       isAlchemical.initialize<float>(cu, paddedNumAtoms, "isAlchemical");
        vdwLambda.initialize<float>(cu, 1, "vdwLambda");
        CHECK_RESULT(cuMemHostAlloc(&vdwLambdaPinnedBuffer, sizeof(float), 0), "Error allocating vdwLambda pinned buffer");
     }
 
     for (int i = 0; i < force.getNumParticles(); i++) {
-        int ivIndex;
+        int ivIndex, type;
         double sigma, epsilon, reductionFactor;
         bool alchemical;
-        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical);
-        sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
+        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical, type);
         isAlchemicalVec[i] = (alchemical) ? 1.0f : 0.0f;
         bondReductionAtomsVec[i] = ivIndex;
         bondReductionFactorsVec[i] = (float) reductionFactor;
         force.getParticleExclusions(i, exclusions[i]);
         exclusions[i].push_back(i);
     }
-    sigmaEpsilon.upload(sigmaEpsilonVec);
     bondReductionAtoms.upload(bondReductionAtomsVec);
     bondReductionFactors.upload(bondReductionFactorsVec);
     if (force.getUseDispersionCorrection())
@@ -2412,7 +2423,8 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
     // this force, so it will have its own neighbor list and interaction kernel.
     
     nonbonded = new CudaNonbondedUtilities(cu);
-    nonbonded->addParameter(CudaNonbondedUtilities::ParameterInfo("sigmaEpsilon", "float", 2, sizeof(float2), sigmaEpsilon.getDevicePointer()));
+    nonbonded->addParameter(CudaNonbondedUtilities::ParameterInfo("atomType", "int", 1, sizeof(int), atomType.getDevicePointer()));
+    nonbonded->addArgument(CudaNonbondedUtilities::ParameterInfo("sigmaEpsilon", "float", 2, sizeof(float2), sigmaEpsilon.getDevicePointer()));
 
     if (hasAlchemical) {
        isAlchemical.upload(isAlchemicalVec);
@@ -2426,32 +2438,11 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
     // Create the interaction kernel.
     
     map<string, string> replacements;
-    string sigmaCombiningRule = force.getSigmaCombiningRule();
-    if (sigmaCombiningRule == "ARITHMETIC")
-        replacements["SIGMA_COMBINING_RULE"] = "1";
-    else if (sigmaCombiningRule == "GEOMETRIC")
-        replacements["SIGMA_COMBINING_RULE"] = "2";
-    else if (sigmaCombiningRule == "CUBIC-MEAN")
-        replacements["SIGMA_COMBINING_RULE"] = "3";
-    else
-        throw OpenMMException("Illegal combining rule for sigma: "+sigmaCombiningRule);
-    string epsilonCombiningRule = force.getEpsilonCombiningRule();
-    if (epsilonCombiningRule == "ARITHMETIC")
-        replacements["EPSILON_COMBINING_RULE"] = "1";
-    else if (epsilonCombiningRule == "GEOMETRIC")
-        replacements["EPSILON_COMBINING_RULE"] = "2";
-    else if (epsilonCombiningRule =="HARMONIC")
-        replacements["EPSILON_COMBINING_RULE"] = "3";
-    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);
-
     replacements["VDW_ALCHEMICAL_METHOD"] = cu.intToString(force.getAlchemicalMethod()); 
     replacements["VDW_SOFTCORE_POWER"] = cu.intToString(force.getSoftcorePower());
     replacements["VDW_SOFTCORE_ALPHA"] = cu.doubleToString(force.getSoftcoreAlpha()); 
+    replacements["POTENTIAL_FUNCTION"] = cu.intToString(force.getPotentialFunction());
+    replacements["NUM_TYPES"] = cu.intToString(numTypes);
 
     double cutoff = force.getCutoffDistance();
     double taperCutoff = cutoff*0.9;
@@ -2467,7 +2458,7 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
     // Create the other kernels.
     
     map<string, string> defines;
-    defines["PADDED_NUM_ATOMS"] = cu.intToString(cu.getPaddedNumAtoms());
+    defines["PADDED_NUM_ATOMS"] = cu.intToString(paddedNumAtoms);
     CUmodule module = cu.createModule(CudaAmoebaKernelSources::amoebaVdwForce1, defines);
     prepareKernel = cu.getKernel(module, "prepareToComputeForce");
     spreadKernel = cu.getKernel(module, "spreadForces");
@@ -2512,22 +2503,33 @@ void CudaCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& context,
     if (force.getNumParticles() != cu.getNumAtoms())
         throw OpenMMException("updateParametersInContext: The number of particles has changed");
     
+    vector<int> atomTypeVec;
+    vector<vector<double> > sigmaMatrix, epsilonMatrix;
+    AmoebaVdwForceImpl::createParameterMatrix(force, atomTypeVec, sigmaMatrix, epsilonMatrix);
+    atomTypeVec.resize(cu.getPaddedNumAtoms(), 0);
+    int numTypes = sigmaMatrix.size();
+    if (sigmaEpsilon.getSize() != numTypes*numTypes)
+        throw OpenMMException("updateParametersInContext: The number of particle types has changed");
+    vector<float2> sigmaEpsilonVec(sigmaEpsilon.getSize());
+    for (int i = 0; i < numTypes; i++)
+        for (int j = 0; j < numTypes; j++)
+            sigmaEpsilonVec[i*numTypes+j] = make_float2((float) sigmaMatrix[i][j], (float) epsilonMatrix[i][j]);
+    atomType.upload(atomTypeVec);
+    sigmaEpsilon.upload(sigmaEpsilonVec);
+
     // Record the per-particle parameters.
-    vector<float2> sigmaEpsilonVec(cu.getPaddedNumAtoms(), make_float2(0, 1));
     vector<float> isAlchemicalVec(cu.getPaddedNumAtoms(), 0);
     vector<int> bondReductionAtomsVec(cu.getPaddedNumAtoms(), 0);
     vector<float> bondReductionFactorsVec(cu.getPaddedNumAtoms(), 0);
     for (int i = 0; i < force.getNumParticles(); i++) {
-        int ivIndex;
+        int ivIndex, type;
         double sigma, epsilon, reductionFactor;
         bool alchemical;
-        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical);
-        sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
+        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical, type);
         isAlchemicalVec[i] = (alchemical) ? 1.0f : 0.0f;
         bondReductionAtomsVec[i] = ivIndex;
         bondReductionFactorsVec[i] = (float) reductionFactor;
     }
-    sigmaEpsilon.upload(sigmaEpsilonVec);
     if (hasAlchemical) isAlchemical.upload(isAlchemicalVec);
     bondReductionAtoms.upload(bondReductionAtomsVec);
     bondReductionFactors.upload(bondReductionFactorsVec);

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.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: Mark Friedrichs, Peter Eastman                                    *
  * Contributors:                                                              *
  *                                                                            *
@@ -584,7 +584,7 @@ private:
     CudaArray isAlchemical;
 
     double dispersionCoefficient;
-    CudaArray sigmaEpsilon;
+    CudaArray sigmaEpsilon, atomType;
     CudaArray bondReductionAtoms;
     CudaArray bondReductionFactors;
     CudaArray tempPosq;

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

@@ -5,34 +5,9 @@
     unsigned int includeInteraction = (!isExcluded);
 #endif
     real tempForce = 0.0f;
-#if SIGMA_COMBINING_RULE == 1
-    real sigma = sigmaEpsilon1.x + sigmaEpsilon2.x;
-#elif SIGMA_COMBINING_RULE == 2
-    real sigma = 2*SQRT(sigmaEpsilon1.x*sigmaEpsilon2.x);
-#else
-    real sigma1_2 = sigmaEpsilon1.x*sigmaEpsilon1.x;
-    real sigma2_2 = sigmaEpsilon2.x*sigmaEpsilon2.x;
-    real sigmasum = sigma1_2+sigma2_2;
-    real sigma = (sigmasum == 0.0f ? (real) 0 : 2*(sigmaEpsilon1.x*sigma1_2 + sigmaEpsilon2.x*sigma2_2)/(sigma1_2+sigma2_2));
-#endif
-#if EPSILON_COMBINING_RULE == 1
-    real epsilon = 0.5f*(sigmaEpsilon1.y + sigmaEpsilon2.y);
-#elif EPSILON_COMBINING_RULE == 2
-    real epsilon = SQRT(sigmaEpsilon1.y*sigmaEpsilon2.y);
-#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));
-#endif
+    float2 pairSigmaEpsilon = sigmaEpsilon[atomType1+atomType2*NUM_TYPES];
+    real sigma = pairSigmaEpsilon.x;
+    real epsilon = pairSigmaEpsilon.y;
     real softcore = 0.0f;
 #if VDW_ALCHEMICAL_METHOD == 1
     if (isAlchemical1 != isAlchemical2) { 
@@ -45,6 +20,15 @@
        softcore = VDW_SOFTCORE_ALPHA * (1.0f - lambda) * (1.0f - lambda);
     }
 #endif
+#if POTENTIAL_FUNCTION == 1
+    real pp1 = sigma / r;
+    real pp2 = pp1 * pp1;
+    real pp3 = pp2 * pp1;
+    real pp6 = pp3 * pp3;
+    real pp12 = pp6 * pp6;
+    real termEnergy = 4 * epsilon * (pp12 - pp6);
+    real deltaE = -24 * epsilon * (2*pp12 - pp6) / r;
+#else
     real dhal = 0.07f;
     real ghal = 0.12f;
     real dhal1 = 1.07f;
@@ -65,6 +49,7 @@
     real dt2 = -7.0f * rho6 * t2 * s2;
     real termEnergy = epsilon * t1 * t2min;
     real deltaE = epsilon * (dt1 * t2min + t1 * dt2) / sigma;
+#endif
 #ifdef USE_CUTOFF
     if (r > TAPER_CUTOFF) {
         real x = r-TAPER_CUTOFF;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaVdwForce.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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,6 +35,7 @@
 
 #include "openmm/internal/AssertionUtilities.h"
 #include "openmm/Context.h"
+#include "openmm/CustomNonbondedForce.h"
 #include "OpenMMAmoeba.h"
 #include "openmm/System.h"
 #include "openmm/AmoebaVdwForce.h"
@@ -138,13 +139,13 @@ void testVdw() {
         positions[ii][2] *= AngstromToNm;
     }
     for (int ii = 0; ii < amoebaVdwForce->getNumParticles();  ii++) {
-        int indexIV;
+        int indexIV, type;
         double sigma, epsilon, reduction;
         bool isAlchemical;
-        amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
+        amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical, type);
         sigma        *= AngstromToNm;
         epsilon      *= CalToJoule;
-        amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
+        amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical, type);
     }
     platformName = "CUDA";
     Context context(system, integrator, Platform::getPlatformByName(platformName));
@@ -169,11 +170,11 @@ void testVdw() {
     // Try changing the particle parameters and make sure it's still correct.
     
     for (int i = 0; i < numberOfParticles; i++) {
-        int indexIV;
+        int indexIV, type;
         double mass, sigma, epsilon, reduction;
         bool isAlchemical;
-        amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction, isAlchemical);
-        amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction, isAlchemical);
+        amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction, isAlchemical, type);
+        amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction, isAlchemical, type);
     }
     LangevinIntegrator integrator2(0.0, 0.1, 0.01);
     Context context2(system, integrator2, Platform::getPlatformByName(platformName));
@@ -2106,6 +2107,94 @@ void testTriclinic() {
     }
 }
 
+void testCompareToCustom(AmoebaVdwForce::PotentialFunction potential) {
+    const int numParticles = 10;
+    const double cutoff = 1.0;
+    System system;
+    AmoebaVdwForce* vdw = new AmoebaVdwForce();
+    vdw->setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
+    vdw->setCutoff(cutoff);
+    vdw->setPotentialFunction(potential);
+    vdw->setSigmaCombiningRule("ARITHMETIC");
+    vdw->setEpsilonCombiningRule("GEOMETRIC");
+    vdw->setUseDispersionCorrection(true);
+    system.addForce(vdw);
+    CustomNonbondedForce* custom;
+    if (potential == AmoebaVdwForce::LennardJones)
+        custom = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6); sigma=sigma1+sigma2; eps=sqrt(eps1*eps2)");
+    else
+        custom = new CustomNonbondedForce("eps*((1.07/(p+0.07))^7 * ((1.12/(p^7+0.12))-2)); p=r/sigma; sigma=sigma1+sigma2; eps=sqrt(eps1*eps2)");
+    custom->addPerParticleParameter("sigma");
+    custom->addPerParticleParameter("eps");
+    custom->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
+    custom->setCutoffDistance(cutoff);
+    custom->setUseLongRangeCorrection(true);
+    custom->setUseSwitchingFunction(true);
+    custom->setSwitchingDistance(0.9*cutoff);
+    custom->setForceGroup(1);
+    system.addForce(custom);
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        if (i%2 == 0) {
+            vdw->addParticle(i, 0.2, 1.0, 0.0);
+            custom->addParticle({0.2, 1.0});
+        }
+        else {
+            vdw->addParticle(i, 0.25, 0.5, 0.0);
+            custom->addParticle({0.25, 0.5});
+        }
+        positions.push_back(2*Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+    }
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+    context.setPositions(positions);
+    State state1 = context.getState(State::Energy | State::Forces, true, 1);
+    State state2 = context.getState(State::Energy | State::Forces, true, 2);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-3);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
+}
+
+void testParticleTypes() {
+    System system;
+    for (int i = 0; i < 4; i++)
+        system.addParticle(1.0);
+    AmoebaVdwForce* vdw = new AmoebaVdwForce();
+    system.addForce(vdw);
+    vdw->setPotentialFunction(AmoebaVdwForce::LennardJones);
+    vdw->setSigmaCombiningRule("ARITHMETIC");
+    vdw->setEpsilonCombiningRule("GEOMETRIC");
+    vdw->addParticle(0, 0, 1.0);
+    vdw->addParticle(1, 2, 1.0);
+    vdw->addParticle(2, 0, 1.0);
+    vdw->addParticle(3, 1, 1.0);
+    vdw->addParticleType(0.3, 1.0);
+    vdw->addParticleType(0.4, 1.1);
+    vdw->addParticleType(0.5, 1.2);
+    vdw->addTypePair(2, 0, 0.6, 1.5);
+    vector<Vec3> positions;
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(1, 0, 0));
+    positions.push_back(Vec3(0, 1, 0));
+    positions.push_back(Vec3(1, 1, 0));
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    vector<double> r = {1.0, 1.0, sqrt(2.0), sqrt(2.0), 1.0, 1.0};
+    vector<double> sigma = {0.6, 0.3+0.3, 0.3+0.4, 0.6, 0.5+0.4, 0.3+0.4};
+    vector<double> epsilon = {1.5, sqrt(1.0*1.0), sqrt(1.0*1.1), 1.5, sqrt(1.1*1.2), sqrt(1.0*1.1)};
+    double expectedEnergy = 0;
+    for (int i = 0; i < 6; i++) {
+        double p = sigma[i]/r[i];
+        expectedEnergy += 4*epsilon[i]*(pow(p, 12) - pow(p, 6));
+    }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         std::cout << "TestCudaAmoebaVdwForce running test..." << std::endl;
@@ -2162,6 +2251,15 @@ int main(int argc, char* argv[]) {
 
         testTriclinic();
         
+        // Compare to an equivalent CustomNonbondedForce.
+
+        testCompareToCustom(AmoebaVdwForce::Buffered147);
+        testCompareToCustom(AmoebaVdwForce::LennardJones);
+        
+        // Test specifying parameters by particle type.
+        
+        testParticleTypes();
+
         // Set lambda and the softcore power (n) to any values, while softcore alpha set to 0. 
         // The energy and forces are equal to scaling those from testVdwAmmoniaCubicMeanHhg by lambda^n;
         int n = 5;

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:                                                              *
  *                                                                            *
@@ -32,7 +32,6 @@
 #include "AmoebaReferenceStretchBendForce.h"
 #include "AmoebaReferenceOutOfPlaneBendForce.h"
 #include "AmoebaReferenceTorsionTorsionForce.h"
-#include "AmoebaReferenceVdwForce.h"
 #include "AmoebaReferenceWcaDispersionForce.h"
 #include "AmoebaReferenceGeneralizedKirkwoodForce.h"
 #include "openmm/internal/AmoebaTorsionTorsionForceImpl.h"
@@ -1001,9 +1000,8 @@ ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(const std::
 }
 
 ReferenceCalcAmoebaVdwForceKernel::~ReferenceCalcAmoebaVdwForceKernel() {
-    if (neighborList) {
+    if (neighborList)
         delete neighborList;
-    } 
 }
 
 void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) {
@@ -1011,102 +1009,41 @@ 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;
-        double sigma, epsilon, reduction;
-        bool alchemical;
-        std::vector<int> exclusions;
-
-        force.getParticleParameters(ii, indexIV, sigma, epsilon, reduction, alchemical);
-        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);
 }
 
 double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
 
     vector<Vec3>& posData   = extractPositions(context);
     vector<Vec3>& forceData = extractForces(context);
-    AmoebaReferenceVdwForce vdwForce(sigmaCombiningRule, epsilonCombiningRule);
-    if (alchemicalMethod == AmoebaVdwForce::Decouple) {
-       vdwForce.setAlchemicalMethod(AmoebaReferenceVdwForce::Decouple);
-    } else if (alchemicalMethod == AmoebaVdwForce::Annihilate) {
-       vdwForce.setAlchemicalMethod(AmoebaReferenceVdwForce::Annihilate);
-    } else {
-       vdwForce.setAlchemicalMethod(AmoebaReferenceVdwForce::None);
-    }
-    vdwForce.setSoftcorePower(n);
-    vdwForce.setSoftcoreAlpha(alpha);
     double energy;
     double lambda = context.getParameter(AmoebaVdwForce::Lambda());
     if (useCutoff) {
-        vdwForce.setCutoff(cutoff);
-        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) {
-            vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::CutoffPeriodic);
             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 {
-            vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::CutoffNonPeriodic);
         }
-    } else {
-        vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::NoCutoff);
-        energy = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, indexIVs, sigmas, epsilons, 
-                                                  reductions, isAlchemical, allExclusions, forceData);
     }
+    else
+        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;
-        double sigma, epsilon, reduction;
-        bool alchemical;
-        force.getParticleParameters(i, indexIV, sigma, epsilon, reduction, alchemical);
-        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:                                                              *
  *                                                                            *
@@ -33,6 +33,7 @@
 #include "openmm/HippoNonbondedForce.h"
 #include "AmoebaReferenceMultipoleForce.h"
 #include "AmoebaReferenceHippoNonbondedForce.h"
+#include "AmoebaReferenceVdwForce.h"
 #include "ReferenceNeighborList.h"
 #include "SimTKOpenMMRealType.h"
 
@@ -501,17 +502,7 @@ private:
     int usePBC;
     double cutoff;
     double dispersionCoefficient;
-    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;
+    AmoebaReferenceVdwForce vdwForce;
     const System& system;
     NeighborList* neighborList;
 };

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

@@ -25,6 +25,7 @@
 #include "AmoebaReferenceForce.h"
 #include "AmoebaReferenceVdwForce.h"
 #include "ReferenceForce.h"
+#include "openmm/internal/AmoebaVdwForceImpl.h"
 #include <algorithm>
 #include <cctype>
 #include <cmath>
@@ -32,55 +33,38 @@
 using std::vector;
 using namespace OpenMM;
 
-AmoebaReferenceVdwForce::AmoebaReferenceVdwForce() : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9), _n(5), _alpha(0.7), _alchemicalMethod(None) {
-
-    setTaperCoefficients(_cutoff);
-    setSigmaCombiningRule("ARITHMETIC");
-    setEpsilonCombiningRule("GEOMETRIC");
-}
-
-
-AmoebaReferenceVdwForce::AmoebaReferenceVdwForce(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule) : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9),  _n(5), _alpha(0.7), _alchemicalMethod(None) {
-
+AmoebaReferenceVdwForce::AmoebaReferenceVdwForce() : _nonbondedMethod(AmoebaVdwForce::NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9), _n(5), _alpha(0.7), _alchemicalMethod(AmoebaVdwForce::None) {
     setTaperCoefficients(_cutoff);
-    setSigmaCombiningRule(sigmaCombiningRule);
-    setEpsilonCombiningRule(epsilonCombiningRule);
-}
-
-AmoebaReferenceVdwForce::NonbondedMethod AmoebaReferenceVdwForce::getNonbondedMethod() const {
-    return _nonbondedMethod;
-}
-
-void AmoebaReferenceVdwForce::setNonbondedMethod(AmoebaReferenceVdwForce::NonbondedMethod nonbondedMethod) {
-    _nonbondedMethod = nonbondedMethod;
-}
-
-AmoebaReferenceVdwForce::AlchemicalMethod AmoebaReferenceVdwForce::getAlchemicalMethod() const {
-    return _alchemicalMethod;
-}
-
-
-void AmoebaReferenceVdwForce::setAlchemicalMethod(AmoebaReferenceVdwForce::AlchemicalMethod alchemicalMethod){
-    _alchemicalMethod = alchemicalMethod;
-}
-
-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::initialize(const AmoebaVdwForce& force) {
+    _alchemicalMethod = force.getAlchemicalMethod();
+    _n = force.getSoftcorePower();
+    _alpha = force.getSoftcoreAlpha();
+    _nonbondedMethod = force.getNonbondedMethod();
+    potentialFunction = force.getPotentialFunction();
+    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]);
+    }
 }
 
-
 void AmoebaReferenceVdwForce::setTaperCoefficients(double cutoff) {
     _taperCutoff = cutoff*_taperCutoffFactor;
     if (_taperCutoff != cutoff) {
@@ -99,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,
@@ -219,14 +119,23 @@ double AmoebaReferenceVdwForce::calculatePairIxn(double combinedSigma, double co
     // get deltaR, R2, and R between 2 atoms
 
     double deltaR[ReferenceForce::LastDeltaRIndex];
-    if (_nonbondedMethod == CutoffPeriodic)
+    if (_nonbondedMethod == AmoebaVdwForce::CutoffPeriodic)
         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;
+    if (potentialFunction == AmoebaVdwForce::LennardJones) {
+        double pp1 = combinedSigma / r_ij;
+        double pp2 = pp1 * pp1;
+        double pp3 = pp2 * pp1;
+        double pp6 = pp3 * pp3;
+        double pp12 = pp6 * pp6;
+        energy = 4 * combinedEpsilon * (pp12 - pp6);
+        dEdR = -24 * combinedEpsilon * (2*pp12 - pp6) / r_ij;
+    }
+    else {
         double rho = r_ij / combinedSigma;
         double rho2 = rho * rho;
         double rho6 = rho2 * rho2 * rho2;
@@ -241,11 +150,12 @@ double AmoebaReferenceVdwForce::calculatePairIxn(double combinedSigma, double co
         double t2min = t2 - 2;
         double dt1 = -7.0 * rhodec * t1 * s1;
         double dt2 = -7.0 * rho6 * t2 * s2;
-    double energy = combinedEpsilon * t1 * t2min;
-    double dEdR = combinedEpsilon * (dt1 * t2min + t1 * dt2) / combinedSigma;
+        energy = combinedEpsilon * t1 * t2min;
+        dEdR = combinedEpsilon * (dt1 * t2min + t1 * dt2) / combinedSigma;
+    }
 
     // tapering
-    if ((_nonbondedMethod == CutoffNonPeriodic || _nonbondedMethod == CutoffPeriodic) && r_ij > _taperCutoff) {
+    if (_nonbondedMethod == AmoebaVdwForce::CutoffPeriodic && r_ij > _taperCutoff) {
         double delta    = r_ij - _taperCutoff;
         double taper    = 1.0 + delta*delta*delta*(_taperCoefficients[C3] + delta*(_taperCoefficients[C4] + delta*_taperCoefficients[C5]));
         double dtaper   = delta*delta*(3.0*_taperCoefficients[C3] + delta*(4.0*_taperCoefficients[C4] + delta*5.0*_taperCoefficients[C5]));
@@ -276,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
@@ -310,31 +214,24 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
     std::vector<unsigned int> exclusions(numParticles, 0);
     for (unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++) {
 
-        double sigmaI      = sigmas[ii];
-        double epsilonI    = epsilons[ii];
         bool isAlchemicalI = isAlchemical[ii];
         for (int jj : allExclusions[ii])
             exclusions[jj] = 1;
 
         for (unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++) {
             if (exclusions[jj] == 0) {
-
-                double combinedSigma   = (this->*_combineSigmas)(sigmaI, sigmas[jj]);
-
-                double combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj], sigmaI, sigmas[jj]);
-
-              
+                double combinedSigma = sigmaMatrix[particleType[ii]][particleType[jj]];
+                double combinedEpsilon = epsilonMatrix[particleType[ii]][particleType[jj]];
                 double softcore = 0.0;
 
-                if (this->_alchemicalMethod == Decouple && (isAlchemicalI != isAlchemical[jj])) {
+                if (this->_alchemicalMethod == AmoebaVdwForce::Decouple && (isAlchemicalI != isAlchemical[jj])) {
                    combinedEpsilon *= pow(lambda, this->_n);
                    softcore = this->_alpha * pow(1.0 - lambda, 2);
-                } else if (this->_alchemicalMethod == Annihilate && (isAlchemicalI || isAlchemical[jj])) {
+                } else if (this->_alchemicalMethod == AmoebaVdwForce::Annihilate && (isAlchemicalI || isAlchemical[jj])) {
                    combinedEpsilon *= pow(lambda, this->_n);
                    softcore = this->_alpha * pow(1.0 - lambda, 2);
                 }
 
-
                 Vec3 force;
                 energy += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
                                            reducedPositions[ii], reducedPositions[jj], force);
@@ -366,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 {
 
@@ -392,23 +284,21 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
         int siteI                   = pair.first;
         int siteJ                   = pair.second;
 
-        double combinedSigma   = (this->*_combineSigmas)(sigmas[siteI], sigmas[siteJ]);
-
-        double combinedEpsilon = (this->*_combineEpsilons)(epsilons[siteI], epsilons[siteJ], sigmas[siteI], sigmas[siteJ]);
+        double combinedSigma = sigmaMatrix[particleType[siteI]][particleType[siteJ]];
+        double combinedEpsilon = epsilonMatrix[particleType[siteI]][particleType[siteJ]];
 
         double softcore        = 0.0;
         int isAlchemicalI      = isAlchemical[siteI];
         int isAlchemicalJ      = isAlchemical[siteJ];
 
-        if (this->_alchemicalMethod == Decouple && (isAlchemicalI != isAlchemicalJ)) {
+        if (this->_alchemicalMethod == AmoebaVdwForce::Decouple && (isAlchemicalI != isAlchemicalJ)) {
            combinedEpsilon *= pow(lambda, this->_n);
            softcore = this->_alpha * pow(1.0 - lambda, 2);
-        } else if (this->_alchemicalMethod == Annihilate && (isAlchemicalI || isAlchemicalJ)) {
+        } else if (this->_alchemicalMethod == AmoebaVdwForce::Annihilate && (isAlchemicalI || isAlchemicalJ)) {
            combinedEpsilon *= pow(lambda, this->_n);
            softcore = this->_alpha * pow(1.0 - lambda, 2);
         }
 
-
         Vec3 force;
         energy += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
                                    reducedPositions[siteI], reducedPositions[siteJ], force);
@@ -427,7 +317,6 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
         } else {
             addReducedForce(siteJ, indexIVs[siteJ], reductions[siteJ], 1.0, force, forces);
         }
-
     }
 
     return energy;

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
@@ -26,7 +26,9 @@
 #define __AmoebaReferenceVdwForce_H__
 
 #include "openmm/Vec3.h"
+#include "openmm/AmoebaVdwForce.h"
 #include "ReferenceNeighborList.h"
+#include <set>
 #include <string>
 #include <vector>
 
@@ -42,48 +44,6 @@ class AmoebaReferenceVdwForce {
 
 public:
  
-    /** 
-     * This is an enumeration of the different methods that may be used for handling long range Vdw forces.
-     */
-    enum NonbondedMethod {
-
-        /**
-         * No cutoff is applied to the interactions.  The full set of N^2 interactions is computed exactly.
-         * This necessarily means that periodic boundary conditions cannot be used.  This is the default.
-         */
-
-        NoCutoff = 0,
-
-        /**
-         * Interactions beyond the cutoff distance are ignored.  
-         */
-        CutoffNonPeriodic = 1,
-        /**
-         * Periodic boundary conditions are used, so that each particle interacts only with the nearest periodic copy of
-         * each other particle.  Interactions beyond the cutoff distance are ignored.  
-         */
-        CutoffPeriodic = 2,
-    };
-
-    /**
-     * This is an enumeration of the different alchemical methods used when applying softcore interactions.
-     */
-    enum AlchemicalMethod {
-        /**
-         * All vdW interactions are treated normally. This is the default.
-         */
-        None = 0,
-        /**
-         * Maintain full strength vdW interactions between two alchemical particles.
-         */
-        Decouple = 1,
-        /**
-         * Interactions between two alchemical particles are turned off at lambda=0.
-         */
-        Annihilate = 2,
-    };
-
- 
     /**---------------------------------------------------------------------------------------
        
        Constructor
@@ -92,77 +52,11 @@ public:
  
     AmoebaReferenceVdwForce();
     
-    /**---------------------------------------------------------------------------------------
-       
-       Constructor
-       
-       --------------------------------------------------------------------------------------- */
- 
-    AmoebaReferenceVdwForce(const std::string& sigmaCombiningRule,
-                            const std::string& epsilonCombiningRule);
+    void initialize(const AmoebaVdwForce& force);
  
     /**---------------------------------------------------------------------------------------
     
-       Destructor
-       
-       --------------------------------------------------------------------------------------- */
- 
-    ~AmoebaReferenceVdwForce() {};
- 
-    /**---------------------------------------------------------------------------------------
-    
-       Get nonbonded method
-    
-       @return nonbonded method
-    
-       --------------------------------------------------------------------------------------- */
-    
-    NonbondedMethod getNonbondedMethod() const;
-
-    /**---------------------------------------------------------------------------------------
-    
-       Set nonbonded method
-    
-       @param nonbonded method
-    
-       --------------------------------------------------------------------------------------- */
-    
-    void setNonbondedMethod(NonbondedMethod nonbondedMethod);
-
-    /**---------------------------------------------------------------------------------------
-
-       Get alchemical method
-
-       @return alchemical method
-
-       --------------------------------------------------------------------------------------- */
-
-    AlchemicalMethod getAlchemicalMethod() const;
-
-    /**---------------------------------------------------------------------------------------
-
-       Set alchemical method
-
-       @param alchemical method
-
-       --------------------------------------------------------------------------------------- */
-
-    void setAlchemicalMethod(AlchemicalMethod alchemicalMethod);
-
-
-    /**---------------------------------------------------------------------------------------
-    
-       Get cutoff
-    
-       @return cutoff
-    
-       --------------------------------------------------------------------------------------- */
-    
-    double getCutoff() const;
-
-    /**---------------------------------------------------------------------------------------
-    
-       Set cutof
+       Set cutoff
     
        @param cutoff
     
@@ -172,94 +66,21 @@ public:
 
     /**---------------------------------------------------------------------------------------
     
-       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();
 
     /**---------------------------------------------------------------------------------------
     
@@ -268,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
@@ -281,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;
          
     /**---------------------------------------------------------------------------------------
@@ -295,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
     
@@ -308,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
@@ -321,29 +121,23 @@ private:
     static const int C4=1;
     static const int C5=2;
 
-    std::string _sigmaCombiningRule;
-    std::string _epsilonCombiningRule;
-    NonbondedMethod _nonbondedMethod;
-    AlchemicalMethod _alchemicalMethod;
+    AmoebaVdwForce::NonbondedMethod _nonbondedMethod;
+    AmoebaVdwForce::AlchemicalMethod _alchemicalMethod;
+    AmoebaVdwForce::PotentialFunction potentialFunction;
     double _n;
     double _alpha;
     double _cutoff;
     double _taperCutoffFactor;
     double _taperCutoff;
     double _taperCoefficients[3];
+    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;
-
 
     /**---------------------------------------------------------------------------------------
     

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

@@ -35,6 +35,7 @@
 
 #include "openmm/internal/AssertionUtilities.h"
 #include "openmm/Context.h"
+#include "openmm/CustomNonbondedForce.h"
 #include "OpenMMAmoeba.h"
 #include "openmm/System.h"
 #include "openmm/AmoebaVdwForce.h"
@@ -140,13 +141,13 @@ void testVdw() {
         positions[ii][2] *= AngstromToNm;
     }
     for (int ii = 0; ii < amoebaVdwForce->getNumParticles();  ii++) {
-        int indexIV;
+        int indexIV, type;
         double sigma, epsilon, reduction;
         bool isAlchemical;
-        amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
+        amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical, type);
         sigma        *= AngstromToNm;
         epsilon      *= CalToJoule;
-        amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
+        amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical, type);
     }
     platformName = "Reference";
     Context context(system, integrator, Platform::getPlatformByName(platformName));
@@ -172,11 +173,11 @@ void testVdw() {
     // Try changing the particle parameters and make sure it's still correct.
     
     for (int i = 0; i < numberOfParticles; i++) {
-        int indexIV;
+        int indexIV, type;
         double mass, sigma, epsilon, reduction;
         bool isAlchemical;
-        amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction, isAlchemical);
-        amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction, isAlchemical);
+        amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction, isAlchemical, type);
+        amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction, isAlchemical, type);
     }
     LangevinIntegrator integrator2(0.0, 0.1, 0.01);
     Context context2(system, integrator2, Platform::getPlatformByName(platformName));
@@ -2116,6 +2117,94 @@ void testTriclinic() {
     }
 }
 
+void testCompareToCustom(AmoebaVdwForce::PotentialFunction potential) {
+    const int numParticles = 10;
+    const double cutoff = 1.0;
+    System system;
+    AmoebaVdwForce* vdw = new AmoebaVdwForce();
+    vdw->setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
+    vdw->setCutoff(cutoff);
+    vdw->setPotentialFunction(potential);
+    vdw->setSigmaCombiningRule("ARITHMETIC");
+    vdw->setEpsilonCombiningRule("GEOMETRIC");
+    vdw->setUseDispersionCorrection(true);
+    system.addForce(vdw);
+    CustomNonbondedForce* custom;
+    if (potential == AmoebaVdwForce::LennardJones)
+        custom = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6); sigma=sigma1+sigma2; eps=sqrt(eps1*eps2)");
+    else
+        custom = new CustomNonbondedForce("eps*((1.07/(p+0.07))^7 * ((1.12/(p^7+0.12))-2)); p=r/sigma; sigma=sigma1+sigma2; eps=sqrt(eps1*eps2)");
+    custom->addPerParticleParameter("sigma");
+    custom->addPerParticleParameter("eps");
+    custom->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
+    custom->setCutoffDistance(cutoff);
+    custom->setUseLongRangeCorrection(true);
+    custom->setUseSwitchingFunction(true);
+    custom->setSwitchingDistance(0.9*cutoff);
+    custom->setForceGroup(1);
+    system.addForce(custom);
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        if (i%2 == 0) {
+            vdw->addParticle(i, 0.2, 1.0, 0.0);
+            custom->addParticle({0.2, 1.0});
+        }
+        else {
+            vdw->addParticle(i, 0.25, 0.5, 0.0);
+            custom->addParticle({0.25, 0.5});
+        }
+        positions.push_back(2*Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+    }
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    context.setPositions(positions);
+    State state1 = context.getState(State::Energy | State::Forces, true, 1);
+    State state2 = context.getState(State::Energy | State::Forces, true, 2);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-3);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
+}
+
+void testParticleTypes() {
+    System system;
+    for (int i = 0; i < 4; i++)
+        system.addParticle(1.0);
+    AmoebaVdwForce* vdw = new AmoebaVdwForce();
+    system.addForce(vdw);
+    vdw->setPotentialFunction(AmoebaVdwForce::LennardJones);
+    vdw->setSigmaCombiningRule("ARITHMETIC");
+    vdw->setEpsilonCombiningRule("GEOMETRIC");
+    vdw->addParticle(0, 0, 1.0);
+    vdw->addParticle(1, 2, 1.0);
+    vdw->addParticle(2, 0, 1.0);
+    vdw->addParticle(3, 1, 1.0);
+    vdw->addParticleType(0.3, 1.0);
+    vdw->addParticleType(0.4, 1.1);
+    vdw->addParticleType(0.5, 1.2);
+    vdw->addTypePair(2, 0, 0.6, 1.5);
+    vector<Vec3> positions;
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(1, 0, 0));
+    positions.push_back(Vec3(0, 1, 0));
+    positions.push_back(Vec3(1, 1, 0));
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    vector<double> r = {1.0, 1.0, sqrt(2.0), sqrt(2.0), 1.0, 1.0};
+    vector<double> sigma = {0.6, 0.3+0.3, 0.3+0.4, 0.6, 0.5+0.4, 0.3+0.4};
+    vector<double> epsilon = {1.5, sqrt(1.0*1.0), sqrt(1.0*1.1), 1.5, sqrt(1.1*1.2), sqrt(1.0*1.1)};
+    double expectedEnergy = 0;
+    for (int i = 0; i < 6; i++) {
+        double p = sigma[i]/r[i];
+        expectedEnergy += 4*epsilon[i]*(pow(p, 12) - pow(p, 6));
+    }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
+}
+
 int main(int numberOfArguments, char* argv[]) {
 
     try {
@@ -2170,6 +2259,15 @@ int main(int numberOfArguments, char* argv[]) {
         
         testTriclinic();
         
+        // Compare to an equivalent CustomNonbondedForce.
+
+        testCompareToCustom(AmoebaVdwForce::Buffered147);
+        testCompareToCustom(AmoebaVdwForce::LennardJones);
+        
+        // Test specifying parameters by particle type.
+        
+        testParticleTypes();
+
         // Set lambda and the softcore power (n) to any values (softcore alpha set to 0). 
         // The energy and forces are equal to scaling testVdwAmmoniaCubicMeanHhg by lambda^n;
         int n = 5;
@@ -2188,7 +2286,6 @@ int main(int numberOfArguments, char* argv[]) {
         lambda = 0.0;
         alpha = 0.7;
         testVdwAlchemical(n, alpha, lambda, method);
-    
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/amoeba/serialization/src/AmoebaVdwForceProxy.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) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -42,8 +42,9 @@ AmoebaVdwForceProxy::AmoebaVdwForceProxy() : SerializationProxy("AmoebaVdwForce"
 }
 
 void AmoebaVdwForceProxy::serialize(const void* object, SerializationNode& node) const {
-    node.setIntProperty("version", 3);
+    node.setIntProperty("version", 4);
     const AmoebaVdwForce& force = *reinterpret_cast<const AmoebaVdwForce*>(object);
+    bool useTypes = force.getUseParticleTypes();
 
     node.setIntProperty("forceGroup", force.getForceGroup());
     node.setStringProperty("SigmaCombiningRule", force.getSigmaCombiningRule());
@@ -54,31 +55,49 @@ void AmoebaVdwForceProxy::serialize(const void* object, SerializationNode& node)
     node.setDoubleProperty("n", force.getSoftcorePower());
     node.setDoubleProperty("alpha", force.getSoftcoreAlpha());
     node.setIntProperty("alchemicalMethod", (int) force.getAlchemicalMethod());
+    node.setIntProperty("potentialFunction", (int) force.getPotentialFunction());
+    node.setBoolProperty("useTypes", useTypes);
 
     SerializationNode& particles = node.createChildNode("VdwParticles");
-    for (unsigned int ii = 0; ii < static_cast<unsigned int>(force.getNumParticles()); ii++) {
-
-        int ivIndex;
+    for (int i = 0; i < force.getNumParticles(); i++) {
+        int ivIndex, typeIndex;
         double sigma, epsilon, reductionFactor;
         bool isAlchemical;
-        force.getParticleParameters(ii, ivIndex, sigma, epsilon, reductionFactor, isAlchemical);
-
+        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, isAlchemical, typeIndex);
         SerializationNode& particle = particles.createChildNode("Particle");
-        particle.setIntProperty("ivIndex", ivIndex).setDoubleProperty("sigma", sigma).setDoubleProperty("epsilon", epsilon).setDoubleProperty("reductionFactor", reductionFactor).setBoolProperty("isAlchemical",isAlchemical);
+        if (useTypes)
+            particle.setIntProperty("ivIndex", ivIndex).setIntProperty("type", typeIndex).setDoubleProperty("reductionFactor", reductionFactor).setBoolProperty("isAlchemical", isAlchemical);
+        else
+            particle.setIntProperty("ivIndex", ivIndex).setDoubleProperty("sigma", sigma).setDoubleProperty("epsilon", epsilon).setDoubleProperty("reductionFactor", reductionFactor).setBoolProperty("isAlchemical", isAlchemical); 
 
         std::vector< int > exclusions;
-        force.getParticleExclusions(ii,  exclusions);
-
+        force.getParticleExclusions(i,  exclusions);
         SerializationNode& particleExclusions = particle.createChildNode("ParticleExclusions");
-        for (unsigned int jj = 0; jj < exclusions.size(); jj++) {
-            particleExclusions.createChildNode("excl").setIntProperty("index", exclusions[jj]);
+        for (int j = 0; j < exclusions.size(); j++)
+            particleExclusions.createChildNode("excl").setIntProperty("index", exclusions[j]);
+    }
+    if (useTypes) {
+        SerializationNode& types = node.createChildNode("ParticleTypes");
+        for (int i = 0; i < force.getNumParticleTypes(); i++) {
+            double sigma, epsilon;
+            force.getParticleTypeParameters(i, sigma, epsilon);
+            SerializationNode& type = types.createChildNode("Type");
+            type.setDoubleProperty("sigma", sigma).setDoubleProperty("epsilon", epsilon);
+        }
+        SerializationNode& pairs = node.createChildNode("TypePairs");
+        for (int i = 0; i < force.getNumParticleTypes(); i++) {
+            int type1, type2;
+            double sigma, epsilon;
+            force.getTypePairParameters(i, type1, type2, sigma, epsilon);
+            SerializationNode& pair = pairs.createChildNode("Pair");
+            pair.setIntProperty("type1", type1).setIntProperty("type2", type2).setDoubleProperty("sigma", sigma).setDoubleProperty("epsilon", epsilon);
         }
     }
 }
 
 void* AmoebaVdwForceProxy::deserialize(const SerializationNode& node) const {
     int version = node.getIntProperty("version");
-    if (version < 1 || version > 3)
+    if (version < 1 || version > 4)
         throw OpenMMException("Unsupported version number");
     AmoebaVdwForce* force = new AmoebaVdwForce();
     try {
@@ -95,13 +114,22 @@ void* AmoebaVdwForceProxy::deserialize(const SerializationNode& node) const {
            force->setSoftcoreAlpha(node.getDoubleProperty("alpha"));
         }
         
+        bool useTypes = false;
+        if (version > 3) {
+            useTypes = node.getBoolProperty("useTypes");
+           force->setPotentialFunction((AmoebaVdwForce::PotentialFunction) node.getIntProperty("potentialFunction"));
+        }
+
         const SerializationNode& particles = node.getChildNode("VdwParticles");
-        for (unsigned int ii = 0; ii < particles.getChildren().size(); ii++) {
-            const SerializationNode& particle = particles.getChildren()[ii];
+        for (int i = 0; i < particles.getChildren().size(); i++) {
+            const SerializationNode& particle = particles.getChildren()[i];
 
             if (version < 3) 
                force->addParticle(particle.getIntProperty("ivIndex"), particle.getDoubleProperty("sigma"), 
                                   particle.getDoubleProperty("epsilon"), particle.getDoubleProperty("reductionFactor"));
+            else if (useTypes)
+               force->addParticle(particle.getIntProperty("ivIndex"), particle.getIntProperty("type"), 
+                                  particle.getDoubleProperty("reductionFactor"), particle.getBoolProperty("isAlchemical"));
             else
                force->addParticle(particle.getIntProperty("ivIndex"), particle.getDoubleProperty("sigma"), 
                                   particle.getDoubleProperty("epsilon"), particle.getDoubleProperty("reductionFactor"), 
@@ -110,13 +138,19 @@ void* AmoebaVdwForceProxy::deserialize(const SerializationNode& node) const {
             // exclusions
 
             const SerializationNode& particleExclusions = particle.getChildNode("ParticleExclusions");
-            std::vector< int > exclusions;
-            for (unsigned int jj = 0; jj < particleExclusions.getChildren().size(); jj++) {
-                exclusions.push_back(particleExclusions.getChildren()[jj].getIntProperty("index"));
+            std::vector<int> exclusions;
+            for (int j = 0; j < particleExclusions.getChildren().size(); j++)
+                exclusions.push_back(particleExclusions.getChildren()[j].getIntProperty("index"));
+            force->setParticleExclusions(i, exclusions);
         }
-            force->setParticleExclusions(ii, exclusions);
+        if (useTypes) {
+            const SerializationNode& types = node.getChildNode("ParticleTypes");
+            for (const auto& type : types.getChildren())
+                force->addParticleType(type.getDoubleProperty("sigma"), type.getDoubleProperty("epsilon"));
+            const SerializationNode& pairs = node.getChildNode("TypePairs");
+            for (const auto& pair : pairs.getChildren())
+                force->addTypePair(pair.getIntProperty("type1"), pair.getIntProperty("type2"), pair.getDoubleProperty("sigma"), pair.getDoubleProperty("epsilon"));
         }
-
     }
     catch (...) {
         delete force;

plugins/amoeba/serialization/tests/TestSerializeAmoebaVdwForce.cpp

@@ -51,16 +51,17 @@ void testSerialization() {
     force1.setCutoff(0.9);
     force1.setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
     force1.setAlchemicalMethod(AmoebaVdwForce::None);
+    force1.setPotentialFunction(AmoebaVdwForce::Buffered147);
 
     force1.addParticle(0, 1.0, 2.0, 0.9, false);
     force1.addParticle(1, 1.1, 2.1, 0.9, true);
     force1.addParticle(2, 1.3, 4.1, 0.9, false);
-    for (unsigned int ii = 0; ii < 3; ii++) {
+    for (int i = 0; i < 3; i++) {
         std::vector< int > exclusions;
-        exclusions.push_back(ii);
-        exclusions.push_back(ii + 1);
-        exclusions.push_back(ii + 10);
-        force1.setParticleExclusions(ii, exclusions);
+        exclusions.push_back(i);
+        exclusions.push_back(i + 1);
+        exclusions.push_back(i + 10);
+        force1.setParticleExclusions(i, exclusions);
     }
 
     // Serialize and then deserialize it.
@@ -78,13 +79,14 @@ void testSerialization() {
     ASSERT_EQUAL(force1.getCutoff(),                force2.getCutoff());
     ASSERT_EQUAL(force1.getNonbondedMethod(),       force2.getNonbondedMethod());
     ASSERT_EQUAL(force1.getAlchemicalMethod(),      force2.getAlchemicalMethod());
+    ASSERT_EQUAL(force1.getPotentialFunction(),     force2.getPotentialFunction());
 
     ASSERT_EQUAL(force1.getNumParticles(),          force2.getNumParticles());
 
-    for (unsigned int ii = 0; ii < static_cast<unsigned int>(force1.getNumParticles()); ii++) {
+    for (int i = 0; i < force1.getNumParticles(); i++) {
 
-        int ivIndex1;
-        int ivIndex2;
+        int ivIndex1, type1;
+        int ivIndex2, type2;
 
         double sigma1, epsilon1, reductionFactor1;
         double sigma2, epsilon2, reductionFactor2;
@@ -92,38 +94,111 @@ void testSerialization() {
         bool isAlchemical1;
         bool isAlchemical2;
 
-        force1.getParticleParameters(ii, ivIndex1, sigma1, epsilon1, reductionFactor1, isAlchemical1);
-        force2.getParticleParameters(ii, ivIndex2, sigma2, epsilon2, reductionFactor2, isAlchemical2);
+        force1.getParticleParameters(i, ivIndex1, sigma1, epsilon1, reductionFactor1, isAlchemical1, type1);
+        force2.getParticleParameters(i, ivIndex2, sigma2, epsilon2, reductionFactor2, isAlchemical2, type2);
 
         ASSERT_EQUAL(ivIndex1,          ivIndex2);
         ASSERT_EQUAL(sigma1,            sigma2);
         ASSERT_EQUAL(epsilon1,          epsilon2);
         ASSERT_EQUAL(reductionFactor1,  reductionFactor2);
         ASSERT_EQUAL(isAlchemical1,     isAlchemical2);
+        ASSERT_EQUAL(type1,             type2);
     }
-    for (unsigned int ii = 0; ii < static_cast<unsigned int>(force1.getNumParticles()); ii++) {
+    for (int i = 0; i < force1.getNumParticles(); i++) {
 
         std::vector< int > exclusions1;
         std::vector< int > exclusions2;
 
-        force1.getParticleExclusions(ii, exclusions1);
-        force2.getParticleExclusions(ii, exclusions2);
+        force1.getParticleExclusions(i, exclusions1);
+        force2.getParticleExclusions(i, exclusions2);
 
         ASSERT_EQUAL(exclusions1.size(), exclusions2.size());
-        for (unsigned int jj = 0; jj < exclusions1.size(); jj++) {
+        for (int j = 0; j < exclusions1.size(); j++) {
             int hit = 0;
-            for (unsigned int kk = 0; kk < exclusions2.size(); kk++) {
-                if (exclusions2[jj] == exclusions1[kk])hit++;
+            for (int kk = 0; kk < exclusions2.size(); kk++) {
+                if (exclusions2[j] == exclusions1[kk])hit++;
             }
             ASSERT_EQUAL(hit, 1);
         }
     }
 }
 
+void testSerializeTypes() {
+    // Create a Force that specifies parameters by type.
+
+    AmoebaVdwForce force1;
+    force1.setPotentialFunction(AmoebaVdwForce::LennardJones);
+
+    force1.addParticle(0, 2, 1.0, false);
+    force1.addParticle(1, 2, 0.9, true);
+    force1.addParticle(2, 0, 1.0, false);
+    force1.addParticle(3, 1, 0.9, false);
+    force1.addParticleType(1.1, 2.0);
+    force1.addParticleType(1.2, 2.1);
+    force1.addParticleType(1.3, 2.2);
+    force1.addTypePair(0, 2, 1.5, 2.5);
+
+    // Serialize and then deserialize it.
+
+    stringstream buffer;
+    XmlSerializer::serialize<AmoebaVdwForce>(&force1, "Force", buffer);
+    AmoebaVdwForce* copy = XmlSerializer::deserialize<AmoebaVdwForce>(buffer);
+
+    // Compare the two forces to see if they are identical.  
+    AmoebaVdwForce& force2 = *copy;
+
+    ASSERT_EQUAL(force1.getPotentialFunction(), force2.getPotentialFunction());
+    ASSERT_EQUAL(force1.getNumParticles(),      force2.getNumParticles());
+    ASSERT_EQUAL(force1.getNumParticleTypes(),  force2.getNumParticleTypes());
+    ASSERT_EQUAL(force1.getNumTypePairs(),      force2.getNumTypePairs());
+
+    for (int i = 0; i < force1.getNumParticles(); i++) {
+        int ivIndex1, type1;
+        int ivIndex2, type2;
+
+        double sigma1, epsilon1, reductionFactor1;
+        double sigma2, epsilon2, reductionFactor2;
+
+        bool isAlchemical1;
+        bool isAlchemical2;
+
+        force1.getParticleParameters(i, ivIndex1, sigma1, epsilon1, reductionFactor1, isAlchemical1, type1);
+        force2.getParticleParameters(i, ivIndex2, sigma2, epsilon2, reductionFactor2, isAlchemical2, type2);
+
+        ASSERT_EQUAL(ivIndex1,          ivIndex2);
+        ASSERT_EQUAL(sigma1,            sigma2);
+        ASSERT_EQUAL(epsilon1,          epsilon2);
+        ASSERT_EQUAL(reductionFactor1,  reductionFactor2);
+        ASSERT_EQUAL(isAlchemical1,     isAlchemical2);
+        ASSERT_EQUAL(type1,             type2);
+    }
+    for (int i = 0; i < force1.getNumParticleTypes(); i++) {
+        double sigma1, epsilon1;
+        double sigma2, epsilon2;
+        force1.getParticleTypeParameters(i, sigma1, epsilon1);
+        force2.getParticleTypeParameters(i, sigma2, epsilon2);
+        ASSERT_EQUAL(sigma1, sigma2);
+        ASSERT_EQUAL(epsilon1, epsilon2);
+    }
+    for (int i = 0; i < force1.getNumTypePairs(); i++) {
+        int type11, type21;
+        int type12, type22;
+        double sigma1, epsilon1;
+        double sigma2, epsilon2;
+        force1.getTypePairParameters(i, type11, type21, sigma1, epsilon1);
+        force2.getTypePairParameters(i, type12, type22, sigma2, epsilon2);
+        ASSERT_EQUAL(type11, type12);
+        ASSERT_EQUAL(type21, type22);
+        ASSERT_EQUAL(sigma1, sigma2);
+        ASSERT_EQUAL(epsilon1, epsilon2);
+    }
+}
+
 int main() {
     try {
         registerAmoebaSerializationProxies();
         testSerialization();
+        testSerializeTypes();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

plugins/amoeba/wrappers/generateAmoebaWrappers.py

@@ -257,6 +257,7 @@ class CHeaderGenerator(WrapperGenerator):
             methodNames[methodNode] = shortMethodDefinition.split()[-1]
         
         # Write other methods
+        nameCount = {}
         for methodNode in methodList:
             methodName = methodNames[methodNode]
             if methodName in (shortClassName, destructorName):
@@ -268,6 +269,13 @@ class CHeaderGenerator(WrapperGenerator):
                 # There are two identical methods that differ only in whether they are const.  Skip the const one.
                 continue
             returnType = self.getType(getText("type", methodNode))
+            if methodName in nameCount:
+                # There are multiple methods with the same name.
+                count = nameCount[methodName]
+                methodName = "%s_%d" % (methodName, count)
+                nameCount[methodName] = count+1
+            else:
+                nameCount[methodName] = 1
             self.out.write("extern OPENMM_EXPORT_AMOEBA %s %s_%s(" % (returnType, typeName, methodName))
             isInstanceMethod = (methodNode.attrib['static'] != 'yes')
             if isInstanceMethod:
@@ -438,6 +446,7 @@ class CSourceGenerator(WrapperGenerator):
             methodNames[methodNode] = shortMethodDefinition.split()[-1]
         
         # Write other methods
+        nameCount = {}
         for methodNode in methodList:
             methodName = methodNames[methodNode]
             if methodName in (shortClassName, destructorName):
@@ -448,6 +457,13 @@ class CSourceGenerator(WrapperGenerator):
             if isConstMethod and any(methodNames[m] == methodName and m.attrib['const'] == 'no' for m in methodList):
                 # There are two identical methods that differ only in whether they are const.  Skip the const one.
                 continue
+            if methodName in nameCount:
+                # There are multiple methods with the same name.
+                count = nameCount[methodName]
+                methodName = "%s_%d" % (methodName, count)
+                nameCount[methodName] = count+1
+            else:
+                nameCount[methodName] = 1
             methodType = getText("type", methodNode)
             returnType = self.getType(methodType)
             if methodType in self.classesByShortName:
@@ -474,7 +490,7 @@ class CSourceGenerator(WrapperGenerator):
                 self.out.write('%s*>(target)->' % className)
             else:
                 self.out.write('%s::' % className)
-            self.out.write('%s(' % methodName)
+            self.out.write('%s(' % methodNames[methodNode])
             self.writeInvocationArguments(methodNode, False)
             self.out.write(');\n')
             if returnType != 'void':
@@ -758,6 +774,7 @@ class FortranHeaderGenerator(WrapperGenerator):
             methodNames[methodNode] = shortMethodDefinition.split()[-1]
         
         # Write other methods
+        nameCount = {}
         for methodNode in methodList:
             methodName = methodNames[methodNode]
             if methodName in (shortClassName, destructorName):
@@ -768,6 +785,13 @@ class FortranHeaderGenerator(WrapperGenerator):
             if isConstMethod and any(methodNames[m] == methodName and m.attrib['const'] == 'no' for m in methodList):
                 # There are two identical methods that differ only in whether they are const.  Skip the const one.
                 continue
+            if methodName in nameCount:
+                # There are multiple methods with the same name.
+                count = nameCount[methodName]
+                methodName = "%s_%d" % (methodName, count)
+                nameCount[methodName] = count+1
+            else:
+                nameCount[methodName] = 1
             returnType = self.getType(getText("type", methodNode))
             hasReturnValue = (returnType in ('integer*4', 'real*8'))
             hasReturnArg = not (hasReturnValue or returnType == 'void')
@@ -973,6 +997,7 @@ class FortranSourceGenerator(WrapperGenerator):
             methodNames[methodNode] = shortMethodDefinition.split()[-1]
         
         # Write other methods
+        nameCount = {}
         for methodNode in methodList:
             methodName = methodNames[methodNode]
             if methodName in (shortClassName, destructorName):
@@ -985,6 +1010,13 @@ class FortranSourceGenerator(WrapperGenerator):
             if isConstMethod and any(methodNames[m] == methodName and m.attrib['const'] == 'no' for m in methodList):
                 # There are two identical methods that differ only in whether they are const.  Skip the const one.
                 continue
+            if methodName in nameCount:
+                # There are multiple methods with the same name.
+                count = nameCount[methodName]
+                methodName = "%s_%d" % (methodName, count)
+                nameCount[methodName] = count+1
+            else:
+                nameCount[methodName] = 1
             functionName = "%s_%s" % (typeName, methodName)
             self.writeOneMethod(classNode, methodNode, functionName, functionName.lower()+'_')
             self.writeOneMethod(classNode, methodNode, functionName, functionName.upper())

wrappers/generateWrappers.py

@@ -267,6 +267,7 @@ class CHeaderGenerator(WrapperGenerator):
             methodNames[methodNode] = shortMethodDefinition.split()[-1]
         
         # Write other methods
+        nameCount = {}
         for methodNode in methodList:
             methodName = methodNames[methodNode]
             if methodName in (shortClassName, destructorName):
@@ -278,6 +279,13 @@ class CHeaderGenerator(WrapperGenerator):
                 # There are two identical methods that differ only in whether they are const.  Skip the const one.
                 continue
             returnType = self.getType(getText("type", methodNode))
+            if methodName in nameCount:
+                # There are multiple methods with the same name.
+                count = nameCount[methodName]
+                methodName = "%s_%d" % (methodName, count)
+                nameCount[methodName] = count+1
+            else:
+                nameCount[methodName] = 1
             self.out.write("extern OPENMM_EXPORT %s %s_%s(" % (returnType, typeName, methodName))
             isInstanceMethod = (methodNode.attrib['static'] != 'yes')
             if isInstanceMethod:
@@ -526,6 +534,7 @@ class CSourceGenerator(WrapperGenerator):
             methodNames[methodNode] = shortMethodDefinition.split()[-1]
         
         # Write other methods
+        nameCount = {}
         for methodNode in methodList:
             methodName = methodNames[methodNode]
             if methodName in (shortClassName, destructorName):
@@ -536,6 +545,13 @@ class CSourceGenerator(WrapperGenerator):
             if isConstMethod and any(methodNames[m] == methodName and m.attrib['const'] == 'no' for m in methodList):
                 # There are two identical methods that differ only in whether they are const.  Skip the const one.
                 continue
+            if methodName in nameCount:
+                # There are multiple methods with the same name.
+                count = nameCount[methodName]
+                methodName = "%s_%d" % (methodName, count)
+                nameCount[methodName] = count+1
+            else:
+                nameCount[methodName] = 1
             methodType = getText("type", methodNode)
             returnType = self.getType(methodType)
             if methodType in self.classesByShortName:
@@ -562,7 +578,7 @@ class CSourceGenerator(WrapperGenerator):
                 self.out.write('%s*>(target)->' % className)
             else:
                 self.out.write('%s::' % className)
-            self.out.write('%s(' % methodName)
+            self.out.write('%s(' % methodNames[methodNode])
             self.writeInvocationArguments(methodNode, False)
             self.out.write(');\n')
             if returnType != 'void':
@@ -982,6 +998,7 @@ class FortranHeaderGenerator(WrapperGenerator):
             methodNames[methodNode] = shortMethodDefinition.split()[-1]
         
         # Write other methods
+        nameCount = {}
         for methodNode in methodList:
             methodName = methodNames[methodNode]
             if methodName in (shortClassName, destructorName):
@@ -992,6 +1009,13 @@ class FortranHeaderGenerator(WrapperGenerator):
             if isConstMethod and any(methodNames[m] == methodName and m.attrib['const'] == 'no' for m in methodList):
                 # There are two identical methods that differ only in whether they are const.  Skip the const one.
                 continue
+            if methodName in nameCount:
+                # There are multiple methods with the same name.
+                count = nameCount[methodName]
+                methodName = "%s_%d" % (methodName, count)
+                nameCount[methodName] = count+1
+            else:
+                nameCount[methodName] = 1
             returnType = self.getType(getText("type", methodNode))
             hasReturnValue = (returnType in ('integer*4', 'real*8'))
             hasReturnArg = not (hasReturnValue or returnType == 'void')
@@ -1539,6 +1563,7 @@ class FortranSourceGenerator(WrapperGenerator):
             methodNames[methodNode] = shortMethodDefinition.split()[-1]
         
         # Write other methods
+        nameCount = {}
         for methodNode in methodList:
             methodName = methodNames[methodNode]
             if methodName in (shortClassName, destructorName):
@@ -1551,6 +1576,13 @@ class FortranSourceGenerator(WrapperGenerator):
             if isConstMethod and any(methodNames[m] == methodName and m.attrib['const'] == 'no' for m in methodList):
                 # There are two identical methods that differ only in whether they are const.  Skip the const one.
                 continue
+            if methodName in nameCount:
+                # There are multiple methods with the same name.
+                count = nameCount[methodName]
+                methodName = "%s_%d" % (methodName, count)
+                nameCount[methodName] = count+1
+            else:
+                nameCount[methodName] = 1
             functionName = "%s_%s" % (typeName, methodName)
             truncatedName = functionName[:63]
             self.writeOneMethod(classNode, methodNode, functionName, truncatedName.lower()+'_')

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

@@ -4537,12 +4537,21 @@ class AmoebaVdwGenerator(object):
 
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
 
+        potentialMap = {'BUFFERED-14-7':0, 'LENNARD-JONES':1}
         sigmaMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'CUBIC-MEAN':1}
         epsilonMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'HARMONIC':1, 'W-H':1, 'HHG':1}
 
         force = mm.AmoebaVdwForce()
         sys.addForce(force)
 
+        # Potential function
+
+        if (self.type.upper() in potentialMap):
+            force.setPotentialFunction(potentialMap[self.type.upper()])
+        else:
+            stringList = ' '.join(str(x) for x in potentialMap.keys())
+            raise ValueError("AmoebaVdwGenerator: potential type %s not recognized; valid values are %s; using default." % (self.type, stringList))
+
         # sigma and epsilon combining rules
 
         if ('sigmaCombiningRule' in args):

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -325,10 +325,13 @@ UNITS = {
 ("AmoebaVdwForce",                        "getEpsilonCombiningRule")                       :  ( None, ()),
 ("AmoebaVdwForce",                        "getParticleExclusions")                         :  ( None, ()),
 ("AmoebaVdwForce",                        "getAlchemicalMethod")                           :  ( None, ()),
+("AmoebaVdwForce",                        "getPotentialFunction")                          :  ( None, ()),
 ("AmoebaVdwForce",                        "getSoftcorePower")                              :  ( None, ()),
 ("AmoebaVdwForce",                        "getSoftcoreAlpha")                              :  ( None, ()),
 ("AmoebaVdwForce",                        "getCutoff")                                     :  ( 'unit.nanometer', ()),
-("AmoebaVdwForce",                        "getParticleParameters")                         :  ( None, (None, 'unit.nanometer', 'unit.kilojoule_per_mole', None, None)),
+("AmoebaVdwForce",                        "getParticleParameters")                         :  ( None, (None, 'unit.nanometer', 'unit.kilojoule_per_mole', None, None, None)),
+("AmoebaVdwForce",                        "getParticleTypeParameters")                     :  ( None, (None, 'unit.nanometer', 'unit.kilojoule_per_mole')),
+("AmoebaVdwForce",                        "getTypePairParameters")                         :  ( None, (None, None, None, 'unit.nanometer', 'unit.kilojoule_per_mole')),
 
 ("AmoebaWcaDispersionForce",              "getParticleParameters")                         :  ( None, ('unit.nanometer', 'unit.kilojoule_per_mole')),
 ("AmoebaWcaDispersionForce",              "getAwater")                                     :  ( '1/(unit.nanometer*unit.nanometer*unit.nanometer)',()),