Continuing to new features in AmoebaVdwForce

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

@@ -34,6 +34,7 @@
 
 #include "openmm/Force.h"
 #include "internal/windowsExportAmoeba.h"
+#include <string>
 #include <vector>
 
 namespace OpenMM {

plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp

@@ -109,6 +109,8 @@ void AmoebaVdwForceImpl::createParameterMatrix(const AmoebaVdwForce& force, vect
             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;
@@ -212,7 +214,7 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
     vector<vector<double> > sigmaMatrix;
     vector<vector<double> > epsilonMatrix;
     createParameterMatrix(force, type, sigmaMatrix, epsilonMatrix);
-    int numTypes = type.size();
+    int numTypes = sigmaMatrix.size();
     vector<int> typeCounts(numTypes, 0);
     for (int i = 0; i < force.getNumParticles(); i++)
         typeCounts[type[i]]++;
@@ -277,14 +279,10 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
     // Double loop over different atom types.
     
     for (int i = 0; i < numTypes; i++) {
-        for (int j = 0; j <= i; j++) {
+        for (int j = 0; j < numTypes; 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];
+            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.
@@ -298,18 +296,18 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
                 r2 = r*r;
                 double r3 = r2 * r;
                 double r6 = r3 * r3;
-                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 tau7 = pow(tau, 7);
-                e = epsilon * rv7 * tau7 * ((ghal + 1.0) * rv7 / rho - 2.0);
+                if (force.getPotentialFunction() == AmoebaVdwForce::LennardJones) {
+                    double p6 = rv6 / r6;
+                    double p12 = p6 * p6;
+                    e = 4 * epsilon * (p12 - p6);
+                }
+                else {
+                    double r7 = r6 * r;
+                    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);
+                }
                 double taper = 0.0;
                 if (r < off) {
                     double r4 = r2 * r2;

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

@@ -2388,10 +2388,10 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
     }
 
     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);
+        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical, type);
         sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
         isAlchemicalVec[i] = (alchemical) ? 1.0f : 0.0f;
         bondReductionAtomsVec[i] = ivIndex;
@@ -2518,10 +2518,10 @@ void CudaCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& context,
     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);
+        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical, type);
         sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
         isAlchemicalVec[i] = (alchemical) ? 1.0f : 0.0f;
         bondReductionAtomsVec[i] = ivIndex;

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

@@ -138,13 +138,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 +169,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));

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp

@@ -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"
@@ -1021,12 +1020,12 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
 
     for (int ii = 0; ii < numParticles; ii++) {
 
-        int indexIV;
+        int indexIV, type;
         double sigma, epsilon, reduction;
         bool alchemical;
         std::vector<int> exclusions;
 
-        force.getParticleParameters(ii, indexIV, sigma, epsilon, reduction, alchemical);
+        force.getParticleParameters(ii, indexIV, sigma, epsilon, reduction, alchemical, type);
         force.getParticleExclusions(ii, exclusions);
         for (unsigned int jj = 0; jj < exclusions.size(); jj++) {
            allExclusions[ii].insert(exclusions[jj]);
@@ -1048,29 +1047,18 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
     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);
         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) {
@@ -1080,14 +1068,11 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
             energy  = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, indexIVs, sigmas, epsilons, 
                                                        reductions, isAlchemical, *neighborList, forceData);
             energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
-        } else {
-            vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::CutoffNonPeriodic);
         }
-    } else {
-        vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::NoCutoff);
+    }
+    else
         energy = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, indexIVs, sigmas, epsilons, 
                                                   reductions, isAlchemical, allExclusions, forceData);
-    }
     return static_cast<double>(energy);
 }
 
@@ -1097,16 +1082,17 @@ void ReferenceCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& con
 
     // Record the values.
     for (int i = 0; i < numParticles; ++i) {
-        int indexIV;
+        int indexIV, type;
         double sigma, epsilon, reduction;
         bool alchemical;
-        force.getParticleParameters(i, indexIV, sigma, epsilon, reduction, alchemical);
+        force.getParticleParameters(i, indexIV, sigma, epsilon, reduction, alchemical, type);
         indexIVs[i] = indexIV;
         sigmas[i] = sigma;
         epsilons[i] = epsilon;
         reductions[i]= reduction;
         isAlchemical[i]= alchemical;
     }
