diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 682a75f201..5c5f059e96 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -57,6 +57,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`dpd/atom ` * :doc:`edpd/temp/atom ` * :doc:`efield/atom ` + * :doc:`efield/wolf/atom ` * :doc:`entropy/atom ` * :doc:`erotate/asphere ` * :doc:`erotate/rigid ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 7d1e423e93..5f48401cbe 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -202,6 +202,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`dpd/atom ` - per-particle values of internal conductive energy, internal mechanical energy, chemical energy, and internal temperature * :doc:`edpd/temp/atom ` - per-atom temperature for each eDPD particle in a group * :doc:`efield/atom ` - electric field at each atom +* :doc:`efield/wolf/atom ` - electric field at each atom * :doc:`entropy/atom ` - pair entropy fingerprint of each atom * :doc:`erotate/asphere ` - rotational energy of aspherical particles * :doc:`erotate/rigid ` - rotational energy of rigid bodies diff --git a/doc/src/compute_efield_wolf_atom.rst b/doc/src/compute_efield_wolf_atom.rst new file mode 100644 index 0000000000..c2679354d6 --- /dev/null +++ b/doc/src/compute_efield_wolf_atom.rst @@ -0,0 +1,127 @@ +.. 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 ` 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 `, 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 + ` 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 ` page for an overview of +LAMMPS output options. + +The vector contains 4 values per atom. The first 3 per-atom vector +values will be the x-, y-, and z-direction electric field components +in force units the 4th per-atom vector value is the field magnitude. + +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 `, +:doc:`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). diff --git a/src/.gitignore b/src/.gitignore index 0fbd54afcd..414dfab4d1 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -521,6 +521,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 diff --git a/src/EXTRA-COMPUTE/compute_efield_wolf_atom.cpp b/src/EXTRA-COMPUTE/compute_efield_wolf_atom.cpp new file mode 100644 index 0000000000..a5e287f2db --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_efield_wolf_atom.cpp @@ -0,0 +1,197 @@ +/* ---------------------------------------------------------------------- + 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 +#include + +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 = 4; + + 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], "limit") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "compute efield/atom/wolf limit", 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() +{ + // 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,4,"efield/atom/wolf:efield"); + array_atom = efield; + } + memset(&efield[0][0], 0, sizeof(double)*nmax*4); + + // 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 < 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; + } + } + } + efield[i][3] = sqrt(efield[i][0]*efield[i][0] + efield[i][1]*efield[i][1] + + efield[i][2]*efield[i][2]); + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeEfieldWolfAtom::memory_usage() +{ + double bytes = 4.0 * nmax * sizeof(double); + return bytes; +} diff --git a/src/EXTRA-COMPUTE/compute_efield_wolf_atom.h b/src/EXTRA-COMPUTE/compute_efield_wolf_atom.h new file mode 100644 index 0000000000..445f79cbd0 --- /dev/null +++ b/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