Lots of cleanup to extrapolated polarization

plugins/amoeba/openmmapi/include/openmm/AmoebaMultipoleForce.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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs, Peter Eastman                                    *
  * Contributors:                                                              *
  *                                                                            *
@@ -72,19 +72,24 @@ public:
     enum PolarizationType {
 
         /**
-         * Mutual polarization
+         * Full mutually induced polarization.  The dipoles are iterated until the converge to the accuracy specified
+         * by getMutualInducedTargetEpsilon().
          */
         Mutual = 0,
 
         /**
-         * Direct polarization
+         * Direct polarization approximation.  The induced dipoles depend only on the fixed multipoles, not on other
+         * induced dipoles.
          */
         Direct = 1,
 
         /**
-         * Optimized perturbation theory
+         * Extrapolated perturbation theory approximation.  The dipoles are iterated a few times, and then an analytic
+         * approximation is used to extrapolate to the fully converged values.  Call setExtrapolationCoefficients()
+         * to set the coefficients used for the extrapolation.  The default coefficients used in this release are
+         * [0, -0.3, 0, 1.3], but be aware that those may change in a future release.
          */
-        OPT = 2
+        Extrapolated = 2
 
     };
 
@@ -305,19 +310,21 @@ public:
     void setMutualInducedTargetEpsilon(double inputMutualInducedTargetEpsilon);
 
     /**
-     * Set the coefficients for the mu_0, mu_1, mu_2, ..., mu_n terms in the perturbation
-     * theory algorithm for induced dipoles.
+     * Set the coefficients for the mu_0, mu_1, mu_2, ..., mu_n terms in the extrapolation
+     * algorithm for induced dipoles.
      *
-     * @param optCoefficients a vector whose mth entry specifies the coefficient for mu_m
+     * @param coefficients      a vector whose mth entry specifies the coefficient for mu_m.  The length of this
+     *                          vector determines how many iterations are performed.
      *
      */
-    void setOPTCoefficients(const std::vector<double> &OPTFullCoefficientsIn);
+    void setExtrapolationCoefficients(const std::vector<double> &coefficients);
 
     /**
-     * Get the coefficients for the mu_0, mu_1, mu_2, ..., mu_n terms in the perturbation
-     * theory algorithm for induced dipoles.
+     * Get the coefficients for the mu_0, mu_1, mu_2, ..., mu_n terms in the extrapolation
+     * algorithm for induced dipoles.  In this release, the default values for the coefficients are
+     * [0, -0.3, 0, 1.3], but be aware that those may change in a future release.
      */
