Added code Amoeba Reference platform to handle polarization type = direct

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

@@ -445,6 +445,7 @@ double ReferenceCalcAmoebaTorsionTorsionForceKernel::execute(ContextImpl& contex
 
 ReferenceCalcAmoebaMultipoleForceKernel::ReferenceCalcAmoebaMultipoleForceKernel(std::string name, const Platform& platform, System& system) : 
          CalcAmoebaMultipoleForceKernel(name, platform), system(system) {
+
 }
 
 ReferenceCalcAmoebaMultipoleForceKernel::~ReferenceCalcAmoebaMultipoleForceKernel() {
@@ -521,6 +522,12 @@ void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, c
     if( nonbondedMethod != 0 && nonbondedMethod != 1 ){
          throw OpenMMException("AmoebaMultipoleForce nonbonded method not recognized.\n");
     }
+
+    polarizationType = static_cast<int>(force.getPolarizationType());
+    if( polarizationType != 0 && polarizationType != 1 ){ 
+         throw OpenMMException("AmoebaMultipoleForce polarization type not recognized.\n");
+    }    
+
 }
 
 double ReferenceCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -535,7 +542,7 @@ double ReferenceCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bo
                                                                                     charges, dipoles, quadrupoles, tholes,
                                                                                     dampingFactors, polarity, axisTypes, 
                                                                                     multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
-                                                                                    multipoleAtomCovalentInfo, forceData);
+                                                                                    multipoleAtomCovalentInfo, polarizationType, forceData);
 
     return static_cast<double>(energy);
 }

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

@@ -386,6 +386,7 @@ public:
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
     int numMultipoles;
+    int polarizationType;
     std::vector<RealOpenMM> charges;
     std::vector<RealOpenMM> dipoles;
     std::vector<RealOpenMM> quadrupoles;

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

@@ -945,7 +945,7 @@ void AmoebaReferenceMultipoleForce::updateInducedDipole(  MultipoleParticleData&
     }
 }
 
