EW3DC for non-neutral systems

Add the corrections to force and energy according to Ballenegger, Arnold, and Cerdà, J. Chem. Phys. 131, 094107 2009 (http://dx.doi.org/10.1063/1.3216473).
PicoCentauri · Sep 27, 2016 · 0f94384 · 0f94384
1 parent 9de0eca
commit 0f94384
Showing 1 changed file with 28 additions and 2 deletions.
diff --git a/src/gromacs/ewald/long-range-correction.cpp b/src/gromacs/ewald/long-range-correction.cpp
@@ -328,6 +328,10 @@ void ewald_LRcorrection(int numAtomsLocal,
                 for (j = 0; (j < DIM); j++)
                 {
                     f[i][j] -= dipcorrA[j]*chargeA[i];
+                    if (ewald_geometry == eewg3DC && j == ZZ)
+                    {
+                      f[i][j] += 2*dipole_coeff*fr->qsum[0]*chargeA[i]*x[i][j];
+                    }
                 }
             }
         }
@@ -476,6 +480,11 @@ void ewald_LRcorrection(int numAtomsLocal,
                 {
                     f[i][j] -= L1_q*dipcorrA[j]*chargeA[i]
                         + lambda_q*dipcorrB[j]*chargeB[i];
+                    if (ewald_geometry == eewg3DC && j == ZZ)
+                    {
+                      f[i][j] += 2*dipole_coeff*(L1_q*fr->qsum[0]*chargeA[i]
+                        + lambda_q*fr->qsum[1]*chargeB[i])*x[i][j];
+                    }
                 }
             }
         }
@@ -509,8 +518,9 @@ void ewald_LRcorrection(int numAtomsLocal,
                 }
             }
 
-            /* Apply surface dipole correction:
-             * correction = dipole_coeff * (dipole)^2
+             /* Apply surface and charged surface dipole correction:
+             * correction = dipole_coeff * ( (dipole)^2
+             *               - qsum*sum_i q_i z_i^2 - qsum^2 * box_z^2 / 12 )
              */
             if (dipole_coeff != 0)
             {
@@ -521,6 +531,22 @@ void ewald_LRcorrection(int numAtomsLocal,
                 else if (ewald_geometry == eewg3DC)
                 {
                     Vdipole[q] = dipole_coeff*mutot[q][ZZ]*mutot[q][ZZ];
+                    real qzs2 = 0;
+                    for (i = start; (i < end); i++)
+                    {
+                        if (i < numAtomsLocal)
+                        {
+                          if (q == 0)
+                          {
+                            qzs2 += chargeA[i]*x[i][ZZ]*x[i][ZZ];
+                          }
+                          else if (q == 1)
+                          {
+                            qzs2 += chargeB[i]*x[i][ZZ]*x[i][ZZ];
+                          }
+                        }
+                      }
+                      Vdipole[q] -= dipole_coeff*(fr->qsum[q]*qzs2+fr->qsum[q]*fr->qsum[q]*box[ZZ][ZZ]*box[ZZ][ZZ]/12.);
                 }
             }
         }