-    const std::vector<double>& getOPTCoefficients() const;
+    const std::vector<double>& getExtrapolationCoefficients() const;
 
     /**
      * Get the error tolerance for Ewald summation.  This corresponds to the fractional error in the forces
@@ -405,7 +412,7 @@ private:
     int pmeBSplineOrder;
     std::vector<int> pmeGridDimension;
     int mutualInducedMaxIterations;
-    std::vector<double>  OPTFullCoefficients;
+    std::vector<double> extrapolationCoefficients;
 
     double mutualInducedTargetEpsilon;
     double scalingDistanceCutoff;

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors:                                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -43,6 +43,10 @@ AmoebaMultipoleForce::AmoebaMultipoleForce() : nonbondedMethod(NoCutoff), polari
                                                mutualInducedTargetEpsilon(1.0e-02), scalingDistanceCutoff(100.0), electricConstant(138.9354558456), aewald(0.0) {
     pmeGridDimension.resize(3);
     pmeGridDimension[0] = pmeGridDimension[1] = pmeGridDimension[2];
+    extrapolationCoefficients.push_back(0.0);
+    extrapolationCoefficients.push_back(-0.3);
+    extrapolationCoefficients.push_back(0.0);
+    extrapolationCoefficients.push_back(1.3);
 }
 
 AmoebaMultipoleForce::NonbondedMethod AmoebaMultipoleForce::getNonbondedMethod() const {
@@ -61,19 +65,14 @@ void AmoebaMultipoleForce::setPolarizationType(AmoebaMultipoleForce::Polarizatio
     polarizationType = type;
 }
 
-void AmoebaMultipoleForce::setOPTCoefficients(const std::vector<double> &OPTFullCoefficientsIn)
-{
-    size_t maxPTOrder = OPTFullCoefficientsIn.size();
-    OPTFullCoefficients.resize(maxPTOrder);
-    std::copy(OPTFullCoefficientsIn.begin(), OPTFullCoefficientsIn.end(), OPTFullCoefficients.begin());
+void AmoebaMultipoleForce::setExtrapolationCoefficients(const std::vector<double> &coefficients) {
+    extrapolationCoefficients = coefficients;
 }
 
-const std::vector<double> & AmoebaMultipoleForce::getOPTCoefficients() const
-{
-    return OPTFullCoefficients;
+const std::vector<double> & AmoebaMultipoleForce::getExtrapolationCoefficients() const {
+    return extrapolationCoefficients;
 }
 
-
 double AmoebaMultipoleForce::getCutoffDistance() const {
     return cutoffDistance;
 }

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

@@ -565,8 +565,8 @@ void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, c
     if (polarizationType == AmoebaMultipoleForce::Mutual) {
         mutualInducedMaxIterations = force.getMutualInducedMaxIterations();
         mutualInducedTargetEpsilon = force.getMutualInducedTargetEpsilon();
-    } else if (polarizationType == AmoebaMultipoleForce::OPT) {
-        OPTFullCoefficients = force.getOPTCoefficients();
+    } else if (polarizationType == AmoebaMultipoleForce::Extrapolated) {
+        extrapolationCoefficients = force.getExtrapolationCoefficients();
     }
 
     // PME
@@ -669,9 +669,9 @@ AmoebaReferenceMultipoleForce* ReferenceCalcAmoebaMultipoleForceKernel::setupAmo
         amoebaReferenceMultipoleForce->setMaximumMutualInducedDipoleIterations(mutualInducedMaxIterations);
     } else if (polarizationType == AmoebaMultipoleForce::Direct) {
         amoebaReferenceMultipoleForce->setPolarizationType(AmoebaReferenceMultipoleForce::Direct);
-    } else if (polarizationType == AmoebaMultipoleForce::OPT) {
-        amoebaReferenceMultipoleForce->setPolarizationType(AmoebaReferenceMultipoleForce::OPT);
-        amoebaReferenceMultipoleForce->setOPTCoefficients(OPTFullCoefficients);
+    } else if (polarizationType == AmoebaMultipoleForce::Extrapolated) {
+        amoebaReferenceMultipoleForce->setPolarizationType(AmoebaReferenceMultipoleForce::Extrapolated);
+        amoebaReferenceMultipoleForce->setExtrapolationCoefficients(extrapolationCoefficients);
     } else {
         throw OpenMMException("Polarization type not recognzied.");
     }

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

@@ -432,7 +432,7 @@ private:
 
     int mutualInducedMaxIterations;
     RealOpenMM mutualInducedTargetEpsilon;
-    std::vector<double> OPTFullCoefficients;
+    std::vector<double> extrapolationCoefficients;
 
     bool usePme;
     RealOpenMM alphaEwald;

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

@@ -161,16 +161,15 @@ RealOpenMM AmoebaReferenceMultipoleForce::getMutualInducedDipoleTargetEpsilon()
     return _mutualInducedDipoleTargetEpsilon;
 }
 
-void AmoebaReferenceMultipoleForce::setOPTCoefficients(const std::vector<RealOpenMM> &OPTFullCoefficients)
+void AmoebaReferenceMultipoleForce::setExtrapolationCoefficients(const std::vector<RealOpenMM> &coefficients)
 {
-    _maxPTOrder = OPTFullCoefficients.size(); // This accounts for the zero-based counting; actual highest order is 1 less
-    _OPTFullCoefficients.resize(_maxPTOrder);
-    _OPTPartCoefficients.resize(_maxPTOrder);
-    std::copy(OPTFullCoefficients.begin(), OPTFullCoefficients.end(), _OPTFullCoefficients.begin());
-    for(int i = 0; i < _maxPTOrder; ++i){
-        _OPTPartCoefficients[i] = 0.0;
-        for(int j = i; j < _maxPTOrder; ++j)
-            _OPTPartCoefficients[i] += _OPTFullCoefficients[j];
+    _maxPTOrder = coefficients.size(); // This accounts for the zero-based counting; actual highest order is 1 less
+    _extrapolationCoefficients = coefficients;
+    _extPartCoefficients.resize(_maxPTOrder);
+    for (int i = 0; i < _maxPTOrder; ++i) {
+        _extPartCoefficients[i] = 0.0;
+        for (int j = i; j < _maxPTOrder; ++j)
+            _extPartCoefficients[i] += _extrapolationCoefficients[j];
     }
 }
 
@@ -568,7 +567,7 @@ void AmoebaReferenceMultipoleForce::formQIRotationMatrix(const RealVec& iPositio
 {
     RealVec vectorZ = (deltaR)/r;
     RealVec vectorX(vectorZ);
-    if ((iPosition[1] != jPosition[1]) || (iPosition[2] != jPosition[2])){
+    if ((iPosition[1] != jPosition[1]) || (iPosition[2] != jPosition[2])) {
         vectorX[0] += 1.0;
     }else{
         vectorX[1] += 1.0;
@@ -813,8 +812,8 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
     RealVec deltaR       = particleJ.position - particleI.position;
     RealOpenMM r         =  SQRT(deltaR.dot(deltaR));
     vector<RealOpenMM> rrI(2);
-    // If we're using the OPT algorithm, we need to compute the field gradient, so ask for one more rrI value.
-    if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT)
+    // If we're using the extrapolation algorithm, we need to compute the field gradient, so ask for one more rrI value.
+    if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated)
         rrI.push_back(0.0);
   
     getAndScaleInverseRs(particleI.dampingFactor, particleJ.dampingFactor,
@@ -826,7 +825,7 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
     for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++) {
         calculateInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, rr3, rr5, deltaR,
                                        *updateInducedDipoleFields[ii].inducedDipoles, updateInducedDipoleFields[ii].inducedDipoleField);
-        if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
+        if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
             // Compute and store the field gradient for later use.
             RealOpenMM dx = deltaR[0];
             RealOpenMM dy = deltaR[1];
@@ -959,13 +958,7 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesBySOR(const vector<Mult
 }
 
 
-void AmoebaReferenceMultipoleForce::convergeInduceDipolesByOPT(const vector<MultipoleParticleData>& particleData, vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleField) {
-    if (_OPTFullCoefficients.empty()){
-        std::stringstream message;
-        message << "An OPT calcultion was requested, but setOPTCoefficients() was not called.";
-        throw OpenMMException(message.str());
-    }
-
+void AmoebaReferenceMultipoleForce::convergeInduceDipolesByExtrapolation(const vector<MultipoleParticleData>& particleData, vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleField) {
     _ptDipoleD.clear();
     _ptDipoleP.clear();
 
@@ -975,7 +968,7 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByOPT(const vector<Mult
     // Start by storing the direct dipoles as PT0
     vector<RealVec> thisDipoleD;
     vector<RealVec> thisDipoleP;
-    for(int atom = 0; atom < _numParticles; ++atom){
+    for (int atom = 0; atom < _numParticles; ++atom) {
         thisDipoleD.push_back((*fieldD.inducedDipoles)[atom]);
         thisDipoleP.push_back((*fieldP.inducedDipoles)[atom]);
     }
@@ -988,13 +981,13 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByOPT(const vector<Mult
     fieldP.inducedDipoleFieldGradient.resize(_numParticles);
 
     // Recursively apply alpha.Tau to the µ_(n) components to generate µ_(n+1), and store the result
-    for(int order = 1; order < _maxPTOrder; ++order){
+    for (int order = 1; order < _maxPTOrder; ++order) {
         std::fill(fieldD.inducedDipoleFieldGradient.begin(), fieldD.inducedDipoleFieldGradient.end(), zeros);
         std::fill(fieldP.inducedDipoleFieldGradient.begin(), fieldP.inducedDipoleFieldGradient.end(), zeros);
         calculateInducedDipoleFields(particleData, updateInducedDipoleField);
         vector<RealVec> thisDipoleD;
         vector<RealVec> thisDipoleP;
-        for(int atom = 0; atom < _numParticles; ++atom){
+        for (int atom = 0; atom < _numParticles; ++atom) {
             (*fieldD.inducedDipoles)[atom] = fieldD.inducedDipoleField[atom] * particleData[atom].polarity;
             (*fieldP.inducedDipoles)[atom] = fieldP.inducedDipoleField[atom] * particleData[atom].polarity;
             thisDipoleD.push_back((*fieldD.inducedDipoles)[atom]);
@@ -1004,8 +997,8 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByOPT(const vector<Mult
         _ptDipoleP.push_back(thisDipoleP);
         vector<RealOpenMM> fieldGradD(6*_numParticles, 0.0);
         vector<RealOpenMM> fieldGradP(6*_numParticles, 0.0);
-        for(int atom = 0; atom < _numParticles; ++atom){
-            for(int component = 0; component < 6; ++component){
+        for (int atom = 0; atom < _numParticles; ++atom) {
+            for (int component = 0; component < 6; ++component) {
                 fieldGradD[6*atom + component] = fieldD.inducedDipoleFieldGradient[atom][component];
                 fieldGradP[6*atom + component] = fieldP.inducedDipoleFieldGradient[atom][component];
             }
@@ -1019,16 +1012,16 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByOPT(const vector<Mult
     std::fill(_inducedDipole.begin(), _inducedDipole.end(), zeroVec);
     std::fill(_inducedDipolePolar.begin(), _inducedDipolePolar.end(), zeroVec);
 
-    for(int order = 0; order < _maxPTOrder; ++order){
-        for(int atom = 0; atom < _numParticles; ++atom){
-            _inducedDipole[atom]      +=  _ptDipoleD[order][atom] * _OPTPartCoefficients[order];
-            _inducedDipolePolar[atom] +=  _ptDipoleP[order][atom] * _OPTPartCoefficients[order];
+    for (int order = 0; order < _maxPTOrder; ++order) {
+        for (int atom = 0; atom < _numParticles; ++atom) {
+            _inducedDipole[atom]      +=  _ptDipoleD[order][atom] * _extPartCoefficients[order];
+            _inducedDipolePolar[atom] +=  _ptDipoleP[order][atom] * _extPartCoefficients[order];
         }
     }
 
     // Copy the combined dipoles over to compute the field
-    for(int order = 0; order < _maxPTOrder; ++order){
-        for(int atom = 0; atom < _numParticles; ++atom){
+    for (int order = 0; order < _maxPTOrder; ++order) {
+        for (int atom = 0; atom < _numParticles; ++atom) {
             (*fieldD.inducedDipoles)[atom] = _inducedDipole[atom];
             (*fieldP.inducedDipoles)[atom] = _inducedDipolePolar[atom];
         }
@@ -1177,10 +1170,10 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipoles(const vector<Multipo
 
     // UpdateInducedDipoleFieldStruct contains induced dipole, fixed multipole fields and fields
     // due to other induced dipoles at each site
-    if (getPolarizationType() == AmoebaReferenceMultipoleForce::Mutual){
+    if (getPolarizationType() == AmoebaReferenceMultipoleForce::Mutual) {
         convergeInduceDipolesByDIIS(particleData, updateInducedDipoleField);
-    }if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
-        convergeInduceDipolesByOPT(particleData, updateInducedDipoleField);
+    }if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
+        convergeInduceDipolesByExtrapolation(particleData, updateInducedDipoleField);
     }
 
 }
@@ -1298,7 +1291,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPairIxn(const Mu
     // The rInvVec array is defined such that the ith element is R^-i, with the
     // dieleectric constant folded in, to avoid conversions later.
     rInvVec[1] = prefac * rInv;
-    for(int i = 2; i < 7; ++i)
+    for (int i = 2; i < 7; ++i)
         rInvVec[i] = rInvVec[i-1] * rInv;
 
     RealOpenMM mScale = scalingFactors[M_SCALE];
@@ -1504,7 +1497,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPairIxn(const Mu
     RealOpenMM fIZ = qiQI[0]*VijR[0];
     RealOpenMM fJZ = qiQJ[0]*VjiR[0];
     RealOpenMM EIX = 0.0, EIY = 0.0, EIZ = 0.0, EJX = 0.0, EJY = 0.0, EJZ = 0.0;
-    for(int i = 1; i < 9; ++i){
+    for (int i = 1; i < 9; ++i) {
         energy += 0.5*(qiQI[i]*Vij[i] + qiQJ[i]*Vji[i]);
         fIZ += qiQI[i]*VijR[i];
         fJZ += qiQJ[i]*VjiR[i];
@@ -1532,7 +1525,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPairIxn(const Mu
     RealOpenMM iEJY = qiUinpJ[0]*Vjip[1] + qiUindJ[0]*Vjid[1] - qiUinpJ[1]*Vjip[0] - qiUindJ[1]*Vjid[0];
 
     // Add in the induced-induced terms, if needed.
-    if(getPolarizationType() == AmoebaReferenceMultipoleForce::Mutual){
+    if(getPolarizationType() == AmoebaReferenceMultipoleForce::Mutual) {
         // Uind-Uind terms (m=0)
         RealOpenMM eCoef = -4.0*rInvVec[3]*uScale*thole_d0;
         RealOpenMM dCoef = 6.0*rInvVec[4]*uScale*dthole_d0;
@@ -1846,13 +1839,13 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostatic(const vector<Mu
             }
         }
     }
-    if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
+    if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
         RealOpenMM prefac = (_electric/_dielectric);
         for (int i = 0; i < _numParticles; i++) {
             // Compute the µ(m) T µ(n) force contributions here
-            for(int l = 0; l < _maxPTOrder-1; ++l) {
-                for(int m = 0; m < _maxPTOrder-1-l; ++m) {
-                    RealOpenMM p = _OPTPartCoefficients[l+m+1];
+            for (int l = 0; l < _maxPTOrder-1; ++l) {
+                for (int m = 0; m < _maxPTOrder-1-l; ++m) {
+                    RealOpenMM p = _extPartCoefficients[l+m+1];
                     if(std::fabs(p) < 1e-6) continue;
                     forces[i][0] += 0.5*p*prefac*(_ptDipoleD[l][i][0]*_ptDipoleFieldGradientP[m][6*i+0]
                                                 + _ptDipoleD[l][i][1]*_ptDipoleFieldGradientP[m][6*i+3]
@@ -6091,7 +6084,7 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
 
     calculateReciprocalSpaceInducedDipoleField(updateInducedDipoleFields);
 
-    if(getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
+    if(getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
         // While we have the reciprocal space (fractional coordinate) field gradient available, add it to the real space
         // terms computed above, after transforming to Cartesian coordinates.  This allows real and reciprocal space
         // dipole response force contributions to be computed together.
@@ -6109,8 +6102,8 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
             };
 
             RealOpenMM Exx = 0.0, Eyy = 0.0, Ezz = 0.0, Exy = 0.0, Exz = 0.0, Eyz = 0.0;
-            for(int k = 0; k < 3; ++k){
-                for(int l = 0; l < 3; ++l){
+            for (int k = 0; k < 3; ++k) {
+                for (int l = 0; l < 3; ++l) {
                     Exx += fracToCart[0][k] * EmatD[k][l] * fracToCart[0][l];
                     Eyy += fracToCart[1][k] * EmatD[k][l] * fracToCart[1][l];
                     Ezz += fracToCart[2][k] * EmatD[k][l] * fracToCart[2][l];
@@ -6133,8 +6126,8 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
             };
 
             Exx = 0.0; Eyy = 0.0; Ezz = 0.0; Exy = 0.0; Exz = 0.0; Eyz = 0.0;
-            for(int k = 0; k < 3; ++k){
-                for(int l = 0; l < 3; ++l){
+            for (int k = 0; k < 3; ++k) {
+                for (int l = 0; l < 3; ++l) {
                     Exx += fracToCart[0][k] * EmatP[k][l] * fracToCart[0][l];
                     Eyy += fracToCart[1][k] * EmatP[k][l] * fracToCart[1][l];
                     Ezz += fracToCart[2][k] * EmatP[k][l] * fracToCart[2][l];
@@ -6260,7 +6253,7 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
         calculateDirectInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, preFactor1, preFactor2, deltaR,
                                             *updateInducedDipoleFields[ii].inducedDipoles,
                                             updateInducedDipoleFields[ii].inducedDipoleField);
-        if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
+        if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
             // Compute and store the field gradient for later use.
             RealOpenMM dx = deltaR[0];
             RealOpenMM dy = deltaR[1];
@@ -6468,13 +6461,13 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeDirectElectrostaticPair
     // The rInvVec array is defined such that the ith element is R^-i, with the
     // dieleectric constant folded in, to avoid conversions later.
     rInvVec[1] = prefac * rInv;
-    for(int i = 2; i < 7; ++i)
+    for (int i = 2; i < 7; ++i)
         rInvVec[i] = rInvVec[i-1] * rInv;
 
     // The alpharVec array is defined such that the ith element is (alpha R)^i,
     // where kappa (alpha in OpenMM parlance) is the Ewald attenuation parameter.
     alphaRVec[1] = _alphaEwald * r;
-    for(int i = 2; i < 8; ++i)
+    for (int i = 2; i < 8; ++i)
         alphaRVec[i] = alphaRVec[i-1] * alphaRVec[1];
 
     RealOpenMM erfAlphaR = erf(alphaRVec[1]);
@@ -6487,7 +6480,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeDirectElectrostaticPair
     int doubleFactorial = 1, facCount = 1;
     RealOpenMM tmp = alphaRVec[1];
     bVec[1] = -erfAlphaR;
-    for(int i=2; i < 5; ++i){
+    for (int i=2; i < 5; ++i) {
         bVec[i] = bVec[i-1] + tmp * X / (RealOpenMM)(doubleFactorial);
         facCount = facCount + 2;
         doubleFactorial = doubleFactorial * facCount;
@@ -6692,7 +6685,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeDirectElectrostaticPair
     RealOpenMM fIZ = qiQI[0]*VijR[0];
     RealOpenMM fJZ = qiQJ[0]*VjiR[0];
     RealOpenMM EIX = 0.0, EIY = 0.0, EIZ = 0.0, EJX = 0.0, EJY = 0.0, EJZ = 0.0;
-    for(int i = 1; i < 9; ++i){
+    for (int i = 1; i < 9; ++i) {
         energy += 0.5*(qiQI[i]*Vij[i] + qiQJ[i]*Vji[i]);
         fIZ += qiQI[i]*VijR[i];
         fJZ += qiQJ[i]*VjiR[i];
@@ -6801,14 +6794,14 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector
     energy += calculatePmeSelfEnergy(particleData);
 
     // Now that both the direct and reciprocal space contributions have been added, we can compute the dipole
-    // response contributions to the forces, if we're using the OPT polarization algorithm.
-    if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
+    // response contributions to the forces, if we're using the extrapolated polarization algorithm.
+    if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
         RealOpenMM prefac = (_electric/_dielectric);
         for (int i = 0; i < _numParticles; i++) {
             // Compute the µ(m) T µ(n) force contributions here
-            for(int l = 0; l < _maxPTOrder-1; ++l) {
-                for(int m = 0; m < _maxPTOrder-1-l; ++m) {
-                    RealOpenMM p = _OPTPartCoefficients[l+m+1];
+            for (int l = 0; l < _maxPTOrder-1; ++l) {
+                for (int m = 0; m < _maxPTOrder-1-l; ++m) {
+                    RealOpenMM p = _extPartCoefficients[l+m+1];
                     if(std::fabs(p) < 1e-6) continue;
                     forces[i][0] += 0.5*p*prefac*(_ptDipoleD[l][i][0]*_ptDipoleFieldGradientP[m][6*i+0]
                                                 + _ptDipoleD[l][i][1]*_ptDipoleFieldGradientP[m][6*i+3]

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

@@ -351,9 +351,9 @@ public:
         Direct = 1,
 
         /**
-         * Optimized perturbation theory
+         * Extrapolated perturbation theory
          */
-        OPT = 2
+        Extrapolated = 2
     };
 
     /**
@@ -428,13 +428,13 @@ public:
     RealOpenMM getMutualInducedDipoleEpsilon() const;
 
     /**
-     * Set the coefficients for the µ_0, µ_1, µ_2, µ_n terms in the pertubation
+     * Set the coefficients for the µ_0, µ_1, µ_2, µ_n terms in the extrapolation
      * theory algorithm for induced dipoles
      *
      * @param optCoefficients a vector whose mth entry specifies the coefficient for µ_m
      *
      */
-    void setOPTCoefficients(const std::vector<RealOpenMM> &OPTFullCoefficients);
+    void setExtrapolationCoefficients(const std::vector<RealOpenMM> &coefficients);
 
     /**
      * Set the target epsilon for converging mutual induced dipoles.
@@ -677,8 +677,8 @@ protected:
     int _mutualInducedDipoleIterations;
     int _maximumMutualInducedDipoleIterations;
     int _maxPTOrder;
-    std::vector<RealOpenMM>  _OPTFullCoefficients;
-    std::vector<RealOpenMM>  _OPTPartCoefficients;
+    std::vector<RealOpenMM>  _extrapolationCoefficients;
+    std::vector<RealOpenMM>  _extPartCoefficients;
     RealOpenMM  _mutualInducedDipoleEpsilon;
     RealOpenMM  _mutualInducedDipoleTargetEpsilon;
     RealOpenMM  _polarSOR;
@@ -945,13 +945,13 @@ protected:
     virtual void calculateInducedDipoleFields(const std::vector<MultipoleParticleData>& particleData,
                                               std::vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields);
     /**
-     * Calculated induced dipoles using Optimized Perturbation Theory.
+     * Calculated induced dipoles using extrapolated perturbation theory.
      *
      * @param particleData              vector of particle positions and parameters (charge, labFrame dipoles, quadrupoles, ...)
      * @param updateInducedDipoleFields vector of UpdateInducedDipoleFieldStruct containing input induced dipoles and output fields
      */