-void AmoebaReferenceMultipoleForce::calculateNoCutoffInducedDipoles( std::vector<MultipoleParticleData>& particleData ){
+void AmoebaReferenceMultipoleForce::calculateNoCutoffInducedDipoles( int polarizationType, std::vector<MultipoleParticleData>& particleData ){
 
     // ---------------------------------------------------------------------------------------
 
@@ -984,6 +984,10 @@ void AmoebaReferenceMultipoleForce::calculateNoCutoffInducedDipoles( std::vector
         particleData[ii].inducedDipolePolar[2]     = particleData[ii].fieldPolar[2]; 
     }
 
+    if( polarizationType == 1 ){
+        return;
+    }
+
     std::vector<RealOpenMM> field( numParticles*3 );
     std::vector<RealOpenMM> fieldPolar( numParticles*3 );
 
@@ -1037,7 +1041,8 @@ void AmoebaReferenceMultipoleForce::calculateNoCutoffInducedDipoles( std::vector
     return;
 }
 
-RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostaticPairIxn( const MultipoleParticleData& particleI,
+RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostaticPairIxn( int polarizationType,
+                                                                                 const MultipoleParticleData& particleI,
                                                                                  const MultipoleParticleData& particleK,
                                                                                  RealOpenMM* scalingFactors, std::vector<RealVec>& forces,
                                                                                  std::vector<Vec3>& torque ) const {
@@ -1075,6 +1080,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostaticPairIxn(
     RealOpenMM gl[9],gli[7],glip[7];
     RealOpenMM sc[10],sci[8],scip[8];
     RealOpenMM gf[7],gfi[6],gti[6];
+    RealOpenMM gfd, fdir[3];
  
     RealOpenMM delta[3];
     getDelta( particleI, particleK, delta );
@@ -1455,8 +1461,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostaticPairIxn(
 
     // correction to convert mutual to direct polarization force
 
-/*
-    if( poltyp .eq. 'DIRECT'){
+    if( polarizationType == 1 ){
        gfd = 0.5 * (rr5*scip[1]*scale3i
              - rr7*(scip[2]*sci[3]+sci[2]*scip[3])*scale5i);
        temp5 = 0.5 * rr5 * scale5i;
@@ -1473,7 +1478,6 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostaticPairIxn(
        ftm2i[1] = ftm2i[1] - fdir[1] + findmp[1];
        ftm2i[2] = ftm2i[2] - fdir[2] + findmp[2];
     }
-*/
 
     // intermediate terms for induced torque on multipoles
 
@@ -1925,6 +1929,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostatic( std::v
                                                                           const std::vector<int>& multipoleAtomZs,
                                                                           const std::vector<int>& multipoleAtomXs,
                                                                           const std::vector<int>& multipoleAtomYs,
+                                                                          int polarizationType,
                                                                           std::vector<RealVec>& forces ) const {
 
     // ---------------------------------------------------------------------------------------
@@ -1977,7 +1982,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostatic( std::v
                 getScaleFactors( ii, jj, scaleFactors);
             }
 
-            energy += calculateNoCutoffElectrostaticPairIxn( particleData[ii], particleData[jj], scaleFactors, forces, torques );
+            energy += calculateNoCutoffElectrostaticPairIxn( polarizationType, particleData[ii], particleData[jj], scaleFactors, forces, torques );
 
             if( jj <= _maxScaleIndex[ii] ){
                 for( unsigned int kk = 0; kk < LAST_SCALE_TYPE_INDEX; kk++ ){
@@ -2049,6 +2054,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffForceAndEnergy( const
                                                                            const std::vector<int>& multipoleAtomZs,
                                                                            const std::vector<int>& multipoleAtomXs,
                                                                            const std::vector<int>& multipoleAtomYs,
+                                                                           int polarizationType,
                                                                            std::vector<RealVec>& forces ){
 
 
@@ -2133,11 +2139,11 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffForceAndEnergy( const
 
     // get induced dipoles
 
-    calculateNoCutoffInducedDipoles( particleData );
+    calculateNoCutoffInducedDipoles( polarizationType, particleData );
 
     // check if induced dipoles converged
 
-    if( !getMutualInducedDipoleConverged() ){
+    if( !getMutualInducedDipoleConverged() && polarizationType == 0 ){
         std::stringstream message;
         message << "Induced dipoles did not converge: ";
         message << " iterations="      << getMutualInducedDipoleIterations();
@@ -2145,7 +2151,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffForceAndEnergy( const
         throw OpenMMException(message.str());
     }
 
-    RealOpenMM energy = calculateNoCutoffElectrostatic( particleData, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, forces );
+    RealOpenMM energy = calculateNoCutoffElectrostatic( particleData, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, polarizationType, forces );
 
     return energy;
 
@@ -2164,6 +2170,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateForceAndEnergy(
                                                          const std::vector<int>& multipoleAtomXs,
                                                          const std::vector<int>& multipoleAtomYs,
                                                          const std::vector< std::vector< std::vector<int> > >& multipoleParticleCovalentInfo,
+                                                         int polarizationType,
                                                          std::vector<RealVec>& forces ){
 
     
@@ -2171,7 +2178,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateForceAndEnergy(
     if( getNonbondedMethod() == NoCutoff || 1 ){
 
         return calculateNoCutoffForceAndEnergy( particlePositions, charges, dipoles, quadrupoles, tholes, dampingFactors,
-                                                polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, forces );
+                                                polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, polarizationType, forces );
 
 
     }    

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

@@ -215,6 +215,7 @@ public:
                                         const std::vector<int>& multipoleAtomXs,
                                         const std::vector<int>& multipoleAtomYs,
                                         const std::vector< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo,
+                                        int polarizationType,
                                         std::vector<OpenMM::RealVec>& forces );
 
 
@@ -456,13 +457,13 @@ private:
                               const std::vector<RealOpenMM>& fieldPolar,
                               RealOpenMM& epsilonD, RealOpenMM& epsilonP ) const;
 
-    void calculateNoCutoffInducedDipoles( std::vector<MultipoleParticleData>& particleData );
+    void calculateNoCutoffInducedDipoles( int polarizationType, std::vector<MultipoleParticleData>& particleData );
 
     void calculateInducedDipoleField( std::vector<MultipoleParticleData>& particleData,
                                       std::vector<RealOpenMM>& field,
                                       std::vector<RealOpenMM>& fieldPolar ) const;
 
-    RealOpenMM calculateNoCutoffElectrostaticPairIxn( const MultipoleParticleData& particleI,
+    RealOpenMM calculateNoCutoffElectrostaticPairIxn( int polarizationType, const MultipoleParticleData& particleI,
                                                       const MultipoleParticleData& particleK,
                                                       RealOpenMM* scalingFactors, std::vector<OpenMM::RealVec>& forces, std::vector<Vec3>& torque ) const;
 
@@ -495,6 +496,7 @@ private:
                                                const std::vector<int>& multipoleAtomZs,
                                                const std::vector<int>& multipoleAtomXs,
                                                const std::vector<int>& multipoleAtomYs,
+                                               int polarizationType,
                                                std::vector<OpenMM::RealVec>& forces ) const;
 
     /**---------------------------------------------------------------------------------------
@@ -526,6 +528,7 @@ private:
                                                 const std::vector<int>& multipoleAtomZs,
                                                 const std::vector<int>& multipoleAtomXs,
                                                 const std::vector<int>& multipoleAtomYs,
+                                                int polarizationType,
                                                 std::vector<OpenMM::RealVec>& forces );
 
 };