Merge pull request #3551 from akohlmey/compute-efield-atom-wolf

Add new compute efield/wolf/atom command

doc/src/Commands_compute.rst

@@ -57,6 +57,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`dpd/atom <compute_dpd_atom>`
    * :doc:`edpd/temp/atom <compute_edpd_temp_atom>`
    * :doc:`efield/atom <compute_efield_atom>`
+   * :doc:`efield/wolf/atom <compute_efield_wolf_atom>`
    * :doc:`entropy/atom <compute_entropy_atom>`
    * :doc:`erotate/asphere <compute_erotate_asphere>`
    * :doc:`erotate/rigid <compute_erotate_rigid>`

doc/src/compute.rst

@@ -211,6 +211,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`dpd/atom <compute_dpd_atom>` - per-particle values of internal conductive energy, internal mechanical energy, chemical energy, and internal temperature
 * :doc:`edpd/temp/atom <compute_edpd_temp_atom>` - per-atom temperature for each eDPD particle in a group
 * :doc:`efield/atom <compute_efield_atom>` - electric field at each atom
+* :doc:`efield/wolf/atom <compute_efield_wolf_atom>` - electric field at each atom
 * :doc:`entropy/atom <compute_entropy_atom>` - pair entropy fingerprint of each atom
 * :doc:`erotate/asphere <compute_erotate_asphere>` - rotational energy of aspherical particles
 * :doc:`erotate/rigid <compute_erotate_rigid>` - rotational energy of rigid bodies

doc/src/compute_efield_wolf_atom.rst

