Added hooks for the Optimized Perturbation Theory Reference impl.

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

@@ -79,7 +79,13 @@ public:
         /**
          * Direct polarization
          */
-        Direct = 1
+        Direct = 1,
+
+        /**
+         * Optimized perturbation theory
+         */
+        OPT = 2
+
     };
 
     enum MultipoleAxisTypes { ZThenX = 0, Bisector = 1, ZBisect = 2, ThreeFold = 3, ZOnly = 4, NoAxisType = 5, LastAxisTypeIndex = 6 };
@@ -296,6 +302,24 @@ public:
      */
     void setMutualInducedTargetEpsilon(double inputMutualInducedTargetEpsilon);
 
+    /**
+     * Set the coefficients for the µ_0, µ_1, µ_2, µ_n terms in the pertubation
+     * theory algorithm for induced dipoles
+     *
+     * @param optCoefficients a vector whose mth entry specifies the coefficient for µ_m
+     *
+     */
+    void setOPTCoefficients(const std::vector<double> &OPTFullCoefficientsIn);
+
+    /**
+     * Set the coefficients for the µ_0, µ_1, µ_2, µ_n terms in the pertubation
+     * theory algorithm for induced dipoles
+     *
+     * @param optCoefficients a vector whose mth entry specifies the coefficient for µ_m
+     *
+     */
+    const std::vector<double>& getOPTCoefficients() const;
+
     /**
      * Get the error tolerance for Ewald summation.  This corresponds to the fractional error in the forces
      * which is acceptable.  This value is used to select the grid dimensions and separation (alpha)
@@ -382,6 +406,8 @@ private:
     int pmeBSplineOrder;
     std::vector<int> pmeGridDimension;
     int mutualInducedMaxIterations;
+    std::vector<double>  OPTFullCoefficients;
+
     double mutualInducedTargetEpsilon;
     double scalingDistanceCutoff;
     double electricConstant;

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -61,6 +61,19 @@ 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());
+}
+
+const std::vector<double> & AmoebaMultipoleForce::getOPTCoefficients() const
+{
+    return OPTFullCoefficients;
+}
+
+
 double AmoebaMultipoleForce::getCutoffDistance() const {
     return cutoffDistance;
 }

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

@@ -565,6 +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();
     }
 
     // PME
@@ -667,6 +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 {
         throw OpenMMException("Polarization type not recognzied.");
     }

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

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

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

@@ -162,6 +162,19 @@ RealOpenMM AmoebaReferenceMultipoleForce::getMutualInducedDipoleTargetEpsilon()
     return _mutualInducedDipoleTargetEpsilon;
 }
 
+void AmoebaReferenceMultipoleForce::setOPTCoefficients(const std::vector<RealOpenMM> &OPTFullCoefficients)
+{
+    _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];
+    }
+}
+
 void AmoebaReferenceMultipoleForce::setMutualInducedDipoleTargetEpsilon(RealOpenMM mutualInducedDipoleTargetEpsilon)
 {
     _mutualInducedDipoleTargetEpsilon = mutualInducedDipoleTargetEpsilon;

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

@@ -348,7 +348,12 @@ public:
         /**
          * Direct polarization
          */
-        Direct = 1 
+        Direct = 1,
+
+        /**
+         * Optimized perturbation theory
+         */
+        OPT = 2
     };
 
     /**
@@ -422,6 +427,15 @@ public:
      */
     RealOpenMM getMutualInducedDipoleEpsilon() const;
 
+    /**
+     * Set the coefficients for the µ_0, µ_1, µ_2, µ_n terms in the pertubation
+     * theory algorithm for induced dipoles
+     *
+     * @param optCoefficients a vector whose mth entry specifies the coefficient for µ_m
+     *
+     */
+    void setOPTCoefficients(const std::vector<RealOpenMM> &OPTFullCoefficients);
+
     /**
      * Set the target epsilon for converging mutual induced dipoles.
      *
@@ -655,6 +669,9 @@ protected:
     int _mutualInducedDipoleConverged;
     int _mutualInducedDipoleIterations;
     int _maximumMutualInducedDipoleIterations;
+    int _maxPTOrder;
+    std::vector<RealOpenMM>  _OPTFullCoefficients;
+    std::vector<RealOpenMM>  _OPTPartCoefficients;
     RealOpenMM  _mutualInducedDipoleEpsilon;
     RealOpenMM  _mutualInducedDipoleTargetEpsilon;
     RealOpenMM  _polarSOR;