diff --git a/doc/src/fix_ave_grid.rst b/doc/src/fix_ave_grid.rst index 5760bb4508..da0d1f957b 100644 --- a/doc/src/fix_ave_grid.rst +++ b/doc/src/fix_ave_grid.rst @@ -129,12 +129,12 @@ overlays the simulation box. For 2d simulations, *Nz* must be 1. The very coarse compared to the particle count, or very fine. If one or more of the values = 1, then bins are 2d planes or 1d slices of the simulation domain. Note that if the total number of grid cells is -small, it may be more efficient to use the doc:`fix ave/chunk +small, it may be more efficient to use the :doc:`fix ave/chunk ` command which can treat a grid defined by the :doc:`compute chunk/atom ` command as a global grid where each processor owns a copy of all the grid cells. If *Nx* = *Ny* = *Nz* = 1 is used, the same calculation would be more -efficiently performed by the doc:`fix ave/atom ` +efficiently performed by the :doc:`fix ave/atom ` command. If the simulation box size or shape changes during a simulation, the diff --git a/doc/src/fix_polarize.rst b/doc/src/fix_polarize.rst index dd8a0c7d75..f35e34bc67 100644 --- a/doc/src/fix_polarize.rst +++ b/doc/src/fix_polarize.rst @@ -130,10 +130,23 @@ Fix *polarize/functional* employs the energy functional variation approach as described in :ref:`(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) `. - Restart, fix_modify, output, run start/stop, minimize info """""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/doc/src/variable.rst b/doc/src/variable.rst index 67600af62d..a9d21bab9b 100644 --- a/doc/src/variable.rst +++ b/doc/src/variable.rst @@ -1441,7 +1441,7 @@ timestep that the variable needs the tallies. An input script can also request variables be evaluated before or after or in between runs, e.g. by including them in a :doc:`print ` command. -LAMMPS keeps track of all of this as it performs a doc:`run ` or +LAMMPS keeps track of all of this as it performs a :doc:`run ` or :doc:`minimize ` simulation, as well as in between simulations. An error will be generated if you attempt to evaluate a variable when LAMMPS knows it cannot produce accurate values. For @@ -1453,7 +1453,7 @@ command, then an error will occur. However, there are two special cases to be aware when a variable requires invocation of a compute (directly or indirectly). The first -is if the variable is evaluated before the first doc:`run ` or +is if the variable is evaluated before the first :doc:`run ` or :doc:`minimize ` command in the input script. In this case, LAMMPS will generate an error. This is because many computes require initializations which have not yet taken place. One example is the @@ -1463,10 +1463,10 @@ energy or virial quantities; these values are not tallied until the first simulation begins. The second special case is when a variable that depends on a compute -is evaluated in between doc:`run ` or :doc:`minimize ` +is evaluated in between :doc:`run ` or :doc:`minimize ` commands. It is possible for other input script commands issued following the previous run, but before the variable is evaluated, to -change the system. For example, the doc:`delete_atoms ` +change the system. For example, the :doc:`delete_atoms ` command could be used to remove atoms. Since the compute will not re-initialize itself until the next simulation or it may depend on energy/virial computations performed before the system was changed, it diff --git a/src/DIELECTRIC/fix_polarize_bem_gmres.cpp b/src/DIELECTRIC/fix_polarize_bem_gmres.cpp index b11ac7e482..3945b76641 100644 --- a/src/DIELECTRIC/fix_polarize_bem_gmres.cpp +++ b/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; - } } /* ---------------------------------------------------------------------- */ diff --git a/src/DIELECTRIC/fix_polarize_bem_icc.cpp b/src/DIELECTRIC/fix_polarize_bem_icc.cpp index 697d45f416..4dc7dd2669 100644 --- a/src/DIELECTRIC/fix_polarize_bem_icc.cpp +++ b/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]; } } } diff --git a/src/DIELECTRIC/fix_polarize_functional.cpp b/src/DIELECTRIC/fix_polarize_functional.cpp index a0dc034818..ab039f9c48 100644 --- a/src/DIELECTRIC/fix_polarize_functional.cpp +++ b/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 diff --git a/src/DIELECTRIC/pppm_dielectric.cpp b/src/DIELECTRIC/pppm_dielectric.cpp index b9ae1d65ff..e308cb0826 100644 --- a/src/DIELECTRIC/pppm_dielectric.cpp +++ b/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]); } }