@@ -0,0 +1,126 @@
+.. index:: compute efield/wolf/atom
+
+compute efield/wolf/atom command
+================================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   compute ID group-ID efield/wolf/atom alpha keyword val
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* efield/atom/wolf = style name of this compute command
+* alpha = damping parameter (inverse distance units)
+* zero or more keyword/value pairs may be appended
+* keyword = *limit* or *cutoff*
+
+  .. parsed-literal::
+
+     *limit* group2-ID = limit computing the electric field contributions to a group (default: all)
+     *cutoff* value = set cutoff for computing contributions to this value (default: maximum cutoff of pair style)
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute 1 all efield/wolf/atom 0.2
+   compute 1 mols efield/wolf/atom 0.25 limit water cutoff 10.0
+
+Description
+"""""""""""
+
+.. versionadded:: TBD
+
+Define a computation that approximates the electric field at each atom in a group.
+
+.. math::
+
+   \vec{E}_i = \frac{\vec{F}coul_i}{q_i} = \sum_{j \neq i} \frac{q_j}{r_{ij}^2} \qquad r < r_c
+
+The electric field at the position of the atom *i* is the coulomb force
+on a unit charge at that point, which is equivalent to dividing the
+Coulomb force by the charge of the individual atom.
+
+In this compute the electric field is approximated as the derivative of
+the potential energy using the Wolf summation method, described in
+:ref:`Wolf <Wolf4>`, given by:
+
+.. math::
+   E_i = \frac{1}{2} \sum_{j \neq i}
+   \frac{q_i q_j {\rm erfc}(\alpha r_{ij})}{r_{ij}} +
+   \frac{1}{2} \sum_{j \neq i}
+   \frac{q_i q_j {\rm erf}(\alpha r_{ij})}{r_{ij}} \qquad r < r_c
+
+where :math:`\alpha` is the damping parameter, and *erf()* and *erfc()*
+are error-function and complementary error-function terms.  This
+potential is essentially a short-range, spherically-truncated,
+charge-neutralized, force-shifted, pairwise *1/r* summation.  With a
+manipulation of adding and subtracting a self term (for i = j) to the
+first and second term on the right-hand-side, respectively, and a small
+enough :math:`\alpha` damping parameter, the second term shrinks and the
+potential becomes a rapidly-converging real-space summation.  With a
+long enough cutoff and small enough :math:`\alpha` parameter, the
+electric field calculated by the Wolf summation method approaches that
+computed using the Ewald sum.
+
+The value of the electric field components will be 0.0 for atoms not in
+the specified compute group.
+
+When the *limit* keyword is used, only contributions from atoms in the
+selected group will be considered, otherwise contributions from all
+atoms within the cutoff are included.
+
+When the *cutoff* keyword is used, the cutoff used for the electric
+field approximation can be set explicitly.  By default it is the largest
+cutoff of any pair style force computation.
+
+.. admonition:: Computational Efficiency
+   :class: note
+
+   This compute will loop over a full neighbor list just like a pair
+   style does when computing forces, thus it can be quite time consuming
+   and slow down a calculation significantly when its data is used in
+   every time step.  The :doc:`compute efield/atom
+   <compute_efield_atom>` command of the DIELECTRIC package is more
+   efficient in comparison, since the electric field data is collected
+   and stored as part of the force computation at next to no extra
+   computational cost.
+
+Output info
+"""""""""""
+
+This compute calculates a per-atom vector, which can be accessed by
+any command that uses per-atom values from a compute as input.  See
+the :doc:`Howto output <Howto_output>` page for an overview of
+LAMMPS output options.
+
+The vector contains 3 values per atom which are the x-, y-, and
+z-direction electric field components in force units.
+
+Restrictions
+""""""""""""
+
+This compute is part of the EXTRA-COMPUTE package. It is only enabled if
+LAMMPS was built with that package.
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_style coul/wolf <pair_coul>`,
+:doc:`compute efield/atom <compute_efield_atom>`
+
+Default
+"""""""
+
+The option defaults are *limit* = all and *cutoff* = largest cutoff
+for pair styles.
+
+----------
+
+.. _Wolf4:
+
+**(Wolf)** D. Wolf, P. Keblinski, S. R. Phillpot, J. Eggebrecht, J Chem
+Phys, 110, 8254 (1999).

src/.gitignore

@@ -534,6 +534,8 @@
 /compute_dpd.h
 /compute_dpd_atom.cpp
 /compute_dpd_atom.h
+/compute_efield_wolf_atom.cpp
+/compute_efield_wolf_atom.h
 /compute_entropy_atom.cpp
 /compute_entropy_atom.h
 /compute_erotate_asphere.cpp

src/EXTRA-COMPUTE/compute_efield_wolf_atom.cpp

@@ -0,0 +1,200 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_efield_wolf_atom.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "group.h"
+#include "math_const.h"
+#include "memory.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "neighbor.h"
+#include "pair.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using MathConst::MY_PIS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeEfieldWolfAtom::ComputeEfieldWolfAtom(LAMMPS *lmp, int narg, char **arg) :
+    Compute(lmp, narg, arg), group2(nullptr), efield(nullptr)
+
+{
+  if (narg < 4) utils::missing_cmd_args(FLERR, "compute efield/atom/wolf", error);
+
+  peratom_flag = 1;
+  size_peratom_cols = 3;
+
+  nmax = -1;
+  group2 = utils::strdup("all");
+  jgroupbit = group->bitmask[0];
+  cutoff_flag = 0;
+  cutoff = 0.0;
+
+  alpha = utils::numeric(FLERR, arg[3], false, lmp);
+
+  int iarg = 4;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg], "group") == 0) {
+      if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "compute efield/atom/wolf group", error);
+      delete[] group2;
+      group2 = utils::strdup(arg[iarg + 1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg], "cutoff") == 0) {
+      if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "compute efield/atom/wolf cutoff", error);
+      cutoff = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
+      cutoff_flag = 1;
+      iarg += 2;
+    } else
+      error->all(FLERR, "Unknown compute {} keyword: {}", style, arg[iarg]);
+  }
+
+  // sanity checks
+
+  if (alpha <= 0.0) error->all(FLERR, "Compute efield/atom/wolf alpha value {} is invalid", alpha);
+  if (cutoff_flag && cutoff <= 0.0)
+    error->all(FLERR, "Compute efield/atom/wolf cutoff {} is invalid", cutoff);
+
+  jgroup = group->find(group2);
+  if (jgroup < 0) error->all(FLERR, "Compute efield/atom/wolf group {} does not exist", group2);
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeEfieldWolfAtom::~ComputeEfieldWolfAtom()
+{
+  memory->destroy(efield);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeEfieldWolfAtom::init()
+{
+  if (!atom->q_flag)
+    error->all(FLERR,"Compute efield/wolf/atom requires atom attribute q");
+  if (atom->mu_flag && (comm->me == 0))
+    error->warning(FLERR,"Compute efield/wolf/atom does not support per-atom dipoles");
+
+  // need an occasional full neighbor list
+  auto req = neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
+  if (cutoff_flag) req->set_cutoff(cutoff);
+
+  jgroup = group->find(group2);
+  if (jgroup < 0) error->all(FLERR, "Compute efield/atom/wolf group {} does not exist", group2);
+  jgroupbit = group->bitmask[jgroup];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeEfieldWolfAtom::init_list(int /*id*/, NeighList *ptr)
+{
+  list = ptr;
+}
+
+// clang-format off
+/* ---------------------------------------------------------------------- */
+
+void ComputeEfieldWolfAtom::compute_peratom()
+{
+  invoked_peratom = update->ntimestep;
+
+  // grow result array if necessary and clear
+
+  if (atom->nmax > nmax) {
+    memory->destroy(efield);
+    nmax = atom->nmax;
+    memory->create(efield,nmax,3,"efield/atom/wolf:efield");
+    array_atom = efield;
+  }
+  memset(&efield[0][0], 0, sizeof(double)*nmax*3);
+
+  // invoke neighbor list build (will copy or build if necessary)
+
+  neighbor->build_one(list);
+
+  const auto inum = list->inum;
+  const auto ilist = list->ilist;
+  const auto numneigh = list->numneigh;
+  const auto firstneigh = list->firstneigh;
+
+  // compute coulomb force according to Wolf sum approximation
+  const double * const * const x = atom->x;
+  const int * const mask = atom->mask;
+  const double * const q = atom->q;
+  const double * const special_coul = force->special_coul;
+
+  const double qqrd2e = force->qqrd2e;
+
+  if (!cutoff_flag && force->pair) cutoff = force->pair->cutforce;
+  const double cutsq = cutoff*cutoff;
+  const double e_shift = erfc(alpha * cutoff) / cutoff;
+  const double f_shift = -(e_shift + 2.0 * alpha / MY_PIS * exp(-alpha * alpha * cutsq)) / cutoff;
+
+#if defined(_OPENMP)
+#pragma omp parallel for
+#endif
+  for (int ii = 0; ii < inum; ii++) {
+    const int i = ilist[ii];
+    if (mask[i] & groupbit) {
+      const double xtmp = x[i][0];
+      const double ytmp = x[i][1];
+      const double ztmp = x[i][2];
+      const auto jlist = firstneigh[i];
+      const auto jnum = numneigh[i];
+
+      for (int jj = 0; jj < jnum; jj++) {
+        int j = jlist[jj];
+        const double factor_coul = special_coul[sbmask(j)];
+        j &= NEIGHMASK;
+        if (mask[j] & jgroupbit) {
+          const double delx = xtmp - x[j][0];
+          const double dely = ytmp - x[j][1];
+          const double delz = ztmp - x[j][2];
+          const double rsq = delx*delx + dely*dely + delz*delz;
+          if ((rsq > 0.0) && (rsq < cutsq)) {
+            const double r = sqrt(rsq);
+            double prefactor = qqrd2e * q[j] / r;
+            double erfcc = erfc(alpha * r);
+            double erfcd = exp(-alpha * alpha * r * r);
+            double dvdrr = (erfcc / rsq + 2.0 * alpha / MY_PIS * erfcd / r) + f_shift;
+            double forcecoul = dvdrr * rsq * prefactor;
+            if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
+            forcecoul /= rsq;
+
+            efield[i][0] += delx * forcecoul;
+            efield[i][1] += dely * forcecoul;
+            efield[i][2] += delz * forcecoul;
+          }
+        }
+      }
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeEfieldWolfAtom::memory_usage()
+{
+  double bytes = 3.0 * nmax * sizeof(double);
+  return bytes;
+}

src/EXTRA-COMPUTE/compute_efield_wolf_atom.h

@@ -0,0 +1,53 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+// clang-format off
+ComputeStyle(efield/wolf/atom,ComputeEfieldWolfAtom);
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_EFIELD_WOLF_ATOM_H
+#define LMP_COMPUTE_EFIELD_WOLF_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeEfieldWolfAtom : public Compute {
+ public:
+  ComputeEfieldWolfAtom(class LAMMPS *, int, char **);
+  ~ComputeEfieldWolfAtom() override;
+
+  void init() override;
+  void init_list(int, class NeighList *) override;
+
+  void compute_peratom() override;
+  double memory_usage() override;
+
+ private:
+  int nmax, cutoff_flag;
+  double cutoff;
+  double alpha;
+  class NeighList *list;
+
+  char *group2;
+  int jgroup, jgroupbit;
+
+  double **efield;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif