diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index c2753b2afc..f073212524 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -52,6 +52,8 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`dilatation/atom ` * :doc:`dipole ` * :doc:`dipole/chunk ` + * :doc:`dipole/tip4p ` + * :doc:`dipole/tip4p/chunk ` * :doc:`displace/atom ` * :doc:`dpd ` * :doc:`dpd/atom ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 3d74b3884e..880f60a8a6 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -206,6 +206,8 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`dilatation/atom ` - Peridynamic dilatation for each atom * :doc:`dipole ` - dipole vector and total dipole * :doc:`dipole/chunk ` - dipole vector and total dipole for each chunk +* :doc:`dipole/tip4p ` - dipole vector and total dipole with TIP4P pair style +* :doc:`dipole/tip4p/chunk ` - dipole vector and total dipole for each chunk with TIP4P pair style * :doc:`displace/atom ` - displacement of each atom * :doc:`dpd ` - total values of internal conductive energy, internal mechanical energy, chemical energy, and harmonic average of internal temperature * :doc:`dpd/atom ` - per-particle values of internal conductive energy, internal mechanical energy, chemical energy, and internal temperature diff --git a/doc/src/compute_dipole.rst b/doc/src/compute_dipole.rst index 95c5e216f0..fa1ca4ce0b 100644 --- a/doc/src/compute_dipole.rst +++ b/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 ` 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 -and the total dipole moment for the atoms in the compute group. -This includes all effects due to atoms passing through periodic boundaries. -For a group with a net charge the resulting dipole is made position independent -by subtracting the position vector of the center of mass or geometric center -times the net charge from the computed dipole vector. Both per-atom charges -and per-atom dipole moments, if present, contribute to the computed dipole. +These computes calculate the x,y,z coordinates of the dipole vector and +the total dipole moment for the atoms in the compute group. This +includes all effects due to atoms passing through periodic boundaries. +For a group with a net charge the resulting dipole is made position +independent by subtracting the position vector of the center of mass or +geometric center times the net charge from 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 -the computed dipole moment and a global vector of length 3 with the -dipole vector. See the :doc:`Howto output ` page for -an overview of LAMMPS output options. +These computes calculate a global scalar containing the magnitude of the +computed dipole moment and a global vector of length 3 with the dipole +vector. See the :doc:`Howto output ` page for an overview +of LAMMPS output options. The computed values are "intensive". The array values will be in dipole units (i.e., charge :doc:`units ` times distance @@ -60,7 +72,12 @@ dipole units (i.e., charge :doc:`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 ` page for more info. + +Compute style *dipole/tip4p* can only be used with tip4p pair styles. Related commands """""""""""""""" diff --git a/doc/src/compute_dipole_chunk.rst b/doc/src/compute_dipole_chunk.rst index 504e6f20d0..4e89bb748a 100644 --- a/doc/src/compute_dipole_chunk.rst +++ b/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 ` command -* dipole/chunk = style name of this compute command +* style = *dipole/chunk* or *dipole/tip4p/chunk* * chunkID = ID of :doc:`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 -and the total dipole moment for each chunk, which includes all effects due -to atoms passing through periodic boundaries. For chunks with a net -charge the resulting dipole is made position independent by subtracting -the position vector of the center of mass or geometric center times the -net charge from the computed dipole vector. Both per-atom charges and -per-atom dipole moments, if present, contribute to the computed dipole. +These computes calculate the :math:`(x,y,z)` coordinates of the dipole +vector and the total dipole moment for each chunk, which includes all +effects due to atoms passing through periodic boundaries. For chunks +with a net charge the resulting dipole is made position independent by +subtracting the position vector of the center of mass or geometric +center times the net charge from 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/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 ` 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 ` command. The number of columns is 4 for -the :math:`(x,y,z)` dipole vector components and the total dipole of each -chunk. These values can be accessed by any command that uses global -array values from a compute as input. See the :doc:`Howto output +chunk/atom ` command. The number of columns is 4 +for the :math:`(x,y,z)` dipole vector components and the total dipole of +each chunk. These values can be accessed by any command that uses +global array values from a compute as input. See the :doc:`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 ` 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 ` page for more info. + +Compute style *dipole/tip4p/chunk* can only be used with tip4p pair +styles. Related commands """""""""""""""" diff --git a/src/.gitignore b/src/.gitignore index df02ab3397..0e3a69bdc7 100644 --- a/src/.gitignore +++ b/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 diff --git a/src/EXTRA-COMPUTE/compute_dipole_tip4p.cpp b/src/EXTRA-COMPUTE/compute_dipole_tip4p.cpp new file mode 100644 index 0000000000..c81cd3787c --- /dev/null +++ b/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 +#include + +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); +} diff --git a/src/EXTRA-COMPUTE/compute_dipole_tip4p.h b/src/EXTRA-COMPUTE/compute_dipole_tip4p.h new file mode 100644 index 0000000000..75784b83ec --- /dev/null +++ b/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 diff --git a/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp b/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp new file mode 100644 index 0000000000..92e0d6fab2 --- /dev/null +++ b/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 +#include + +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(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(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); +} diff --git a/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h b/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h new file mode 100644 index 0000000000..63d6e8401e --- /dev/null +++ b/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