add compute dipole/tip4p and compute dipole/tip4p/chunk

doc/src/Commands_compute.rst

@@ -52,6 +52,8 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`dilatation/atom <compute_dilatation_atom>`
    * :doc:`dipole <compute_dipole>`
    * :doc:`dipole/chunk <compute_dipole_chunk>`
+   * :doc:`dipole/tip4p <compute_dipole>`
+   * :doc:`dipole/tip4p/chunk <compute_dipole_chunk>`
    * :doc:`displace/atom <compute_displace_atom>`
    * :doc:`dpd <compute_dpd>`
    * :doc:`dpd/atom <compute_dpd_atom>`

doc/src/compute.rst

@@ -206,6 +206,8 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`dilatation/atom <compute_dilatation_atom>` - Peridynamic dilatation for each atom
 * :doc:`dipole <compute_dipole>` - dipole vector and total dipole
 * :doc:`dipole/chunk <compute_dipole_chunk>` - dipole vector and total dipole for each chunk
+* :doc:`dipole/tip4p <compute_dipole>` - dipole vector and total dipole with TIP4P pair style
+* :doc:`dipole/tip4p/chunk <compute_dipole_chunk>` - dipole vector and total dipole for each chunk with TIP4P pair style
 * :doc:`displace/atom <compute_displace_atom>` - displacement of each atom
 * :doc:`dpd <compute_dpd>` - total values of internal conductive energy, internal mechanical energy, chemical energy, and harmonic average of internal temperature
 * :doc:`dpd/atom <compute_dpd_atom>` - per-particle values of internal conductive energy, internal mechanical energy, chemical energy, and internal temperature

doc/src/compute_dipole.rst

@@ -1,6 +1,10 @@
 .. index:: compute dipole
+.. index:: compute dipole/tip4p
 
 compute dipole command
+======================
+
+compute dipole/tip4p command
 ============================
 
 Syntax
@@ -8,10 +12,10 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   compute ID group-ID dipole arg
+   compute ID group-ID style arg
 
 * ID, group-ID are documented in :doc:`compute <compute>` command
-* dipole = style name of this compute command
+* style = *dipole* or *dipole/tip4p*
 * arg = *mass* or *geometry* = use COM or geometric center for charged chunk correction (optional)
 
 Examples
@@ -21,6 +25,7 @@ Examples
 
    compute 1 fluid dipole
    compute dw water dipole geometry
+   compute dw water dipole/tip4p
 
 Description
 """""""""""
@@ -28,13 +33,20 @@ Description
 Define a computation that calculates the dipole vector and total dipole
 for a group of atoms.
 
-This compute calculates the x,y,z coordinates of the dipole vector
+These computes calculate the x,y,z coordinates of the dipole vector and
-and the total dipole moment for the atoms in the compute group.
+the total dipole moment for the atoms in the compute group.  This
-This includes all effects due to atoms passing through periodic boundaries.
+includes all effects due to atoms passing through periodic boundaries.
-For a group with a net charge the resulting dipole is made position independent
+For a group with a net charge the resulting dipole is made position
-by subtracting the position vector of the center of mass or geometric center
+independent by subtracting the position vector of the center of mass or
-times the net charge from the computed dipole vector. Both per-atom charges
+geometric center times the net charge from the computed dipole
-and per-atom dipole moments, if present, contribute to the computed dipole.
+vector.  Both per-atom charges and per-atom dipole moments, if present,
+contribute to the computed dipole.
+
+.. versionadded:: TBD
+
+Compute *dipole/tip4p* includes adjustments for the charge carrying
+point M in molecules with TIP4P water geometry.  The corresponding
+parameters are extracted from the pair style.
 
 .. note::
 
@@ -49,10 +61,10 @@ and per-atom dipole moments, if present, contribute to the computed dipole.
 Output info
 """""""""""
 