+    vdwForce.initialize(force);
 }
 
 /* -------------------------------------------------------------------------- *

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h

@@ -33,6 +33,7 @@
 #include "openmm/HippoNonbondedForce.h"
 #include "AmoebaReferenceMultipoleForce.h"
 #include "AmoebaReferenceHippoNonbondedForce.h"
+#include "AmoebaReferenceVdwForce.h"
 #include "ReferenceNeighborList.h"
 #include "SimTKOpenMMRealType.h"
 
@@ -501,6 +502,7 @@ private:
     int usePBC;
     double cutoff;
     double dispersionCoefficient;
+    AmoebaReferenceVdwForce vdwForce;
     AmoebaVdwForce::AlchemicalMethod alchemicalMethod;
     int n;
     double alpha;

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,38 +33,31 @@
 using std::vector;
 using namespace OpenMM;
 
-AmoebaReferenceVdwForce::AmoebaReferenceVdwForce() : _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("ARITHMETIC");
     setEpsilonCombiningRule("GEOMETRIC");
 }
 
+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);
+}
 
-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(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule) : _nonbondedMethod(AmoebaVdwForce::NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9),  _n(5), _alpha(0.7), _alchemicalMethod(AmoebaVdwForce::None) {
 
     setTaperCoefficients(_cutoff);
     setSigmaCombiningRule(sigmaCombiningRule);
     setEpsilonCombiningRule(epsilonCombiningRule);
 }
 
-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;
 }
@@ -219,7 +213,7 @@ 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);
@@ -227,25 +221,37 @@ double AmoebaReferenceVdwForce::calculatePairIxn(double combinedSigma, double co
     double r_ij_2       = deltaR[ReferenceForce::R2Index];
     double r_ij         = deltaR[ReferenceForce::RIndex];
 
-    double rho = r_ij / combinedSigma;
-    double rho2 = rho * rho;
-    double rho6 = rho2 * rho2 * rho2;
-    double rhoplus = rho + dhal;
-    double rhodec2 = rhoplus * rhoplus;
-    double rhodec = rhodec2 * rhodec2 * rhodec2;
-    double s1 = 1.0 / (softcore + rhodec * rhoplus);
-    double s2 = 1.0 / (softcore + rho6 * rho + 0.12);
-    double point72 = dhal1 * dhal1;
-    double t1 = dhal1 * point72 * point72 * point72 * s1;
-    double t2 = ghal1 * s2;
-    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;
+    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;
+        double rhoplus = rho + dhal;
+        double rhodec2 = rhoplus * rhoplus;
+        double rhodec = rhodec2 * rhodec2 * rhodec2;
+        double s1 = 1.0 / (softcore + rhodec * rhoplus);
+        double s2 = 1.0 / (softcore + rho6 * rho + 0.12);
+        double point72 = dhal1 * dhal1;
+        double t1 = dhal1 * point72 * point72 * point72 * s1;
+        double t2 = ghal1 * s2;
+        double t2min = t2 - 2;
+        double dt1 = -7.0 * rhodec * t1 * s1;
+        double dt2 = -7.0 * rho6 * t2 * s2;
+        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]));
@@ -309,32 +315,25 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
     double energy = 0.0;
     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);
@@ -392,23 +391,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 +424,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

@@ -26,6 +26,7 @@
 #define __AmoebaReferenceVdwForce_H__
 
 #include "openmm/Vec3.h"
+#include "openmm/AmoebaVdwForce.h"
 #include "ReferenceNeighborList.h"
 #include <string>
 #include <vector>
@@ -41,48 +42,6 @@ typedef double (AmoebaReferenceVdwForce::*CombiningFunctionEpsilon)(double x, do
 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,
-    };
-
  
     /**---------------------------------------------------------------------------------------
        
@@ -91,6 +50,8 @@ public:
        --------------------------------------------------------------------------------------- */
  
     AmoebaReferenceVdwForce();
+    
+    void initialize(const AmoebaVdwForce& force);
  
     /**---------------------------------------------------------------------------------------
        
@@ -109,46 +70,6 @@ public:
  
     ~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);
-
 
     /**---------------------------------------------------------------------------------------
     
@@ -323,14 +244,18 @@ private:
 
     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;
     Vec3 _periodicBoxVectors[3];
     CombiningFunction _combineSigmas;
     double arithmeticSigmaCombiningRule(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,57 @@ 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);
+}
+
 int main(int numberOfArguments, char* argv[]) {
 
     try {
@@ -2169,7 +2221,12 @@ int main(int numberOfArguments, char* argv[]) {
         // test triclinic boxes
         
         testTriclinic();
+        
+        // Compare to an equivalent CustomNonbondedForce.
 
+        testCompareToCustom(AmoebaVdwForce::Buffered147);
+        testCompareToCustom(AmoebaVdwForce::LennardJones);
+        
         // 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;