-    void convergeInduceDipolesByOPT(const std::vector<MultipoleParticleData>& particleData,
-                                    std::vector<UpdateInducedDipoleFieldStruct>& calculateInducedDipoleField);
+    void convergeInduceDipolesByExtrapolation(const std::vector<MultipoleParticleData>& particleData,
+                                              std::vector<UpdateInducedDipoleFieldStruct>& calculateInducedDipoleField);
     /**
      * Converge induced dipoles.
      * 

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

@@ -319,13 +319,13 @@ static void testWaterDimerExPTPolarizationTriclinicPME() {
                                         Vec3(0.2, 2.0, 0.0),
                                         Vec3(0.1, 0.5, 2.0));
     amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::PME);
-    amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::OPT);
+    amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::Extrapolated);
     std::vector<double> coefs;
     coefs.push_back(0.0);  // The mu_0 coefficient
     coefs.push_back(-0.3); // The mu_1 coefficient
     coefs.push_back(0.0);  // The mu_2 coefficient
     coefs.push_back(1.3);  // The mu_3 coefficient
-    amoebaMultipoleForce->setOPTCoefficients(coefs);
+    amoebaMultipoleForce->setExtrapolationCoefficients(coefs);
     amoebaMultipoleForce->setCutoffDistance(9.0*OpenMM::NmPerAngstrom);
     amoebaMultipoleForce->setAEwald(4);
     amoebaMultipoleForce->setEwaldErrorTolerance(1.0e-06);
@@ -368,13 +368,13 @@ static void testWaterDimerExPTPolarizationTriclinicPMENoPolGroups() {
                                         Vec3(0.2, 2.0, 0.0),
                                         Vec3(0.1, 0.5, 2.0));
     amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::PME);
-    amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::OPT);
+    amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::Extrapolated);
     std::vector<double> coefs;
     coefs.push_back(0.0);  // The mu_0 coefficient
     coefs.push_back(-0.3); // The mu_1 coefficient
     coefs.push_back(0.0);  // The mu_2 coefficient
     coefs.push_back(1.3);  // The mu_3 coefficient
-    amoebaMultipoleForce->setOPTCoefficients(coefs);
+    amoebaMultipoleForce->setExtrapolationCoefficients(coefs);
     amoebaMultipoleForce->setCutoffDistance(9.0*OpenMM::NmPerAngstrom);
     amoebaMultipoleForce->setAEwald(4);
     amoebaMultipoleForce->setEwaldErrorTolerance(1.0e-06);
@@ -415,13 +415,13 @@ static void testWaterDimerExPTPolarizationNoCutoff() {
     vector<Vec3> coords = setupWaterDimer(system, amoebaMultipoleForce, true);
 
     amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::NoCutoff);
-    amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::OPT);
+    amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::Extrapolated);
     std::vector<double> coefs;
     coefs.push_back(0.0);  // The mu_0 coefficient
     coefs.push_back(-0.3); // The mu_1 coefficient
     coefs.push_back(0.0);  // The mu_2 coefficient
     coefs.push_back(1.3);  // The mu_3 coefficient
-    amoebaMultipoleForce->setOPTCoefficients(coefs);
+    amoebaMultipoleForce->setExtrapolationCoefficients(coefs);
 
     LangevinIntegrator integrator(0.0, 0.1, 0.01);
     Context context(system, integrator, Platform::getPlatformByName("Reference"));
@@ -456,13 +456,13 @@ static void testWaterDimerExPTPolarizationNoCutoffNoPolGroups() {
     vector<Vec3> coords = setupWaterDimer(system, amoebaMultipoleForce, false);
 
     amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::NoCutoff);
-    amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::OPT);
+    amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::Extrapolated);
     std::vector<double> coefs;
     coefs.push_back(0.0);  // The mu_0 coefficient
     coefs.push_back(-0.3); // The mu_1 coefficient
     coefs.push_back(0.0);  // The mu_2 coefficient
     coefs.push_back(1.3);  // The mu_3 coefficient
-    amoebaMultipoleForce->setOPTCoefficients(coefs);
+    amoebaMultipoleForce->setExtrapolationCoefficients(coefs);
 
     LangevinIntegrator integrator(0.0, 0.1, 0.01);
     Context context(system, integrator, Platform::getPlatformByName("Reference"));
@@ -492,7 +492,6 @@ static void testWaterDimerExPTPolarizationNoCutoffNoPolGroups() {
 int main(int numberOfArguments, char* argv[]) {
 
     try {
-        std::cout << "TestReferenceAmoebaPTPolarization running test..." << std::endl;
         registerAmoebaReferenceKernelFactories();
 
         /*

plugins/amoeba/serialization/src/AmoebaMultipoleForceProxy.cpp

@@ -68,7 +68,7 @@ void loadCovalentMap(const SerializationNode& map, std::vector< int >& covalentM
 }
 
 void AmoebaMultipoleForceProxy::serialize(const void* object, SerializationNode& node) const {
-    node.setIntProperty("version", 2);
+    node.setIntProperty("version", 3);
     const AmoebaMultipoleForce& force = *reinterpret_cast<const AmoebaMultipoleForce*>(object);
 
     node.setIntProperty("nonbondedMethod",                  force.getNonbondedMethod());
@@ -87,6 +87,14 @@ void AmoebaMultipoleForceProxy::serialize(const void* object, SerializationNode&
     force.getPmeGridDimensions(gridDimensions);
     SerializationNode& gridDimensionsNode  = node.createChildNode("MultipoleParticleGridDimension");
     gridDimensionsNode.setIntProperty("d0", gridDimensions[0]).setIntProperty("d1", gridDimensions[1]).setIntProperty("d2", gridDimensions[2]); 
+    
+    SerializationNode& coefficients = node.createChildNode("ExtrapolationCoefficients");
+    vector<double> coeff = force.getExtrapolationCoefficients();
+    for (int i = 0; i < coeff.size(); i++) {
+        stringstream key;
+        key << "c" << i;
+        coefficients.setDoubleProperty(key.str(), coeff[i]);
+    }
 
     std::vector<std::string> covalentTypes;
     getCovalentTypes(covalentTypes);
@@ -124,16 +132,16 @@ void AmoebaMultipoleForceProxy::serialize(const void* object, SerializationNode&
 }
 
 void* AmoebaMultipoleForceProxy::deserialize(const SerializationNode& node) const {
-    if (node.getIntProperty("version") > 2)
+    int version = node.getIntProperty("version");
+    if (version < 0 || version > 3)
         throw OpenMMException("Unsupported version number");
     AmoebaMultipoleForce* force = new AmoebaMultipoleForce();
 
     try {
 
         force->setNonbondedMethod(static_cast<AmoebaMultipoleForce::NonbondedMethod>(node.getIntProperty("nonbondedMethod")));
-        if (node.getIntProperty("version") == 2) {
+        if (version >= 2)
             force->setPolarizationType(static_cast<AmoebaMultipoleForce::PolarizationType>(node.getIntProperty("polarizationType")));
-        }
         //force->setPmeBSplineOrder(node.getIntProperty("pmeBSplineOrder"));
         //force->setMutualInducedIterationMethod(static_cast<AmoebaMultipoleForce::MutualInducedIterationMethod>(node.getIntProperty("mutualInducedIterationMethod")));
         force->setMutualInducedMaxIterations(node.getIntProperty("mutualInducedMaxIterations"));
@@ -151,6 +159,18 @@ void* AmoebaMultipoleForceProxy::deserialize(const SerializationNode& node) cons
         gridDimensions.push_back(gridDimensionsNode.getIntProperty("d2"));
         force->setPmeGridDimensions(gridDimensions);
     
+        if (version >= 3) {
+            const SerializationNode& coefficients = node.getChildNode("ExtrapolationCoefficients");
+            vector<double> coeff;
+            for (int i = 0; ; i++) {
+                stringstream key;
+                key << "c" << i;
+                if (coefficients.getProperties().find(key.str()) == coefficients.getProperties().end())
+                    break;
+                coeff.push_back(coefficients.getDoubleProperty(key.str()));
+            }
+            force->setExtrapolationCoefficients(coeff);
+        }
         std::vector<std::string> covalentTypes;
         getCovalentTypes(covalentTypes);
 

plugins/amoeba/serialization/tests/TestSerializeAmoebaMultipoleForce.cpp

@@ -74,6 +74,12 @@ void testSerialization() {
     force1.setMutualInducedTargetEpsilon(1.0e-05); 
     //force1.setElectricConstant(138.93); 
     force1.setEwaldErrorTolerance(1.0e-05); 
+    
+    vector<double> coeff;
+    coeff.push_back(0.0);
+    coeff.push_back(-0.1);
+    coeff.push_back(1.1);
+    force1.setExtrapolationCoefficients(coeff);
 
     std::vector<std::string> covalentTypes;
     getCovalentTypes(covalentTypes);
@@ -125,6 +131,8 @@ void testSerialization() {
         ASSERT_EQUAL(gridDimension1[jj], gridDimension2[jj]);
     }
     
+    ASSERT_EQUAL_CONTAINERS(force1.getExtrapolationCoefficients(), force2.getExtrapolationCoefficients());
+    
     ASSERT_EQUAL(force1.getNumMultipoles(),  force2.getNumMultipoles());
     for (unsigned int ii = 0; ii < static_cast<unsigned int>(force1.getNumMultipoles()); ii++) {