-This compute calculations a global scalar containing the magnitude of
+These computes calculate a global scalar containing the magnitude of the
-the computed dipole moment and a global vector of length 3 with the
+computed dipole moment and a global vector of length 3 with the dipole
-dipole vector.  See the :doc:`Howto output <Howto_output>` page for
+vector.  See the :doc:`Howto output <Howto_output>` page for an overview
-an overview of LAMMPS output options.
+of LAMMPS output options.
 
 The computed values are "intensive".  The array values will be in
 dipole units (i.e., charge :doc:`units <units>` times distance
@@ -60,7 +72,12 @@ dipole units (i.e., charge :doc:`units <units>` times distance
 
 Restrictions
 """"""""""""
- none
+
+Compute style *dipole/tip4p* is part of the EXTRA-COMPUTE package. It is
+only enabled if LAMMPS was built with that package.  See the :doc:`Build
+package <Build_package>` page for more info.
+
+Compute style *dipole/tip4p* can only be used with tip4p pair styles.
 
 Related commands
 """"""""""""""""

doc/src/compute_dipole_chunk.rst

@@ -1,17 +1,21 @@
 .. index:: compute dipole/chunk
+.. index:: compute dipole/tip4p/chunk
 
 compute dipole/chunk command
 ============================
 
+compute dipole/tip4p/chunk command
+==================================
+
 Syntax
 """"""
 
 .. code-block:: LAMMPS
 
-   compute ID group-ID dipole/chunk chunkID arg
+   compute ID group-ID style chunkID arg
 
 * ID, group-ID are documented in :doc:`compute <compute>` command
-* dipole/chunk = style name of this compute command
+* style = *dipole/chunk* or *dipole/tip4p/chunk*
 * chunkID = ID of :doc:`compute chunk/atom <compute_chunk_atom>` command
 * arg = *mass* or *geometry* = use COM or geometric center for charged chunk correction (optional)
 
@@ -38,13 +42,20 @@ or atoms in a spatial bin.  See the :doc:`compute chunk/atom
 details of how chunks can be defined and examples of how they can be
 used to measure properties of a system.
 
-This compute calculates the :math:`(x,y,z)` coordinates of the dipole vector
+These computes calculate the :math:`(x,y,z)` coordinates of the dipole
-and the total dipole moment for each chunk, which includes all effects due
+vector and the total dipole moment for each chunk, which includes all
-to atoms passing through periodic boundaries. For chunks with a net
+effects due to atoms passing through periodic boundaries.  For chunks
-charge the resulting dipole is made position independent by subtracting
+with a net charge the resulting dipole is made position independent by
-the position vector of the center of mass or geometric center times the
+subtracting the position vector of the center of mass or geometric
-net charge from the computed dipole vector. Both per-atom charges and
+center times the net charge from the computed dipole vector.  Both
-per-atom dipole moments, if present, contribute to the computed dipole.
+per-atom charges and per-atom dipole moments, if present, contribute to
+the computed dipole.
+
+.. versionadded:: TBD
+
+Compute *dipole/tip4p/chunk* includes adjustments for the charge
+carrying point M in molecules with TIP4P water geometry.  The
+corresponding parameters are extracted from the pair style.
 
 Note that only atoms in the specified group contribute to the
 calculation.  The :doc:`compute chunk/atom <compute_chunk_atom>` command
@@ -78,12 +89,12 @@ command, for example:
 Output info
 """""""""""
 
-This compute calculates a global array where the number of rows = the
+These computes calculate a global array where the number of rows = the
 number of chunks *Nchunk* as calculated by the specified :doc:`compute
-chunk/atom <compute_chunk_atom>` command.  The number of columns is 4 for
+chunk/atom <compute_chunk_atom>` command.  The number of columns is 4
-the :math:`(x,y,z)` dipole vector components and the total dipole of each
+for the :math:`(x,y,z)` dipole vector components and the total dipole of
-chunk. These values can be accessed by any command that uses global
+each chunk.  These values can be accessed by any command that uses
-array values from a compute as input.  See the :doc:`Howto output
+global array values from a compute as input.  See the :doc:`Howto output
 <Howto_output>` page for an overview of LAMMPS output options.
 
 The array values are "intensive".  The array values will be in
@@ -92,7 +103,13 @@ dipole units (i.e., charge :doc:`units <units>` times distance
 
 Restrictions
 """"""""""""
- none
+
+Compute style *dipole/tip4p/chunk* is part of the EXTRA-COMPUTE
+package. It is only enabled if LAMMPS was built with that package.  See
+the :doc:`Build package <Build_package>` page for more info.
+
+Compute style *dipole/tip4p/chunk* can only be used with tip4p pair
+styles.
 
 Related commands
 """"""""""""""""

src/.gitignore

@@ -534,6 +534,10 @@
 /compute_damage_atom.h
 /compute_dilatation_atom.cpp
 /compute_dilatation_atom.h
+/compute_dipole_tip4p.cpp
+/compute_dipole_tip4p.h
+/compute_dipole_tip4p_chunk.cpp
+/compute_dipole_tip4p_chunk.h
 /compute_dpd.cpp
 /compute_dpd.h
 /compute_dpd_atom.cpp

src/EXTRA-COMPUTE/compute_dipole_tip4p.cpp

@@ -0,0 +1,219 @@
+/* ----------------------------------------------------------------------
+   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_dipole_tip4p.h"
+
+#include "angle.h"
+#include "atom.h"
+#include "bond.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "pair.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+
+enum { MASSCENTER, GEOMCENTER };
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDipoleTIP4P::ComputeDipoleTIP4P(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg)
+{
+  if ((narg < 3) || (narg > 4)) error->all(FLERR, "Illegal compute dipole command");
+
+  scalar_flag = 1;
+  vector_flag = 1;
+  size_vector = 3;
+  extscalar = 0;
+  extvector = 0;
+
+  vector = new double[size_vector];
+  vector[0] = vector[1] = vector[2] = 0.0;
+  usecenter = MASSCENTER;
+
+  if (narg == 4) {
+    if (utils::strmatch(arg[3], "^geom"))
+      usecenter = GEOMCENTER;
+    else if (strcmp(arg[3], "mass") == 0)
+      usecenter = MASSCENTER;
+    else
+      error->all(FLERR, "Illegal compute dipole command");
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDipoleTIP4P::~ComputeDipoleTIP4P()
+{
+  delete[] vector;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4P::init()
+{
+  int itmp;
+  double *p_qdist = (double *) force->pair->extract("qdist", itmp);
+  int *p_typeO = (int *) force->pair->extract("typeO", itmp);
+  int *p_typeH = (int *) force->pair->extract("typeH", itmp);
+  int *p_typeA = (int *) force->pair->extract("typeA", itmp);
+  int *p_typeB = (int *) force->pair->extract("typeB", itmp);
+  if (!p_qdist || !p_typeO || !p_typeH || !p_typeA || !p_typeB)
+    error->all(FLERR, "Pair style is incompatible with compute {}", style);
+  typeO = *p_typeO;
+  typeH = *p_typeH;
+  int typeA = *p_typeA;
+  int typeB = *p_typeB;
+
+  if (force->angle == nullptr || force->bond == nullptr || force->angle->setflag == nullptr ||
+      force->bond->setflag == nullptr)
+    error->all(FLERR, "Bond and angle potentials must be defined for compute {}", style);
+  if ((typeA < 1) || (typeA > atom->nangletypes) || (force->angle->setflag[typeA] == 0))
+    error->all(FLERR, "Bad TIP4P angle type for compute {}", style);
+  if ((typeB < 1) || (typeB > atom->nbondtypes) || (force->bond->setflag[typeB] == 0))
+    error->all(FLERR, "Bad TIP4P bond type for PPPM/TIP4P");
+  double theta = force->angle->equilibrium_angle(typeA);
+  double blen = force->bond->equilibrium_distance(typeB);
+  alpha = *p_qdist / (cos(0.5 * theta) * blen);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4P::compute_vector()
+{
+  invoked_vector = update->ntimestep;
+
+  const auto x = atom->x;
+  const auto mask = atom->mask;
+  const auto type = atom->type;
+  const auto image = atom->image;
+  const auto mass = atom->mass;
+  const auto rmass = atom->rmass;
+  const auto q = atom->q;
+  const auto mu = atom->mu;
+  const auto nlocal = atom->nlocal;
+
+  double dipole[3] = {0.0, 0.0, 0.0};
+  double comproc[3] = {0.0, 0.0, 0.0};
+  double com[3] = {0.0, 0.0, 0.0};
+  double masstotal = 0.0;
+  double chrgtotal = 0.0;
+  double massproc = 0.0;
+  double chrgproc = 0.0;
+  double unwrap[3], xM[3];
+  double *xi;
+
+  for (int i = 0; i < nlocal; ++i) {
+    if (mask[i] & groupbit) {
+      double massone = 1.0;    // for usecenter == GEOMCENTER
+      if (usecenter == MASSCENTER) {
+        if (rmass)
+          massone = rmass[i];
+        else
+          massone = mass[type[i]];
+      }
+      massproc += massone;
+      if (atom->q_flag) chrgproc += q[i];
+      domain->unmap(x[i], image[i], unwrap);
+      comproc[0] += unwrap[0] * massone;
+      comproc[1] += unwrap[1] * massone;
+      comproc[2] += unwrap[2] * massone;
+    }
+  }
+  MPI_Allreduce(&massproc, &masstotal, 1, MPI_DOUBLE, MPI_SUM, world);
+  MPI_Allreduce(&chrgproc, &chrgtotal, 1, MPI_DOUBLE, MPI_SUM, world);
+  MPI_Allreduce(comproc, com, 3, MPI_DOUBLE, MPI_SUM, world);
+
+  if (masstotal > 0.0) {
+    com[0] /= masstotal;
+    com[1] /= masstotal;
+    com[2] /= masstotal;
+  }
+
+  // compute dipole moment
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      if (type[i] == typeO) {
+        find_M(i,xM);
+        xi = xM;
+      } else {
+        xi = x[i];
+      }
+      domain->unmap(xi, image[i], unwrap);
+      if (atom->q_flag) {
+        dipole[0] += q[i] * unwrap[0];
+        dipole[1] += q[i] * unwrap[1];
+        dipole[2] += q[i] * unwrap[2];
+      }
+      if (atom->mu_flag) {
+        dipole[0] += mu[i][0];
+        dipole[1] += mu[i][1];
+        dipole[2] += mu[i][2];
+      }
+    }
+  }
+
+  MPI_Allreduce(dipole, vector, 3, MPI_DOUBLE, MPI_SUM, world);
+
+  // correct for position dependence with a net charged group
+  vector[0] -= chrgtotal * com[0];
+  vector[1] -= chrgtotal * com[1];
+  vector[2] -= chrgtotal * com[2];
+}
+
+/* ---------------------------------------------------------------------- */
+
+double ComputeDipoleTIP4P::compute_scalar()
+{
+  if (invoked_vector != update->ntimestep) compute_vector();
+
+  invoked_scalar = update->ntimestep;
+  scalar = sqrt(vector[0] * vector[0] + vector[1] * vector[1] + vector[2] * vector[2]);
+  return scalar;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4P::find_M(int i, double *xM)
+{
+  double **x = atom->x;
+
+  int iH1 = atom->map(atom->tag[i] + 1);
+  int iH2 = atom->map(atom->tag[i] + 2);
+
+  if ((iH1 == -1) || (iH2 == -1)) error->one(FLERR,"TIP4P hydrogen is missing");
+  if ((atom->type[iH1] != typeH) || (atom->type[iH2] != typeH))
+    error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+
+  // set iH1,iH2 to index of closest image to O
+
+  iH1 = domain->closest_image(i,iH1);
+  iH2 = domain->closest_image(i,iH2);
+
+  double delx1 = x[iH1][0] - x[i][0];
+  double dely1 = x[iH1][1] - x[i][1];
+  double delz1 = x[iH1][2] - x[i][2];
+
+  double delx2 = x[iH2][0] - x[i][0];
+  double dely2 = x[iH2][1] - x[i][1];
+  double delz2 = x[iH2][2] - x[i][2];
+
+  xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2);
+  xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2);
+  xM[2] = x[i][2] + alpha * 0.5 * (delz1 + delz2);
+}

src/EXTRA-COMPUTE/compute_dipole_tip4p.h

@@ -0,0 +1,46 @@
+/* -*- 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(dipole/tip4p,ComputeDipoleTIP4P);
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_DIPOLE_TIP4P_H
+#define LMP_COMPUTE_DIPOLE_TIP4P_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeDipoleTIP4P : public Compute {
+ public:
+  ComputeDipoleTIP4P(class LAMMPS *, int, char **);
+  ~ComputeDipoleTIP4P() override;
+  void init() override;
+  void compute_vector() override;
+  double compute_scalar() override;
+
+ private:
+  int usecenter;
+
+  int typeO, typeH;
+  double alpha;
+  void find_M(int i, double *xM);
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp

@@ -0,0 +1,364 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   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_dipole_tip4p_chunk.h"
+
+#include "angle.h"
+#include "atom.h"
+#include "bond.h"
+#include "comm.h"
+#include "compute_chunk_atom.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "pair.h"
+#include "math_special.h"
+#include "memory.h"
+#include "modify.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace MathSpecial;
+
+enum { MASSCENTER, GEOMCENTER };
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDipoleTIP4PChunk::ComputeDipoleTIP4PChunk(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg),
+  idchunk(nullptr), massproc(nullptr), masstotal(nullptr), chrgproc(nullptr),
+  chrgtotal(nullptr), com(nullptr),
+  comall(nullptr), dipole(nullptr), dipoleall(nullptr)
+{
+  if ((narg != 4) && (narg != 5))
+    error->all(FLERR,"Illegal compute dipole/chunk command");
+
+  array_flag = 1;
+  size_array_cols = 4;
+  size_array_rows = 0;
+  size_array_rows_variable = 1;
+  extarray = 0;
+
+  // ID of compute chunk/atom
+
+  idchunk = utils::strdup(arg[3]);
+
+  usecenter = MASSCENTER;
+
+  if (narg == 5) {
+    if (strncmp(arg[4],"geom",4) == 0) usecenter = GEOMCENTER;
+    else if (strcmp(arg[4],"mass") == 0) usecenter = MASSCENTER;
+    else error->all(FLERR,"Illegal compute dipole/chunk command");
+  }
+
+  ComputeDipoleTIP4PChunk::init();
+
+  // chunk-based data
+
+  nchunk = 1;
+  maxchunk = 0;
+  allocate();
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeDipoleTIP4PChunk::~ComputeDipoleTIP4PChunk()
+{
+  delete[] idchunk;
+  memory->destroy(massproc);
+  memory->destroy(masstotal);
+  memory->destroy(chrgproc);
+  memory->destroy(chrgtotal);
+  memory->destroy(com);
+  memory->destroy(comall);
+  memory->destroy(dipole);
+  memory->destroy(dipoleall);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4PChunk::init()
+{
+  int icompute = modify->find_compute(idchunk);
+  if (icompute < 0)
+    error->all(FLERR,"Chunk/atom compute does not exist for compute dipole/chunk");
+  cchunk = dynamic_cast<ComputeChunkAtom *>(modify->compute[icompute]);
+  if (strcmp(cchunk->style,"chunk/atom") != 0)
+    error->all(FLERR,"Compute dipole/chunk does not use chunk/atom compute");
+
+  int itmp;
+  double *p_qdist = (double *) force->pair->extract("qdist", itmp);
+  int *p_typeO = (int *) force->pair->extract("typeO", itmp);
+  int *p_typeH = (int *) force->pair->extract("typeH", itmp);
+  int *p_typeA = (int *) force->pair->extract("typeA", itmp);
+  int *p_typeB = (int *) force->pair->extract("typeB", itmp);
+  if (!p_qdist || !p_typeO || !p_typeH || !p_typeA || !p_typeB)
+    error->all(FLERR, "Pair style is incompatible with compute {}", style);
+  typeO = *p_typeO;
+  typeH = *p_typeH;
+  int typeA = *p_typeA;
+  int typeB = *p_typeB;
+
+  if (force->angle == nullptr || force->bond == nullptr || force->angle->setflag == nullptr ||
+      force->bond->setflag == nullptr)
+    error->all(FLERR, "Bond and angle potentials must be defined for compute {}", style);
+  if ((typeA < 1) || (typeA > atom->nangletypes) || (force->angle->setflag[typeA] == 0))
+    error->all(FLERR, "Bad TIP4P angle type for compute {}", style);
+  if ((typeB < 1) || (typeB > atom->nbondtypes) || (force->bond->setflag[typeB] == 0))
+    error->all(FLERR, "Bad TIP4P bond type for PPPM/TIP4P");
+  double theta = force->angle->equilibrium_angle(typeA);
+  double blen = force->bond->equilibrium_distance(typeB);
+  alpha = *p_qdist / (cos(0.5 * theta) * blen);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4PChunk::compute_array()
+{
+  int i,index;
+  double massone;
+  double unwrap[3];
+
+  invoked_array = update->ntimestep;
+
+  // compute chunk/atom assigns atoms to chunk IDs
+  // extract ichunk index vector from compute
+  // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
+
+  nchunk = cchunk->setup_chunks();
+  cchunk->compute_ichunk();
+  int *ichunk = cchunk->ichunk;
+
+  if (nchunk > maxchunk) allocate();
+  size_array_rows = nchunk;
+
+  // zero local per-chunk values
+
+  for (i = 0; i < nchunk; i++) {
+    massproc[i] = chrgproc[i] = 0.0;
+    com[i][0] = com[i][1] = com[i][2] = 0.0;
+    dipole[i][0] = dipole[i][1] = dipole[i][2] = dipole[i][3] = 0.0;
+  }
+
+  // compute COM for each chunk
+
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  imageint *image = atom->image;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  double *q = atom->q;
+  double **mu = atom->mu;
+  double xM[3];
+  double *xi;
+
+  int nlocal = atom->nlocal;
+
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      index = ichunk[i]-1;
+      if (index < 0) continue;
+      if (usecenter == MASSCENTER) {
+        if (rmass) massone = rmass[i];
+        else massone = mass[type[i]];
+      } else massone = 1.0;     // usecenter == GEOMCENTER
+      domain->unmap(x[i],image[i],unwrap);
+      massproc[index] += massone;
+      if (atom->q_flag) chrgproc[index] += q[i];
+      com[index][0] += unwrap[0] * massone;
+      com[index][1] += unwrap[1] * massone;
+      com[index][2] += unwrap[2] * massone;
+    }
+
+  MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(chrgproc,chrgtotal,nchunk,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world);
+
+  for (i = 0; i < nchunk; i++) {
+    if (masstotal[i] > 0.0) {
+      comall[i][0] /= masstotal[i];
+      comall[i][1] /= masstotal[i];
+      comall[i][2] /= masstotal[i];
+    }
+  }
+
+  // compute dipole for each chunk
+
+  for (i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      index = ichunk[i]-1;
+      if (index < 0) continue;
+
+      if (type[i] == typeO) {
+        find_M(i,xM);
+        xi = xM;
+      } else {
+        xi = x[i];
+      }
+      domain->unmap(xi,image[i],unwrap);
+      if (atom->q_flag) {
+        dipole[index][0] += q[i]*unwrap[0];
+        dipole[index][1] += q[i]*unwrap[1];
+        dipole[index][2] += q[i]*unwrap[2];
+      }
+      if (atom->mu_flag) {
+        dipole[index][0] += mu[i][0];
+        dipole[index][1] += mu[i][1];
+        dipole[index][2] += mu[i][2];
+      }
+    }
+  }
+
+  MPI_Allreduce(&dipole[0][0],&dipoleall[0][0],4*nchunk,MPI_DOUBLE,MPI_SUM,world);
+
+  for (i = 0; i < nchunk; i++) {
+    // correct for position dependence with charged chunks
+    dipoleall[i][0] -= chrgtotal[i]*comall[i][0];
+    dipoleall[i][1] -= chrgtotal[i]*comall[i][1];
+    dipoleall[i][2] -= chrgtotal[i]*comall[i][2];
+    // compute total dipole moment
+    dipoleall[i][3] = sqrt(square(dipoleall[i][0])
+                           + square(dipoleall[i][1])
+                           + square(dipoleall[i][2]));
+  }
+}
+
+/* ----------------------------------------------------------------------
+   lock methods: called by fix ave/time
+   these methods ensure vector/array size is locked for Nfreq epoch
+     by passing lock info along to compute chunk/atom
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   increment lock counter
+------------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4PChunk::lock_enable()
+{
+  cchunk->lockcount++;
+}
+
+/* ----------------------------------------------------------------------
+   decrement lock counter in compute chunk/atom, it if still exists
+------------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4PChunk::lock_disable()
+{
+  int icompute = modify->find_compute(idchunk);
+  if (icompute >= 0) {
+    cchunk = dynamic_cast<ComputeChunkAtom *>(modify->compute[icompute]);
+    cchunk->lockcount--;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   calculate and return # of chunks = length of vector/array
+------------------------------------------------------------------------- */
+
+int ComputeDipoleTIP4PChunk::lock_length()
+{
+  nchunk = cchunk->setup_chunks();
+  return nchunk;
+}
+
+/* ----------------------------------------------------------------------
+   set the lock from startstep to stopstep
+------------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4PChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
+{
+  cchunk->lock(fixptr,startstep,stopstep);
+}
+
+/* ----------------------------------------------------------------------
+   unset the lock
+------------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4PChunk::unlock(Fix *fixptr)
+{
+  cchunk->unlock(fixptr);
+}
+
+/* ----------------------------------------------------------------------
+   free and reallocate per-chunk arrays
+------------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4PChunk::allocate()
+{
+  memory->destroy(massproc);
+  memory->destroy(masstotal);
+  memory->destroy(chrgproc);
+  memory->destroy(chrgtotal);
+  memory->destroy(com);
+  memory->destroy(comall);
+  memory->destroy(dipole);
+  memory->destroy(dipoleall);
+  maxchunk = nchunk;
+  memory->create(massproc,maxchunk,"dipole/chunk:massproc");
+  memory->create(masstotal,maxchunk,"dipole/chunk:masstotal");
+  memory->create(chrgproc,maxchunk,"dipole/chunk:chrgproc");
+  memory->create(chrgtotal,maxchunk,"dipole/chunk:chrgtotal");
+  memory->create(com,maxchunk,3,"dipole/chunk:com");
+  memory->create(comall,maxchunk,3,"dipole/chunk:comall");
+  memory->create(dipole,maxchunk,4,"dipole/chunk:dipole");
+  memory->create(dipoleall,maxchunk,4,"dipole/chunk:dipoleall");
+  array = dipoleall;
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeDipoleTIP4PChunk::memory_usage()
+{
+  double bytes = (bigint) maxchunk * 2 * sizeof(double);
+  bytes += (double)maxchunk * 2*3 * sizeof(double);
+  bytes += (double)maxchunk * 2*4 * sizeof(double);
+  return bytes;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeDipoleTIP4PChunk::find_M(int i, double *xM)
+{
+  double **x = atom->x;
+
+  int iH1 = atom->map(atom->tag[i] + 1);
+  int iH2 = atom->map(atom->tag[i] + 2);
+
+  if ((iH1 == -1) || (iH2 == -1)) error->one(FLERR,"TIP4P hydrogen is missing");
+  if ((atom->type[iH1] != typeH) || (atom->type[iH2] != typeH))
+    error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+
+  // set iH1,iH2 to index of closest image to O
+
+  iH1 = domain->closest_image(i,iH1);
+  iH2 = domain->closest_image(i,iH2);
+
+  double delx1 = x[iH1][0] - x[i][0];
+  double dely1 = x[iH1][1] - x[i][1];
+  double delz1 = x[iH1][2] - x[i][2];
+
+  double delx2 = x[iH2][0] - x[i][0];
+  double dely2 = x[iH2][1] - x[i][1];
+  double delz2 = x[iH2][2] - x[i][2];
+
+  xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2);
+  xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2);
+  xM[2] = x[i][2] + alpha * 0.5 * (delz1 + delz2);
+}

src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h

@@ -0,0 +1,64 @@
+/* -*- 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(dipole/tip4p/chunk,ComputeDipoleTIP4PChunk);
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_DIPOLE_TIP4P_CHUNK_H
+#define LMP_COMPUTE_DIPOLE_TIP4P_CHUNK_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+class Fix;
+
+class ComputeDipoleTIP4PChunk : public Compute {
+ public:
+  ComputeDipoleTIP4PChunk(class LAMMPS *, int, char **);
+  ~ComputeDipoleTIP4PChunk() override;
+  void init() override;
+  void compute_array() override;
+
+  void lock_enable() override;
+  void lock_disable() override;
+  int lock_length() override;
+  void lock(Fix *, bigint, bigint) override;
+  void unlock(Fix *) override;
+
+  double memory_usage() override;
+
+ private:
+  int nchunk, maxchunk;
+  char *idchunk;
+  class ComputeChunkAtom *cchunk;
+
+  double *massproc, *masstotal;
+  double *chrgproc, *chrgtotal;
+  double **com, **comall;
+  double **dipole, **dipoleall;
+  int usecenter;
+
+  int typeO, typeH;
+  double alpha;
+  void find_M(int i, double *xM);
+
+  void allocate();
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif