Merge pull request #3845 from ndtrung81/dielectric-updates2

Update and bugfix for the DIELECTRIC package

doc/src/fix_polarize.rst

@@ -130,10 +130,23 @@ Fix *polarize/functional* employs the energy functional variation
 approach as described in :ref:`(Jadhao) <Jadhao>` to solve
 :math:`\sigma_b`.
 
+The induced charges computed by these fixes are stored in the *q_scaled* field,
+and can be accessed as in the following example:
+
+.. code-block:: LAMMPS
+
+  compute qs all property/atom q_scaled
+  dump 1 all custom 1000 all.txt id type q x y z c_qs
+
+Note that the *q* field is the regular atom charges, which do not change
+during the simulation. For interface particles, *q_scaled* is the sum
+of the real charge, divided by the local dielectric constant *epsilon*,
+and their induced charges. For non-interface particles, *q_scaled* is
+the real charge, divided by the local dielectric constant *epsilon*.
+
 More details on the implementation of these fixes and their recommended
 use are described in :ref:`(NguyenTD) <NguyenTD>`.
 
-
 Restart, fix_modify, output, run start/stop, minimize info
 """"""""""""""""""""""""""""""""""""""""""""""""""""""""""
 

src/DIELECTRIC/fix_polarize_bem_gmres.cpp

@@ -445,31 +445,7 @@ void FixPolarizeBEMGMRES::compute_induced_charges()
 
   comm->forward_comm(this);
 
-  // compute the total induced charges of the interface particles
-  // for interface particles: set the charge to be the sum of unscaled (free) charges and induced charges
-
-  double tmp = 0;
-  for (int i = 0; i < nlocal; i++) {
-    if (!(mask[i] & groupbit)) continue;
-
-    double q_bound = q_scaled[i] - q[i];
-    tmp += q_bound;
-    q[i] = q_scaled[i];
-  }
-
   if (first) first = 0;
-
-  // ensure sum of all induced charges being zero
-
-  int ncount = group->count(igroup);
-  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

@@ -361,30 +361,6 @@ void FixPolarizeBEMICC::compute_induced_charges()
   }
 
   iterations = itr;
-
-  // compute the total induced charges of the interface particles
-  // for interface particles: set the charge to be the sum of unscaled (free) charges and induced charges
-
-  double tmp = 0;
-  for (int i = 0; i < nlocal; i++) {
-    if (!(mask[i] & groupbit)) continue;
-
-    double q_bound = q_scaled[i] - q[i];
-    tmp += q_bound;
-    q[i] = q_scaled[i];
-  }
-
-  // ensure sum of all induced charges being zero
-
-  int ncount = group->count(igroup);
-  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;
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -479,6 +455,7 @@ void FixPolarizeBEMICC::set_dielectric_params(double ediff, double emean, double
   double *ed = atom->ed;
   double *em = atom->em;
   double *q = atom->q;
+  double *q_scaled = atom->q_scaled;
   double *epsilon = atom->epsilon;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
@@ -490,6 +467,7 @@ void FixPolarizeBEMICC::set_dielectric_params(double ediff, double emean, double
       if (areai > 0) area[i] = areai;
       if (epsiloni > 0) epsilon[i] = epsiloni;
       if (set_charge) q[i] = qvalue;
+      q_scaled[i] = q[i] / epsilon[i];
     }
   }
 }

src/DIELECTRIC/fix_polarize_functional.cpp

@@ -386,18 +386,7 @@ void FixPolarizeFunctional::update_induced_charges()
   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);
-    q_scaled[i] = q[i] / epsilon[i];
-    tmp += q_scaled[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_scaled[i] -=  qboundave;
+    q_scaled[i] = -induced_charges[idx] / (4 * MY_PI);
   }
 
   // revert to scaled charges to calculate forces

src/DIELECTRIC/pppm_dielectric.cpp

@@ -575,7 +575,7 @@ void PPPMDielectric::slabcorr()
   int nlocal = atom->nlocal;
 
   double dipole = 0.0;
-  for (int i = 0; i < nlocal; i++) dipole += q[i]*x[i][2];
+  for (int i = 0; i < nlocal; i++) dipole += eps[i]*q[i]*x[i][2];
 
   // sum local contributions to get global dipole moment
 
@@ -588,7 +588,7 @@ void PPPMDielectric::slabcorr()
   double dipole_r2 = 0.0;
   if (eflag_atom || fabs(qsum) > SMALL) {
     for (int i = 0; i < nlocal; i++)
-      dipole_r2 += q[i]*x[i][2]*x[i][2];
+      dipole_r2 += eps[i]*q[i]*x[i][2]*x[i][2];
 
     // sum local contributions
 
@@ -621,6 +621,6 @@ void PPPMDielectric::slabcorr()
 
   for (int i = 0; i < nlocal; i++) {
     f[i][2] += ffact * eps[i]*q[i]*(dipole_all - qsum*x[i][2]);
-    efield[i][2] += ffact * eps[i]*(dipole_all - qsum*x[i][2]);
+    efield[i][2] += ffact * (dipole_all - qsum*x[i][2]);
   }
 }