Began implementing new features in AmoebaVdwForce

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:                                                              *
  *                                                                            *
@@ -39,12 +39,25 @@
 namespace OpenMM {
 
 /**
- * This class implements a buffered 14-7 potential used to model van der Waals forces.
- *
- * 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 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.
  *
+ * 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
  * particle to another one.  This is typically done for hydrogens to place the interaction site slightly
@@ -91,6 +104,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 +148,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 +172,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 +187,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 +206,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
      *
@@ -210,6 +329,13 @@ public:
     void setUseDispersionCorrection(bool useCorrection) {
         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 +394,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 +462,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 +477,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 +487,7 @@ private:
  */
 class AmoebaVdwForce::VdwInfo {
 public:
-    int parentIndex;
+    int parentIndex, typeIndex;
     double reductionFactor, sigma, epsilon, cutoff;
     bool isAlchemical;
     VdwInfo() {
@@ -353,8 +497,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 = pairs[typeIndex].sigma;
+    epsilon = pairs[typeIndex].epsilon;
+}
+
+void AmoebaVdwForce::setParticleTypeParameters(int typeIndex, double sigma, double epsilon) {
+    ASSERT_VALID_INDEX(typeIndex, types);
+    pairs[typeIndex].sigma = sigma;
+    pairs[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,131 @@ 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())
+                typeForParams[params] = typeForParams.size();
+            type[i] = typeForParams[params];
+        }
+        numTypes = typeForParams.size();
+        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 {
+              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;
+            }
+            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 {
+              double epsilonS = sqrt(iEpsilon)+sqrt(jEpsilon);
+              if (epsilonS != 0.0)
+                epsilon = 4*(iEpsilon*jEpsilon) / (epsilonS*epsilonS);
+              else
+                epsilon = 0.0;
+            }
+            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 +205,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 = type.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 +257,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 +273,18 @@ 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 <= i; j++) {
+            double sigma = sigmaMatrix[i][j];
+            double epsilon = epsilonMatrix[i][j];
+            int count;
+            if (i == j)
+                count = typeCounts[i]*(typeCounts[i]+1)/2;
+            else
+                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.
@@ -255,9 +320,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;

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 {
@@ -94,15 +113,24 @@ void* AmoebaVdwForceProxy::deserialize(const SerializationNode& node) const {
            force->setSoftcorePower(node.getDoubleProperty("n"));
            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 
+            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"), 
                                   particle.getBoolProperty("isAlchemical"));
@@ -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"));
-            }
-            force->setParticleExclusions(ii, exclusions);
+            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);
+        }
+        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(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;