Enforced net induced charges to be zero

src/DIELECTRIC/fix_polarize_bem_gmres.cpp

@@ -41,6 +41,7 @@
 #include "comm.h"
 #include "error.h"
 #include "force.h"
+#include "group.h"
 #include "kspace.h"
 #include "math_const.h"
 #include "memory.h"
@@ -347,6 +348,7 @@ void FixPolarizeBEMGMRES::compute_induced_charges()
   double *ed = atom->ed;
   double *em = atom->em;
   double *epsilon = atom->epsilon;
+  int *mask = atom->mask;
   int nlocal = atom->nlocal;
   int eflag = 0;
   int vflag = 0;
@@ -438,6 +440,26 @@ void FixPolarizeBEMGMRES::compute_induced_charges()
   comm->forward_comm(this);
 
   if (first) first = 0;
+
+  // ensure sum of all induced charges being zero
+
+  double tmp = 0;
+  int ncount = group->count(igroup);
+  for (int i = 0; i < nlocal; i++) {
+    if (!(mask[i] & groupbit)) continue;
+
+    double q_bound = q[i] - q_real[i];
+    tmp += q_bound;
+  }
+
+  double sum = 0;
+  MPI_Allreduce(&tmp, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
+  double qboundave = sum/(double)ncount;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (!(mask[i] & groupbit)) continue;
+    q[i] -=  qboundave;
+  }
 }
 
 /* ---------------------------------------------------------------------- */

src/DIELECTRIC/fix_polarize_bem_icc.cpp

@@ -351,6 +351,26 @@ void FixPolarizeBEMICC::compute_induced_charges()
   }
 
   iterations = itr;
+
+  // ensure sum of all induced charges being zero
+
+  double tmp = 0;
+  int ncount = group->count(igroup);
+  for (int i = 0; i < nlocal; i++) {
+    if (!(mask[i] & groupbit)) continue;
+
+    double q_bound = q[i] - q_real[i];
+    tmp += q_bound;
+  }
+
+  double sum = 0;
+  MPI_Allreduce(&tmp, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
+  double qboundave = sum/(double)ncount;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (!(mask[i] & groupbit)) continue;
+    q[i] -=  qboundave;
+  }
 }
 
 /* ---------------------------------------------------------------------- */

src/DIELECTRIC/fix_polarize_functional.cpp

@@ -377,11 +377,23 @@ void FixPolarizeFunctional::update_induced_charges()
 
   double *q = atom->q;
   int nlocal = atom->nlocal;
+  double tmp = 0;
+
 
   for (int i = 0; i < nlocal; i++) {
     if (induced_charge_idx[i] < 0) continue;
     int idx = induced_charge_idx[i];
     q[i] = -induced_charges[idx] / (4 * MY_PI);
+    tmp += q[i];
+  }
+
+  double sum = 0;
+  MPI_Allreduce(&tmp, &sum, 1, MPI_DOUBLE, MPI_SUM, world);
+  double qboundave = sum/(double)num_induced_charges;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (induced_charge_idx[i] < 0) continue;
+    q[i] -=  qboundave;
   }
 
   // revert to scaled charges to calculate forces