Support GK with extrapolated polarization

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

@@ -1145,12 +1145,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::Extrapolated) {
+    else if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated)
         convergeInduceDipolesByExtrapolation(particleData, updateInducedDipoleField);
-    }
-
 }
 
 RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPairIxn(const MultipoleParticleData& particleI,
@@ -2555,7 +2553,10 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipoles(c
     updateInducedDipoleField.push_back(UpdateInducedDipoleFieldStruct(_gkField, _inducedDipoleS, _ptDipoleDS, _ptDipoleFieldGradientDS));
     updateInducedDipoleField.push_back(UpdateInducedDipoleFieldStruct(gkFieldPolar, _inducedDipolePolarS, _ptDipolePS, _ptDipoleFieldGradientPS));
 
-    convergeInduceDipolesByDIIS(particleData, updateInducedDipoleField);
+    if (getPolarizationType() == AmoebaReferenceMultipoleForce::Mutual)
+        convergeInduceDipolesByDIIS(particleData, updateInducedDipoleField);
+    else if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated)
+        convergeInduceDipolesByExtrapolation(particleData, updateInducedDipoleField);
 }
 
 RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodPairIxn(const MultipoleParticleData& particleI,
@@ -3810,7 +3811,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodPa
 
     // mutual polarization electrostatic solvation free energy gradient
 
-    if (getPolarizationType() == AmoebaReferenceMultipoleForce::Mutual) {
+    if (getPolarizationType() == AmoebaReferenceMultipoleForce::Mutual || getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
 
         dpdx = dpdx - 0.5 *
                            (_inducedDipoleS[iIndex][0]*(_inducedDipolePolarS[jIndex][0]*gux5+_inducedDipolePolarS[jIndex][1]*gux6+_inducedDipolePolarS[jIndex][2]*gux7)
@@ -4029,6 +4030,55 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateElectrosta
     }
     energy += (_electric/_dielectric)*eDiffEnergy;
 
+    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 = _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]
+                                                + _ptDipoleD[l][i][2]*_ptDipoleFieldGradientP[m][6*i+4]);
+                    forces[i][1] -= 0.5*p*prefac*(_ptDipoleD[l][i][0]*_ptDipoleFieldGradientP[m][6*i+3]
+                                                + _ptDipoleD[l][i][1]*_ptDipoleFieldGradientP[m][6*i+1]
+                                                + _ptDipoleD[l][i][2]*_ptDipoleFieldGradientP[m][6*i+5]);
+                    forces[i][2] -= 0.5*p*prefac*(_ptDipoleD[l][i][0]*_ptDipoleFieldGradientP[m][6*i+4]
+                                                + _ptDipoleD[l][i][1]*_ptDipoleFieldGradientP[m][6*i+5]
+                                                + _ptDipoleD[l][i][2]*_ptDipoleFieldGradientP[m][6*i+2]);
+                    forces[i][0] -= 0.5*p*prefac*(_ptDipoleP[l][i][0]*_ptDipoleFieldGradientD[m][6*i+0]
+                                                + _ptDipoleP[l][i][1]*_ptDipoleFieldGradientD[m][6*i+3]
+                                                + _ptDipoleP[l][i][2]*_ptDipoleFieldGradientD[m][6*i+4]);
+                    forces[i][1] -= 0.5*p*prefac*(_ptDipoleP[l][i][0]*_ptDipoleFieldGradientD[m][6*i+3]
+                                                + _ptDipoleP[l][i][1]*_ptDipoleFieldGradientD[m][6*i+1]
+                                                + _ptDipoleP[l][i][2]*_ptDipoleFieldGradientD[m][6*i+5]);
+                    forces[i][2] -= 0.5*p*prefac*(_ptDipoleP[l][i][0]*_ptDipoleFieldGradientD[m][6*i+4]
+                                                + _ptDipoleP[l][i][1]*_ptDipoleFieldGradientD[m][6*i+5]
+                                                + _ptDipoleP[l][i][2]*_ptDipoleFieldGradientD[m][6*i+2]);
+                    forces[i][0] += 0.5*p*prefac*(_ptDipoleDS[l][i][0]*_ptDipoleFieldGradientPS[m][6*i+0]
+                                                + _ptDipoleDS[l][i][1]*_ptDipoleFieldGradientPS[m][6*i+3]
+                                                + _ptDipoleDS[l][i][2]*_ptDipoleFieldGradientPS[m][6*i+4]);
+                    forces[i][1] += 0.5*p*prefac*(_ptDipoleDS[l][i][0]*_ptDipoleFieldGradientPS[m][6*i+3]
+                                                + _ptDipoleDS[l][i][1]*_ptDipoleFieldGradientPS[m][6*i+1]
+                                                + _ptDipoleDS[l][i][2]*_ptDipoleFieldGradientPS[m][6*i+5]);
+                    forces[i][2] += 0.5*p*prefac*(_ptDipoleDS[l][i][0]*_ptDipoleFieldGradientPS[m][6*i+4]
+                                                + _ptDipoleDS[l][i][1]*_ptDipoleFieldGradientPS[m][6*i+5]
+                                                + _ptDipoleDS[l][i][2]*_ptDipoleFieldGradientPS[m][6*i+2]);
+                    forces[i][0] += 0.5*p*prefac*(_ptDipolePS[l][i][0]*_ptDipoleFieldGradientDS[m][6*i+0]
+                                                + _ptDipolePS[l][i][1]*_ptDipoleFieldGradientDS[m][6*i+3]
+                                                + _ptDipolePS[l][i][2]*_ptDipoleFieldGradientDS[m][6*i+4]);
+                    forces[i][1] += 0.5*p*prefac*(_ptDipolePS[l][i][0]*_ptDipoleFieldGradientDS[m][6*i+3]
+                                                + _ptDipolePS[l][i][1]*_ptDipoleFieldGradientDS[m][6*i+1]
+                                                + _ptDipolePS[l][i][2]*_ptDipoleFieldGradientDS[m][6*i+5]);
+                    forces[i][2] += 0.5*p*prefac*(_ptDipolePS[l][i][0]*_ptDipoleFieldGradientDS[m][6*i+4]
+                                                + _ptDipolePS[l][i][1]*_ptDipoleFieldGradientDS[m][6*i+5]
+                                                + _ptDipolePS[l][i][2]*_ptDipoleFieldGradientDS[m][6*i+2]);
+                }
+            }
+        }
+    }
+
     return energy;
 }
 
