Neutralizing plasma self energy for AMOEBA (#4920)

plugins/amoeba/platforms/common/src/AmoebaCommonKernels.cpp

@@ -233,11 +233,13 @@ void CommonCalcAmoebaMultipoleForceKernel::initialize(const System& system, cons
     vector<float> localDipolesVec;
     vector<float> localQuadrupolesVec;
     vector<mm_int4> multipoleParticlesVec;
+    totalCharge = 0.0;
     for (int i = 0; i < numMultipoles; i++) {
         double charge, thole, damping, polarity;
         int axisType, atomX, atomY, atomZ;
         vector<double> dipole, quadrupole;
         force.getMultipoleParameters(i, charge, dipole, quadrupole, axisType, atomZ, atomX, atomY, thole, damping, polarity);
+        totalCharge += charge;
         if (cc.getUseDoublePrecision())
             posqd[i] = mm_double4(0, 0, 0, charge);
         else
@@ -1228,7 +1230,17 @@ double CommonCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool
     
     cc.getPosq().copyTo(lastPositions);
     multipolesAreValid = true;
-    return 0.0;
+
+    // Correction for the neutralizing plasma.
+
+    if (usePME) {
+        Vec3 a, b, c;
+        cc.getPeriodicBoxVectors(a, b, c);
+        double volume = a[0] * b[1] * c[2];
+        return -totalCharge*totalCharge/(8*EPSILON0*volume*pmeAlpha*pmeAlpha);
+    }
+    else
+        return 0.0;
 }
 
 void CommonCalcAmoebaMultipoleForceKernel::computeInducedField() {
@@ -1664,11 +1676,13 @@ void CommonCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl&
     vector<float> localDipolesVec;
     vector<float> localQuadrupolesVec;
     vector<mm_int4> multipoleParticlesVec;
+    totalCharge = 0.0;
     for (int i = 0; i < force.getNumMultipoles(); i++) {
         double charge, thole, damping, polarity;
         int axisType, atomX, atomY, atomZ;
         vector<double> dipole, quadrupole;
         force.getMultipoleParameters(i, charge, dipole, quadrupole, axisType, atomZ, atomX, atomY, thole, damping, polarity);
+        totalCharge += charge;
         if (cc.getUseDoublePrecision())
             posqd[i].w = charge;
         else
@@ -2424,12 +2438,14 @@ void CommonCalcHippoNonbondedForceKernel::initialize(const System& system, const
     vector<double> localDipolesVec, localQuadrupolesVec;
     vector<mm_int4> multipoleParticlesVec;
     vector<vector<int> > exclusions(numParticles);
+    totalCharge = 0.0;
     for (int i = 0; i < numParticles; i++) {
         double charge, coreCharge, alpha, epsilon, damping, c6, pauliK, pauliQ, pauliAlpha, polarizability;
         int axisType, atomX, atomY, atomZ;
         vector<double> dipole, quadrupole;
         force.getParticleParameters(i, charge, dipole, quadrupole, coreCharge, alpha, epsilon, damping, c6, pauliK, pauliQ, pauliAlpha,
                                     polarizability, axisType, atomZ, atomX, atomY);
+        totalCharge += charge;
         coreChargeVec.push_back(coreCharge);
         valenceChargeVec.push_back(charge-coreCharge);
         alphaVec.push_back(alpha);
@@ -3378,7 +3394,17 @@ double CommonCalcHippoNonbondedForceKernel::execute(ContextImpl& context, bool i
     
     cc.getPosq().copyTo(lastPositions);
     multipolesAreValid = true;
-    return 0.0;
+
+    // Correction for the neutralizing plasma.
+
+    if (usePME) {
+        Vec3 a, b, c;
+        cc.getPeriodicBoxVectors(a, b, c);
+        double volume = a[0] * b[1] * c[2];
+        return -totalCharge*totalCharge/(8*EPSILON0*volume*pmeAlpha*pmeAlpha);
+    }
+    else
+        return 0.0;
 }
 
 void CommonCalcHippoNonbondedForceKernel::computeInducedField(int optOrder) {
@@ -3513,12 +3539,14 @@ void CommonCalcHippoNonbondedForceKernel::copyParametersToContext(ContextImpl& c
     vector<double> coreChargeVec, valenceChargeVec, alphaVec, epsilonVec, dampingVec, c6Vec, pauliKVec, pauliQVec, pauliAlphaVec, polarizabilityVec;
     vector<double> localDipolesVec, localQuadrupolesVec;
     vector<mm_int4> multipoleParticlesVec;
+    totalCharge = 0.0;
     for (int i = 0; i < numParticles; i++) {
         double charge, coreCharge, alpha, epsilon, damping, c6, pauliK, pauliQ, pauliAlpha, polarizability;
         int axisType, atomX, atomY, atomZ;
         vector<double> dipole, quadrupole;
         force.getParticleParameters(i, charge, dipole, quadrupole, coreCharge, alpha, epsilon, damping, c6, pauliK, pauliQ, pauliAlpha,
                                     polarizability, axisType, atomZ, atomX, atomY);
+        totalCharge += charge;
         coreChargeVec.push_back(coreCharge);
         valenceChargeVec.push_back(charge-coreCharge);
         alphaVec.push_back(alpha);

plugins/amoeba/platforms/common/src/AmoebaCommonKernels.h

@@ -168,7 +168,7 @@ protected:
     int numMultipoles, maxInducedIterations, maxExtrapolationOrder;
     int fixedFieldThreads, inducedFieldThreads, electrostaticsThreads;
     int gridSizeX, gridSizeY, gridSizeZ;
-    double pmeAlpha, inducedEpsilon;
+    double pmeAlpha, inducedEpsilon, totalCharge;
     bool usePME, hasQuadrupoles, hasInitializedScaleFactors, multipolesAreValid, hasCreatedEvent;
     AmoebaMultipoleForce::PolarizationType polarizationType;
     ComputeContext& cc;
@@ -507,7 +507,7 @@ protected:
     int numParticles, maxExtrapolationOrder, maxTiles, fieldThreadBlockSize;
     int gridSizeX, gridSizeY, gridSizeZ;
     int dispersionGridSizeX, dispersionGridSizeY, dispersionGridSizeZ;
-    double pmeAlpha, dpmeAlpha, cutoff;
+    double pmeAlpha, dpmeAlpha, cutoff, totalCharge;
     bool usePME, hasInitializedKernels, multipolesAreValid;
     std::vector<double> extrapolationCoefficients;
     ComputeContext& cc;

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

 
-/* Portions copyright (c) 2006-2022 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2025 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -2503,9 +2503,11 @@ double AmoebaReferencePmeHippoNonbondedForce::calculatePmeSelfEnergy(const vecto
     double dii = 0.0;
     double qii = 0.0;
     double c6ii = 0.0;
+    double totalCharge = 0.0;
     for (int i = 0; i < _numParticles; i++) {
         const MultipoleParticleData& particleI = particleData[i];
         double charge = particleI.coreCharge + particleI.valenceCharge;
+        totalCharge += charge;
         cii += charge*charge;
         dii += particleI.dipole.dot(particleI.dipole);
         qii += particleI.quadrupole[QXX]*particleI.quadrupole[QXX] +
@@ -2519,7 +2521,11 @@ double AmoebaReferencePmeHippoNonbondedForce::calculatePmeSelfEnergy(const vecto
     double term = 2*_alphaEwald*_alphaEwald;
     double fterm = -_electric*_alphaEwald/SQRT_PI;
     double alpha3 = _dalphaEwald*_dalphaEwald*_dalphaEwald;
-    return fterm*(cii + term*(dii/3+2*term*qii/5)) + alpha3*alpha3*c6ii/12;
+    double energy = fterm*(cii + term*(dii/3+2*term*qii/5)) + alpha3*alpha3*c6ii/12;
+    // Correction for the neutralizing plasma.
+    double volume = _periodicBoxVectors[0][0] * _periodicBoxVectors[1][1] * _periodicBoxVectors[2][2];
+    energy -= totalCharge*totalCharge/(8*EPSILON0*volume*_alphaEwald*_alphaEwald);
+    return energy;
 }
 
 void AmoebaReferencePmeHippoNonbondedForce::calculatePmeSelfTorque(const vector<MultipoleParticleData>& particleData,

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

 
-/* Portions copyright (c) 2006-222 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2025 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -6370,10 +6370,12 @@ double AmoebaReferencePmeMultipoleForce::calculatePmeSelfEnergy(const vector<Mul
     double cii = 0.0;
     double dii = 0.0;
     double qii = 0.0;
+    double totalCharge = 0.0;
     for (unsigned int ii = 0; ii < _numParticles; ii++) {
 
         const MultipoleParticleData& particleI = particleData[ii];
 
+        totalCharge += particleI.charge;
         cii += particleI.charge*particleI.charge;
 
         Vec3 dipole(particleI.sphericalDipole[1], particleI.sphericalDipole[2], particleI.sphericalDipole[0]);
@@ -6389,6 +6391,9 @@ double AmoebaReferencePmeMultipoleForce::calculatePmeSelfEnergy(const vector<Mul
     double a2 = _alphaEwald * _alphaEwald;
     double a4 = a2*a2;
     double energy = prefac*(cii + twoThirds*a2*dii + fourOverFifteen*a4*qii);
+    // Correction for the neutralizing plasma.
+    double volume = _periodicBoxVectors[0][0] * _periodicBoxVectors[1][1] * _periodicBoxVectors[2][2];
+    energy -= totalCharge*totalCharge/(8*EPSILON0*volume*_alphaEwald*_alphaEwald);
     return energy;
 }
 

plugins/amoeba/tests/TestAmoebaMultipoleForce.h

@@ -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-2025 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -3261,6 +3261,48 @@ void testZOnly() {
     ASSERT_EQUAL_TOL(state2.getPotentialEnergy(), state3.getPotentialEnergy()+norm*delta, 1e-3)
 }
 
+void testNeutralizingPlasmaCorrection() {
+    // Verify that the energy of a system with nonzero charge doesn't depend on alpha.
+
+    System system;
+    AmoebaMultipoleForce* force = new AmoebaMultipoleForce();
+    force->setNonbondedMethod(AmoebaMultipoleForce::PME);
+    system.addForce(force);
+    vector<double> d(3, -.0), q(9, 0.0);
+    for (int i = 0; i < 2; i++) {
+        system.addParticle(1.0);
+        force->addMultipole(1.0, d, q, AmoebaMultipoleForce::NoAxisType, 0, 0, 0, 0.39, 0.33, 0.001);
+    }
+    vector<Vec3> positions(2);
+    positions[0] = Vec3();
+    positions[1] = Vec3(0.3, 0.4, 0.0);
+
+    // Compute the energy.
+
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    double energy1 = context.getState(State::Energy).getPotentialEnergy();
+
+    // Change the cutoff distance, which will change alpha, and see if the energy is the same.
+
+    force->setCutoffDistance(0.7);
+    context.reinitialize(true);
+    double energy2 = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(energy1, energy2, 1e-4);
+
+    // Try changing a particle charge with updateParametersInContext() and make sure the
+    // energy changes by the correct amount.
+
+    force->setMultipoleParameters(0, 2.0, d, q, AmoebaMultipoleForce::NoAxisType, 0, 0, 0, 0.39, 0.33, 0.001);
+    force->updateParametersInContext(context);
+    double energy3 = context.getState(State::Energy).getPotentialEnergy();
+    force->setCutoffDistance(1.0);
+    context.reinitialize(true);
+    double energy4 = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(energy3, energy4, 1e-4);
+}
+
 void setupKernels(int argc, char* argv[]);
 void runPlatformTests();
 
@@ -3318,6 +3360,7 @@ int main(int argc, char* argv[]) {
         // test the ZOnly axis type.
         
         testZOnly();
+        testNeutralizingPlasmaCorrection();
 
         runPlatformTests();
     }

plugins/amoeba/tests/TestHippoNonbondedForce.h

@@ -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) 2019 Stanford University and the Authors.           *
+ * Portions copyright (c) 2019-2025 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -1542,6 +1542,48 @@ void testChangingParameters() {
         ASSERT_EQUAL_VEC(state2.getForces()[i], state3.getForces()[i], 1e-5);
 }
 
+void testNeutralizingPlasmaCorrection() {
+    // Verify that the energy of a system with nonzero charge doesn't depend on alpha.
+
+    System system;
+    HippoNonbondedForce* force = new HippoNonbondedForce();
+    force->setNonbondedMethod(HippoNonbondedForce::PME);
+    system.addForce(force);
+    vector<double> d(3, 0.0), q(9, 0.0);
+    for (int i = 0; i < 2; i++) {
+        system.addParticle(1.0);
+        force->addParticle(1.0, d, q, 6.0, 50.0, 5000.0, 400.0, 0.04, 1.5, -2.4233, 40.0, 0.001, HippoNonbondedForce::NoAxisType, -1, -1, -1);
+    }
+    vector<Vec3> positions(2);
+    positions[0] = Vec3();
+    positions[1] = Vec3(0.3, 0.4, 0.0);
+
+    // Compute the energy.
+
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    double energy1 = context.getState(State::Energy).getPotentialEnergy();
+
+    // Change the cutoff distance, which will change alpha, and see if the energy is the same.
+
+    force->setCutoffDistance(0.7);
+    context.reinitialize(true);
+    double energy2 = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(energy1, energy2, 1e-4);
+
+    // Try changing a particle charge with updateParametersInContext() and make sure the
+    // energy changes by the correct amount.
+
+    force->setParticleParameters(0, 2.0, d, q, 6.0, 50.0, 5000.0, 400.0, 0.04, 1.5, -2.4233, 40.0, 0.001, HippoNonbondedForce::NoAxisType, -1, -1, -1);
+    force->updateParametersInContext(context);
+    double energy3 = context.getState(State::Energy).getPotentialEnergy();
+    force->setCutoffDistance(1.0);
+    context.reinitialize(true);
+    double energy4 = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(energy3, energy4, 1e-4);
+}
+
 int main(int argc, char* argv[]) {
     try {
         setupKernels(argc, argv);
@@ -1562,6 +1604,7 @@ int main(int argc, char* argv[]) {
         testWaterDimer();
         testWaterBox();
         testChangingParameters();
+        testNeutralizingPlasmaCorrection();
     }
     catch (const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;