@@ -4445,7 +4495,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
 
     // correction to convert mutual to direct polarization force
 
-    if (getPolarizationType() == AmoebaReferenceMultipoleForce::Direct) {
+    if (getPolarizationType() != AmoebaReferenceMultipoleForce::Mutual) {
         RealOpenMM gfd      = 0.5*(rr5*scip2*scale3i - rr7*(scip3*sci4+sci3*scip4)*scale5i);
         RealOpenMM fdir1    = gfd*xr + 0.5*rr5*scale5i* (sci4*_inducedDipolePolarS[iIndex][0]+scip4*_inducedDipoleS[iIndex][0] + sci3*_inducedDipolePolarS[jIndex][0]+scip3*_inducedDipoleS[jIndex][0]);
         RealOpenMM fdir2    = gfd*yr + 0.5*rr5*scale5i* (sci4*_inducedDipolePolarS[iIndex][1]+scip4*_inducedDipoleS[iIndex][1] + sci3*_inducedDipolePolarS[jIndex][1]+scip3*_inducedDipoleS[jIndex][1]);
@@ -4708,7 +4758,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
 
     // correction to convert mutual to direct polarization force
 
-    if (getPolarizationType() == AmoebaReferenceMultipoleForce::Direct) {
+    if (getPolarizationType() != AmoebaReferenceMultipoleForce::Mutual) {
 
         RealOpenMM gfd    = 0.5*(rr5*scip2*scale3i- rr7*(scip3*sci4+sci3*scip4)*scale5i);
         RealOpenMM fdir1  = gfd*xr + 0.5*rr5*scale5i* (sci4*_inducedDipolePolar[iIndex][0]+scip4*_inducedDipole[iIndex][0] + sci3*_inducedDipolePolar[jIndex][0]+scip3*_inducedDipole[jIndex][0]);

plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaGeneralizedKirkwoodForce.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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -50,11 +50,36 @@
 using namespace OpenMM;
+using namespace std;
 extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
 const double TOL = 1e-4;
+static void checkFiniteDifferences(vector<Vec3> forces, Context &context, vector<Vec3> positions)
+{
+    // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
+
+    double norm = 0.0;
+    for (int i = 0; i < (int) forces.size(); ++i)
+        norm += forces[i].dot(forces[i]);
+    norm = std::sqrt(norm);
+    const double stepSize = 1e-3;
+    double step = 0.5*stepSize/norm;
+    vector<Vec3> positions2(forces.size()), positions3(forces.size());
+    for (int i = 0; i < (int) positions.size(); ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = forces[i];
+        positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+        positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
+    }
+    context.setPositions(positions2);
+    State state2 = context.getState(State::Energy);
+    context.setPositions(positions3);
+    State state3 = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
+}
+
 // setup for 2 ammonia molecules
 static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce,
@@ -291,6 +316,10 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>&
     State state                      = context.getState(State::Forces | State::Energy);
     forces                           = state.getForces();
     energy                           = state.getPotentialEnergy();
+    
+    // Check that the forces and energy are consistent.
+    
+    checkFiniteDifferences(forces, context, positions);
 }
 // setup for villin
@@ -6974,6 +7003,10 @@ static void setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Polariz
     State state                      = context.getState(State::Forces | State::Energy);
     forces                           = state.getForces();
     energy                           = state.getPotentialEnergy();
+    
+    // Check that the forces and energy are consistent.
+    
+    checkFiniteDifferences(forces, context, positions);
 }
 // compare forces and energies 
@@ -7052,6 +7085,25 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization() {
     compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
 }
+static void testGeneralizedKirkwoodAmmoniaExtrapolatedPolarization() {
+
+    std::string testName      = "testGeneralizedKirkwoodAmmoniaExtrapolatedPolarization";
+
+    int numberOfParticles     = 8;
+    std::vector<Vec3> forces;
+    double energy;
+
+    System system;
+    AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce  = new AmoebaGeneralizedKirkwoodForce();
+    setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+
+    // We don't have reference values for this case, but at least check that force and energy are consistent.
+
+    getForcesEnergyMultipoleAmmonia(context, forces, energy);
+}
+
 // test GK mutual polarization for system comprised of two ammonia molecules
 static void testGeneralizedKirkwoodAmmoniaMutualPolarization() {
@@ -7770,6 +7822,19 @@ static void testGeneralizedKirkwoodVillinDirectPolarization() {
     compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
 }
+static void testGeneralizedKirkwoodVillinExtrapolatedPolarization() {
+
+    std::string testName      = "testGeneralizedKirkwoodVillinExtrapolatedPolarization";
+
+    int numberOfParticles     = 596;
+    std::vector<Vec3> forces;
+    double energy;
+
+    // We don't have reference values for this case, but at least check that force and energy are consistent.
+    
+    setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Extrapolated, 0, forces, energy);
+}
+
 // test GK mutual polarization for villin system
 static void testGeneralizedKirkwoodVillinMutualPolarization() {
@@ -8405,8 +8470,10 @@ int main(int numberOfArguments, char* argv[]) {
         testGeneralizedKirkwoodAmmoniaMutualPolarization();
         testGeneralizedKirkwoodAmmoniaDirectPolarization();
+        testGeneralizedKirkwoodAmmoniaExtrapolatedPolarization();
         testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm();
         testGeneralizedKirkwoodVillinDirectPolarization();
+        testGeneralizedKirkwoodVillinExtrapolatedPolarization();
         testGeneralizedKirkwoodVillinMutualPolarization();
     }

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -247,7 +247,7 @@ UNITS = {
 ("AmoebaMultipoleForce",                 "getPmeBSplineOrder")                            :  ( None,()),
 ("AmoebaMultipoleForce",                 "getMutualInducedMaxIterations")                 :  ( None, ()),
 ("AmoebaMultipoleForce",                 "getMutualInducedTargetEpsilon")                 :  ( None, ()),
-("AmoebaMultipoleForce",                 "getOPTCoefficients")                            :  ( None, ()),
+("AmoebaMultipoleForce",                 "getExtrapolationCoefficients")                            :  ( None, ()),
 ("AmoebaMultipoleForce",                 "getEwaldErrorTolerance")                        :  ( None, ()),
 ("AmoebaMultipoleForce",                 "getPmeGridDimensions")                          :  ( None,()),