diff --git a/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp b/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp index 651811bba5..571f1d562d 100644 --- a/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp +++ b/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -22,11 +21,9 @@ #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 "pair.h" #include #include @@ -39,13 +36,10 @@ 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) + ComputeChunk(lmp, narg, arg), 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/tip4p/chunk command"); + if ((narg != 4) && (narg != 5)) error->all(FLERR, "Illegal compute dipole/tip4p/chunk command"); array_flag = 1; size_array_cols = 4; @@ -53,32 +47,25 @@ ComputeDipoleTIP4PChunk::ComputeDipoleTIP4PChunk(LAMMPS *lmp, int narg, char **a 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/tip4p/chunk command"); + 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/tip4p/chunk command"); } ComputeDipoleTIP4PChunk::init(); - - // chunk-based data - - nchunk = 1; - maxchunk = 0; - allocate(); + ComputeDipoleTIP4PChunk::allocate(); } /* ---------------------------------------------------------------------- */ ComputeDipoleTIP4PChunk::~ComputeDipoleTIP4PChunk() { - delete[] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(chrgproc); @@ -93,14 +80,10 @@ ComputeDipoleTIP4PChunk::~ComputeDipoleTIP4PChunk() void ComputeDipoleTIP4PChunk::init() { - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for compute dipole/tip4p/chunk"); - cchunk = dynamic_cast(modify->compute[icompute]); - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute dipole/tip4p/chunk does not use chunk/atom compute"); + ComputeChunk::init(); - if (!force->pair) error->all(FLERR, "Pair style must be defined for compute dipole/tip4p/chunk"); + if (!force->pair) + error->all(FLERR, "A pair style must be defined for compute dipole/tip4p/chunk"); int itmp; double *p_qdist = (double *) force->pair->extract("qdist", itmp); @@ -130,23 +113,13 @@ void ComputeDipoleTIP4PChunk::init() void ComputeDipoleTIP4PChunk::compute_array() { - int i,index; + 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(); + ComputeChunk::compute_array(); int *ichunk = cchunk->ichunk; - if (nchunk > maxchunk) allocate(); - size_array_rows = nchunk; - // zero local per-chunk values for (i = 0; i < nchunk; i++) { @@ -172,13 +145,16 @@ void ComputeDipoleTIP4PChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + 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); + 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; @@ -186,9 +162,9 @@ void ComputeDipoleTIP4PChunk::compute_array() 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); + 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) { @@ -202,20 +178,20 @@ void ComputeDipoleTIP4PChunk::compute_array() for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; if (type[i] == typeO) { - find_M(i,xM); + find_M(i, xM); xi = xM; } else { xi = x[i]; } - domain->unmap(xi,image[i],unwrap); + 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]; + 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]; @@ -225,82 +201,26 @@ void ComputeDipoleTIP4PChunk::compute_array() } } - MPI_Allreduce(&dipole[0][0],&dipoleall[0][0],4*nchunk,MPI_DOUBLE,MPI_SUM,world); + 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]; + 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])); + 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() { + ComputeChunk::allocate(); memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(chrgproc); @@ -310,14 +230,14 @@ void ComputeDipoleTIP4PChunk::allocate() memory->destroy(dipole); memory->destroy(dipoleall); maxchunk = nchunk; - memory->create(massproc,maxchunk,"dipole/tip4p/chunk:massproc"); - memory->create(masstotal,maxchunk,"dipole/tip4p/chunk:masstotal"); - memory->create(chrgproc,maxchunk,"dipole/tip4p/chunk:chrgproc"); - memory->create(chrgtotal,maxchunk,"dipole/tip4p/chunk:chrgtotal"); - memory->create(com,maxchunk,3,"dipole/tip4p/chunk:com"); - memory->create(comall,maxchunk,3,"dipole/tip4p/chunk:comall"); - memory->create(dipole,maxchunk,4,"dipole/tip4p/chunk:dipole"); - memory->create(dipoleall,maxchunk,4,"dipole/tip4p/chunk:dipoleall"); + memory->create(massproc, maxchunk, "dipole/tip4p/chunk:massproc"); + memory->create(masstotal, maxchunk, "dipole/tip4p/chunk:masstotal"); + memory->create(chrgproc, maxchunk, "dipole/tip4p/chunk:chrgproc"); + memory->create(chrgtotal, maxchunk, "dipole/tip4p/chunk:chrgtotal"); + memory->create(com, maxchunk, 3, "dipole/tip4p/chunk:com"); + memory->create(comall, maxchunk, 3, "dipole/tip4p/chunk:comall"); + memory->create(dipole, maxchunk, 4, "dipole/tip4p/chunk:dipole"); + memory->create(dipoleall, maxchunk, 4, "dipole/tip4p/chunk:dipoleall"); array = dipoleall; } @@ -327,9 +247,9 @@ void ComputeDipoleTIP4PChunk::allocate() 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); + double bytes = (double) maxchunk * 2 * sizeof(double) + ComputeChunk::memory_usage(); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); + bytes += (double) maxchunk * 2 * 4 * sizeof(double); return bytes; } @@ -342,14 +262,14 @@ void ComputeDipoleTIP4PChunk::find_M(int i, double *xM) 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 ((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"); + 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); + 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]; diff --git a/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h b/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h index 63d6e8401e..b3354c9ab9 100644 --- a/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h +++ b/src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h @@ -20,31 +20,21 @@ ComputeStyle(dipole/tip4p/chunk,ComputeDipoleTIP4PChunk); #ifndef LMP_COMPUTE_DIPOLE_TIP4P_CHUNK_H #define LMP_COMPUTE_DIPOLE_TIP4P_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { class Fix; -class ComputeDipoleTIP4PChunk : public Compute { +class ComputeDipoleTIP4PChunk : public ComputeChunk { 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; @@ -55,10 +45,8 @@ class ComputeDipoleTIP4PChunk : public Compute { double alpha; void find_M(int i, double *xM); - void allocate(); + void allocate() override; }; - } // namespace LAMMPS_NS - #endif #endif diff --git a/src/EXTRA-COMPUTE/compute_gyration_shape_chunk.cpp b/src/EXTRA-COMPUTE/compute_gyration_shape_chunk.cpp index 2a4eb78fa2..cc385c4c79 100644 --- a/src/EXTRA-COMPUTE/compute_gyration_shape_chunk.cpp +++ b/src/EXTRA-COMPUTE/compute_gyration_shape_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -33,9 +32,9 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeGyrationShapeChunk::ComputeGyrationShapeChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), id_gyration_chunk(nullptr), shape_parameters(nullptr) + Compute(lmp, narg, arg), id_gyration_chunk(nullptr), shape_parameters(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute gyration/shape/chunk command"); + if (narg != 4) error->all(FLERR, "Illegal compute gyration/shape/chunk command"); // ID of compute gyration id_gyration_chunk = utils::strdup(arg[3]); @@ -58,7 +57,7 @@ ComputeGyrationShapeChunk::ComputeGyrationShapeChunk(LAMMPS *lmp, int narg, char ComputeGyrationShapeChunk::~ComputeGyrationShapeChunk() { - delete [] id_gyration_chunk; + delete[] id_gyration_chunk; memory->destroy(shape_parameters); } @@ -67,22 +66,22 @@ ComputeGyrationShapeChunk::~ComputeGyrationShapeChunk() void ComputeGyrationShapeChunk::init() { // check that the compute gyration command exist - int icompute = modify->find_compute(id_gyration_chunk); - if (icompute < 0) - error->all(FLERR,"Compute gyration/chunk ID does not exist for " - "compute gyration/shape/chunk"); + + c_gyration_chunk = modify->get_compute_by_id(id_gyration_chunk); + if (!c_gyration_chunk) + error->all(FLERR, + "Compute gyration/chunk ID {} does not exist for compute gyration/shape/chunk", + id_gyration_chunk); // check the id_gyration_chunk corresponds really to a compute gyration/chunk command - c_gyration_chunk = (Compute *) modify->compute[icompute]; - if (strcmp(c_gyration_chunk->style,"gyration/chunk") != 0) - error->all(FLERR,"Compute gyration/shape/chunk does not point to " - "gyration compute/chunk"); + + if (strcmp(c_gyration_chunk->style, "gyration/chunk") != 0) + error->all(FLERR, "Compute {} is not of style gyration/chunk", id_gyration_chunk); // check the compute gyration/chunk command computes the whole gyration tensor if (c_gyration_chunk->array_flag == 0) - error->all(FLERR,"Compute gyration/chunk where gyration/shape/chunk points to " - "does not calculate the gyration tensor"); - + error->all(FLERR, "Compute gyration/chunk {} does not calculate the gyration tensor", + id_gyration_chunk); } /* ---------------------------------------------------------------------- */ @@ -108,7 +107,7 @@ void ComputeGyrationShapeChunk::compute_array() invoked_array = update->ntimestep; c_gyration_chunk->compute_array(); - current_nchunks = c_gyration_chunk->size_array_rows; // how to check for the number of chunks in the gyration/chunk? + current_nchunks = c_gyration_chunk->size_array_rows; if (former_nchunks != current_nchunks) allocate(); double **gyration_tensor = c_gyration_chunk->array; @@ -125,18 +124,20 @@ void ComputeGyrationShapeChunk::compute_array() ione[0][2] = ione[2][0] = gyration_tensor[ichunk][4]; ione[1][2] = ione[2][1] = gyration_tensor[ichunk][5]; - int ierror = MathEigen::jacobi3(ione,evalues,evectors); - if (ierror) error->all(FLERR, "Insufficient Jacobi rotations " - "for gyration/shape"); + int ierror = MathEigen::jacobi3(ione, evalues, evectors); + if (ierror) + error->all(FLERR, + "Insufficient Jacobi rotations " + "for gyration/shape"); // sort the eigenvalues according to their size with bubble sort double t; for (int i = 0; i < 3; i++) { - for (int j = 0; j < 2-i; j++) { - if (fabs(evalues[j]) < fabs(evalues[j+1])) { + for (int j = 0; j < 2 - i; j++) { + if (fabs(evalues[j]) < fabs(evalues[j + 1])) { t = evalues[j]; - evalues[j] = evalues[j+1]; - evalues[j+1] = t; + evalues[j] = evalues[j + 1]; + evalues[j + 1] = t; } } } @@ -147,14 +148,14 @@ void ComputeGyrationShapeChunk::compute_array() double sq_eigen_z = MathSpecial::square(evalues[2]); double nominator = sq_eigen_x + sq_eigen_y + sq_eigen_z; - double denominator = MathSpecial::square(evalues[0]+evalues[1]+evalues[2]); + double denominator = MathSpecial::square(evalues[0] + evalues[1] + evalues[2]); shape_parameters[ichunk][0] = evalues[0]; shape_parameters[ichunk][1] = evalues[1]; shape_parameters[ichunk][2] = evalues[2]; - shape_parameters[ichunk][3] = evalues[0] - 0.5*(evalues[1] + evalues[2]); + shape_parameters[ichunk][3] = evalues[0] - 0.5 * (evalues[1] + evalues[2]); shape_parameters[ichunk][4] = evalues[1] - evalues[2]; - shape_parameters[ichunk][5] = 1.5*nominator/denominator - 0.5; + shape_parameters[ichunk][5] = 1.5 * nominator / denominator - 0.5; } } @@ -176,7 +177,7 @@ void ComputeGyrationShapeChunk::allocate() { memory->destroy(shape_parameters); former_nchunks = current_nchunks; - memory->create(shape_parameters,current_nchunks,6,"gyration/shape/chunk:shape_parameters"); + memory->create(shape_parameters, current_nchunks, 6, "gyration/shape/chunk:shape_parameters"); array = shape_parameters; size_array_rows = current_nchunks; } diff --git a/src/compute_angmom_chunk.cpp b/src/compute_angmom_chunk.cpp index 3715bb4def..fcdd20cf07 100644 --- a/src/compute_angmom_chunk.cpp +++ b/src/compute_angmom_chunk.cpp @@ -18,17 +18,13 @@ #include "domain.h" #include "error.h" #include "memory.h" -#include "modify.h" -#include "update.h" - -#include using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeAngmomChunk::ComputeAngmomChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), idchunk(nullptr), massproc(nullptr), masstotal(nullptr), com(nullptr), + ComputeChunk(lmp, narg, arg), massproc(nullptr), masstotal(nullptr), com(nullptr), comall(nullptr), angmom(nullptr), angmomall(nullptr) { if (narg != 4) error->all(FLERR, "Illegal compute angmom/chunk command"); @@ -39,24 +35,14 @@ ComputeAngmomChunk::ComputeAngmomChunk(LAMMPS *lmp, int narg, char **arg) : size_array_rows_variable = 1; extarray = 0; - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); - ComputeAngmomChunk::init(); - - // chunk-based data - - nchunk = 1; - maxchunk = 0; - allocate(); + ComputeAngmomChunk::allocate(); } /* ---------------------------------------------------------------------- */ ComputeAngmomChunk::~ComputeAngmomChunk() { - delete[] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -67,35 +53,15 @@ ComputeAngmomChunk::~ComputeAngmomChunk() /* ---------------------------------------------------------------------- */ -void ComputeAngmomChunk::init() -{ - cchunk = dynamic_cast(modify->get_compute_by_id(idchunk)); - if (!cchunk) error->all(FLERR, "Chunk/atom compute does not exist for compute angmom/chunk"); - if (strcmp(cchunk->style, "chunk/atom") != 0) - error->all(FLERR, "Compute angmom/chunk does not use chunk/atom compute"); -} - -/* ---------------------------------------------------------------------- */ - void ComputeAngmomChunk::compute_array() { + ComputeChunk::compute_array(); + int i, index; double dx, dy, dz, 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++) { @@ -164,59 +130,6 @@ void ComputeAngmomChunk::compute_array() MPI_Allreduce(&angmom[0][0], &angmomall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); } -/* ---------------------------------------------------------------------- - 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 ComputeAngmomChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputeAngmomChunk::lock_disable() -{ - cchunk = dynamic_cast(modify->get_compute_by_id(idchunk)); - if (cchunk) cchunk->lockcount--; -} - -/* ---------------------------------------------------------------------- - calculate and return # of chunks = length of vector/array -------------------------------------------------------------------------- */ - -int ComputeAngmomChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeAngmomChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr, startstep, stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeAngmomChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ @@ -245,7 +158,8 @@ void ComputeAngmomChunk::allocate() double ComputeAngmomChunk::memory_usage() { - double bytes = (bigint) maxchunk * 2 * sizeof(double); + double bytes = ComputeChunk::memory_usage(); + bytes += (bigint) maxchunk * 2 * sizeof(double); bytes += (double) maxchunk * 2 * 3 * sizeof(double); bytes += (double) maxchunk * 2 * 3 * sizeof(double); return bytes; diff --git a/src/compute_angmom_chunk.h b/src/compute_angmom_chunk.h index 04023dbe44..7f8172f863 100644 --- a/src/compute_angmom_chunk.h +++ b/src/compute_angmom_chunk.h @@ -20,39 +20,25 @@ ComputeStyle(angmom/chunk,ComputeAngmomChunk); #ifndef LMP_COMPUTE_ANGMOM_CHUNK_H #define LMP_COMPUTE_ANGMOM_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { - class Fix; -class ComputeAngmomChunk : public Compute { +class ComputeAngmomChunk : public ComputeChunk { public: ComputeAngmomChunk(class LAMMPS *, int, char **); ~ComputeAngmomChunk() 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 **com, **comall; double **angmom, **angmomall; - void allocate(); + void allocate() override; }; - } // namespace LAMMPS_NS - #endif #endif diff --git a/src/compute_chunk.cpp b/src/compute_chunk.cpp new file mode 100644 index 0000000000..0b5204ff94 --- /dev/null +++ b/src/compute_chunk.cpp @@ -0,0 +1,156 @@ +/* ---------------------------------------------------------------------- + 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_chunk.h" + +#include "compute_chunk_atom.h" +#include "error.h" +#include "modify.h" +#include "update.h" + +using namespace LAMMPS_NS; +/* ---------------------------------------------------------------------- */ + +ComputeChunk::ComputeChunk(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), idchunk(nullptr), cchunk(nullptr) +{ + if (narg < 4) utils::missing_cmd_args(FLERR, std::string("compute ") + style, error); + + // ID of compute chunk/atom + + idchunk = utils::strdup(arg[3]); + + ComputeChunk::init(); + + // chunk-based data + + nchunk = 1; + maxchunk = 0; + firstflag = 1; + massneed = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeChunk::~ComputeChunk() +{ + delete[] idchunk; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeChunk::init() +{ + cchunk = dynamic_cast(modify->get_compute_by_id(idchunk)); + if (!cchunk) + error->all(FLERR, "Chunk/atom compute {} does not exist or is incorrect style for compute {}", + idchunk, style); +} + +/* ---------------------------------------------------------------------- + common code used by all /chunk computes +------------------------------------------------------------------------- */ + +void ComputeChunk::compute_vector() +{ + invoked_vector = 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(); + + if (nchunk > maxchunk) allocate(); + size_vector = nchunk; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeChunk::compute_array() +{ + 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(); + + if (nchunk > maxchunk) allocate(); + size_array_rows = nchunk; +} + +/* ---------------------------------------------------------------------- + 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 ComputeChunk::lock_enable() +{ + cchunk->lockcount++; +} + +/* ---------------------------------------------------------------------- + decrement lock counter in compute chunk/atom, it if still exists +------------------------------------------------------------------------- */ + +void ComputeChunk::lock_disable() +{ + cchunk = dynamic_cast(modify->get_compute_by_id(idchunk)); + if (!cchunk) cchunk->lockcount--; +} + +/* ---------------------------------------------------------------------- + calculate and return # of chunks = length of vector/array +------------------------------------------------------------------------- */ + +int ComputeChunk::lock_length() +{ + nchunk = cchunk->setup_chunks(); + return nchunk; +} + +/* ---------------------------------------------------------------------- + set the lock from startstep to stopstep +------------------------------------------------------------------------- */ + +void ComputeChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) +{ + cchunk->lock(fixptr, startstep, stopstep); +} + +/* ---------------------------------------------------------------------- + unset the lock +------------------------------------------------------------------------- */ + +void ComputeChunk::unlock(Fix *fixptr) +{ + cchunk->unlock(fixptr); +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeChunk::memory_usage() +{ + return sizeof(ComputeChunk) + sizeof(ComputeChunkAtom); +} diff --git a/src/compute_chunk.h b/src/compute_chunk.h new file mode 100644 index 0000000000..0c94f00bb0 --- /dev/null +++ b/src/compute_chunk.h @@ -0,0 +1,49 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +#ifndef LMP_COMPUTE_CHUNK_H +#define LMP_COMPUTE_CHUNK_H + +#include "compute.h" + +namespace LAMMPS_NS { +class ComputeChunkAtom; +class Fix; + +class ComputeChunk : public Compute { + public: + char *idchunk; // fields accessed by other classes + + ComputeChunk(class LAMMPS *, int, char **); + ~ComputeChunk() override; + void init() override; + void compute_vector() 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; + + protected: + int nchunk, maxchunk; + int firstflag, massneed; + ComputeChunkAtom *cchunk; + + virtual void allocate(){}; +}; +} // namespace LAMMPS_NS +#endif diff --git a/src/compute_com_chunk.cpp b/src/compute_com_chunk.cpp index 7cd80545ee..582bd4f6fe 100644 --- a/src/compute_com_chunk.cpp +++ b/src/compute_com_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -19,22 +18,18 @@ #include "domain.h" #include "error.h" #include "memory.h" -#include "modify.h" -#include "update.h" - -#include using namespace LAMMPS_NS; -enum{ONCE,NFREQ,EVERY}; +enum { ONCE, NFREQ, EVERY }; /* ---------------------------------------------------------------------- */ ComputeCOMChunk::ComputeCOMChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - idchunk(nullptr), masstotal(nullptr), massproc(nullptr), com(nullptr), comall(nullptr) + ComputeChunk(lmp, narg, arg), masstotal(nullptr), massproc(nullptr), com(nullptr), + comall(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute com/chunk command"); + if (narg != 4) error->all(FLERR, "Illegal compute com/chunk command"); array_flag = 1; size_array_cols = 3; @@ -42,26 +37,14 @@ ComputeCOMChunk::ComputeCOMChunk(LAMMPS *lmp, int narg, char **arg) : size_array_rows_variable = 1; extarray = 0; - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); - ComputeCOMChunk::init(); - - // chunk-based data - - nchunk = 1; - maxchunk = 0; - allocate(); - - firstflag = massneed = 1; + ComputeCOMChunk::allocate(); } /* ---------------------------------------------------------------------- */ ComputeCOMChunk::~ComputeCOMChunk() { - delete [] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -70,18 +53,6 @@ ComputeCOMChunk::~ComputeCOMChunk() /* ---------------------------------------------------------------------- */ -void ComputeCOMChunk::init() -{ - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for compute com/chunk"); - cchunk = dynamic_cast(modify->compute[icompute]); - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute com/chunk does not use chunk/atom compute"); -} - -/* ---------------------------------------------------------------------- */ - void ComputeCOMChunk::setup() { // one-time calculation of per-chunk mass @@ -101,23 +72,12 @@ void ComputeCOMChunk::compute_array() 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(); + ComputeChunk::compute_array(); int *ichunk = cchunk->ichunk; - if (nchunk > maxchunk) allocate(); - size_array_rows = nchunk; - // zero local per-chunk values - for (int i = 0; i < nchunk; i++) - com[i][0] = com[i][1] = com[i][2] = 0.0; + for (int i = 0; i < nchunk; i++) com[i][0] = com[i][1] = com[i][2] = 0.0; if (massneed) for (int i = 0; i < nchunk; i++) massproc[i] = 0.0; @@ -133,86 +93,32 @@ void ComputeCOMChunk::compute_array() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - domain->unmap(x[i],image[i],unwrap); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + domain->unmap(x[i], image[i], unwrap); com[index][0] += unwrap[0] * massone; com[index][1] += unwrap[1] * massone; com[index][2] += unwrap[2] * massone; if (massneed) massproc[index] += massone; } - MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); - if (massneed) - MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&com[0][0], &comall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); + if (massneed) MPI_Allreduce(massproc, masstotal, nchunk, MPI_DOUBLE, MPI_SUM, world); for (int 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]; - } else comall[i][0] = comall[i][1] = comall[i][2] = 0.0; + } else + comall[i][0] = comall[i][1] = comall[i][2] = 0.0; } } -/* ---------------------------------------------------------------------- - 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 ComputeCOMChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputeCOMChunk::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 ComputeCOMChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeCOMChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeCOMChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ @@ -224,10 +130,10 @@ void ComputeCOMChunk::allocate() memory->destroy(com); memory->destroy(comall); maxchunk = nchunk; - memory->create(massproc,maxchunk,"com/chunk:massproc"); - memory->create(masstotal,maxchunk,"com/chunk:masstotal"); - memory->create(com,maxchunk,3,"com/chunk:com"); - memory->create(comall,maxchunk,3,"com/chunk:comall"); + memory->create(massproc, maxchunk, "com/chunk:massproc"); + memory->create(masstotal, maxchunk, "com/chunk:masstotal"); + memory->create(com, maxchunk, 3, "com/chunk:com"); + memory->create(comall, maxchunk, 3, "com/chunk:comall"); array = comall; } @@ -237,7 +143,8 @@ void ComputeCOMChunk::allocate() double ComputeCOMChunk::memory_usage() { - double bytes = (bigint) maxchunk * 2 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); + double bytes = ComputeChunk::memory_usage(); + bytes += (bigint) maxchunk * 2 * sizeof(double); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); return bytes; } diff --git a/src/compute_com_chunk.h b/src/compute_com_chunk.h index c617e181fe..17731eb4b5 100644 --- a/src/compute_com_chunk.h +++ b/src/compute_com_chunk.h @@ -20,42 +20,26 @@ ComputeStyle(com/chunk,ComputeCOMChunk); #ifndef LMP_COMPUTE_COM_CHUNK_H #define LMP_COMPUTE_COM_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { -class Fix; - -class ComputeCOMChunk : public Compute { +class ComputeCOMChunk : public ComputeChunk { public: - char *idchunk; // fields accessed by other classes double *masstotal; ComputeCOMChunk(class LAMMPS *, int, char **); ~ComputeCOMChunk() override; - void init() override; void setup() 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; - int firstflag, massneed; - class ComputeChunkAtom *cchunk; - double *massproc; double **com, **comall; - void allocate(); + void allocate() override; }; - } // namespace LAMMPS_NS - #endif #endif diff --git a/src/compute_dipole_chunk.cpp b/src/compute_dipole_chunk.cpp index 582a9f8652..8068ac583c 100644 --- a/src/compute_dipole_chunk.cpp +++ b/src/compute_dipole_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -22,8 +21,6 @@ #include "force.h" #include "math_special.h" #include "memory.h" -#include "modify.h" -#include "update.h" #include #include @@ -36,13 +33,10 @@ enum { MASSCENTER, GEOMCENTER }; /* ---------------------------------------------------------------------- */ ComputeDipoleChunk::ComputeDipoleChunk(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) + ComputeChunk(lmp, narg, arg), 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"); + if ((narg != 4) && (narg != 5)) error->all(FLERR, "Illegal compute dipole/chunk command"); array_flag = 1; size_array_cols = 4; @@ -50,32 +44,25 @@ ComputeDipoleChunk::ComputeDipoleChunk(LAMMPS *lmp, int narg, char **arg) : 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"); + 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"); } ComputeDipoleChunk::init(); - - // chunk-based data - - nchunk = 1; - maxchunk = 0; - allocate(); + ComputeDipoleChunk::allocate(); } /* ---------------------------------------------------------------------- */ ComputeDipoleChunk::~ComputeDipoleChunk() { - delete [] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(chrgproc); @@ -90,40 +77,23 @@ ComputeDipoleChunk::~ComputeDipoleChunk() void ComputeDipoleChunk::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"); + ComputeChunk::init(); - if ((force->pair_match("/tip4p/",0) != nullptr) && (comm->me == 0)) - error->warning(FLERR,"Computed dipole moments may be incorrect when " - "using a tip4p pair style"); + if ((force->pair_match("tip4p/", 0) != nullptr) && (comm->me == 0)) + error->warning(FLERR, "Dipole moments may be incorrect when sing a TIP4P pair style"); } /* ---------------------------------------------------------------------- */ void ComputeDipoleChunk::compute_array() { - int i,index; + 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(); + ComputeChunk::compute_array(); int *ichunk = cchunk->ichunk; - if (nchunk > maxchunk) allocate(); - size_array_rows = nchunk; - // zero local per-chunk values for (i = 0; i < nchunk; i++) { @@ -147,14 +117,17 @@ void ComputeDipoleChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + 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 + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + } else + massone = 1.0; // usecenter == GEOMCENTER - domain->unmap(x[i],image[i],unwrap); + domain->unmap(x[i], image[i], unwrap); massproc[index] += massone; if (atom->q_flag) chrgproc[index] += atom->q[i]; com[index][0] += unwrap[0] * massone; @@ -162,9 +135,9 @@ void ComputeDipoleChunk::compute_array() 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); + 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) { @@ -178,13 +151,13 @@ void ComputeDipoleChunk::compute_array() for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - domain->unmap(x[i],image[i],unwrap); + domain->unmap(x[i], 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]; + 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]; @@ -194,83 +167,25 @@ void ComputeDipoleChunk::compute_array() } } - MPI_Allreduce(&dipole[0][0],&dipoleall[0][0],4*nchunk, - MPI_DOUBLE,MPI_SUM,world); + 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]; + 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])); + 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 ComputeDipoleChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputeDipoleChunk::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 ComputeDipoleChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeDipoleChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeDipoleChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ void ComputeDipoleChunk::allocate() { + ComputeChunk::allocate(); memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(chrgproc); @@ -280,14 +195,14 @@ void ComputeDipoleChunk::allocate() 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"); + 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; } @@ -297,8 +212,9 @@ void ComputeDipoleChunk::allocate() double ComputeDipoleChunk::memory_usage() { - double bytes = (bigint) maxchunk * 2 * sizeof(double); - bytes += (double)maxchunk * 2*3 * sizeof(double); - bytes += (double)maxchunk * 2*4 * sizeof(double); + double bytes = ComputeChunk::memory_usage(); + bytes += (bigint) maxchunk * 2 * sizeof(double); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); + bytes += (double) maxchunk * 2 * 4 * sizeof(double); return bytes; } diff --git a/src/compute_dipole_chunk.h b/src/compute_dipole_chunk.h index 42198fe478..603e6a4353 100644 --- a/src/compute_dipole_chunk.h +++ b/src/compute_dipole_chunk.h @@ -20,38 +20,28 @@ ComputeStyle(dipole/chunk,ComputeDipoleChunk); #ifndef LMP_COMPUTE_DIPOLE_CHUNK_H #define LMP_COMPUTE_DIPOLE_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { class Fix; -class ComputeDipoleChunk : public Compute { +class ComputeDipoleChunk : public ComputeChunk { public: ComputeDipoleChunk(class LAMMPS *, int, char **); ~ComputeDipoleChunk() 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; - + protected: double *massproc, *masstotal; double *chrgproc, *chrgtotal; double **com, **comall; double **dipole, **dipoleall; int usecenter; - void allocate(); + void allocate() override; }; } // namespace LAMMPS_NS diff --git a/src/compute_gyration_chunk.cpp b/src/compute_gyration_chunk.cpp index 4492d81d18..a02eb7fa4e 100644 --- a/src/compute_gyration_chunk.cpp +++ b/src/compute_gyration_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,31 +13,23 @@ #include "compute_gyration_chunk.h" -#include -#include #include "atom.h" -#include "update.h" -#include "modify.h" #include "compute_chunk_atom.h" #include "domain.h" -#include "memory.h" #include "error.h" +#include "memory.h" + +#include +#include using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeGyrationChunk::ComputeGyrationChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - idchunk(nullptr), massproc(nullptr), masstotal(nullptr), com(nullptr), comall(nullptr), - rg(nullptr), rgall(nullptr), rgt(nullptr), rgtall(nullptr) + ComputeChunk(lmp, narg, arg), massproc(nullptr), masstotal(nullptr), com(nullptr), + comall(nullptr), rg(nullptr), rgall(nullptr), rgt(nullptr), rgtall(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute gyration/chunk command"); - - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); - ComputeGyrationChunk::init(); // optional args @@ -46,10 +37,11 @@ ComputeGyrationChunk::ComputeGyrationChunk(LAMMPS *lmp, int narg, char **arg) : tensor = 0; int iarg = 4; while (iarg < narg) { - if (strcmp(arg[iarg],"tensor") == 0) { + if (strcmp(arg[iarg], "tensor") == 0) { tensor = 1; iarg++; - } else error->all(FLERR,"Illegal compute gyration/chunk command"); + } else + error->all(FLERR, "Illegal compute gyration/chunk command"); } if (tensor) { @@ -65,18 +57,13 @@ ComputeGyrationChunk::ComputeGyrationChunk(LAMMPS *lmp, int narg, char **arg) : extvector = 0; } - // chunk-based data - - nchunk = 1; - maxchunk = 0; - allocate(); + ComputeGyrationChunk::allocate(); } /* ---------------------------------------------------------------------- */ ComputeGyrationChunk::~ComputeGyrationChunk() { - delete [] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -89,26 +76,13 @@ ComputeGyrationChunk::~ComputeGyrationChunk() /* ---------------------------------------------------------------------- */ -void ComputeGyrationChunk::init() -{ - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for " - "compute gyration/chunk"); - cchunk = dynamic_cast(modify->compute[icompute]); - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute gyration/chunk does not use chunk/atom compute"); -} - -/* ---------------------------------------------------------------------- */ - void ComputeGyrationChunk::compute_vector() { - int i,index; - double dx,dy,dz,massone; + int i, index; + double dx, dy, dz, massone; double unwrap[3]; - invoked_array = update->ntimestep; + ComputeChunk::compute_vector(); com_chunk(); int *ichunk = cchunk->ichunk; @@ -127,33 +101,34 @@ void ComputeGyrationChunk::compute_vector() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - domain->unmap(x[i],image[i],unwrap); + domain->unmap(x[i], image[i], unwrap); dx = unwrap[0] - comall[index][0]; dy = unwrap[1] - comall[index][1]; dz = unwrap[2] - comall[index][2]; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - rg[index] += (dx*dx + dy*dy + dz*dz) * massone; + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + rg[index] += (dx * dx + dy * dy + dz * dz) * massone; } - MPI_Allreduce(rg,rgall,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(rg, rgall, nchunk, MPI_DOUBLE, MPI_SUM, world); for (i = 0; i < nchunk; i++) - if (masstotal[i] > 0.0) - rgall[i] = sqrt(rgall[i]/masstotal[i]); + if (masstotal[i] > 0.0) rgall[i] = sqrt(rgall[i] / masstotal[i]); } /* ---------------------------------------------------------------------- */ void ComputeGyrationChunk::compute_array() { - int i,j,index; - double dx,dy,dz,massone; + int i, j, index; + double dx, dy, dz, massone; double unwrap[3]; - invoked_array = update->ntimestep; + ComputeChunk::compute_array(); com_chunk(); int *ichunk = cchunk->ichunk; @@ -171,34 +146,33 @@ void ComputeGyrationChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - domain->unmap(x[i],image[i],unwrap); + domain->unmap(x[i], image[i], unwrap); dx = unwrap[0] - comall[index][0]; dy = unwrap[1] - comall[index][1]; dz = unwrap[2] - comall[index][2]; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - rgt[index][0] += dx*dx * massone; - rgt[index][1] += dy*dy * massone; - rgt[index][2] += dz*dz * massone; - rgt[index][3] += dx*dy * massone; - rgt[index][4] += dx*dz * massone; - rgt[index][5] += dy*dz * massone; + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + rgt[index][0] += dx * dx * massone; + rgt[index][1] += dy * dy * massone; + rgt[index][2] += dz * dz * massone; + rgt[index][3] += dx * dy * massone; + rgt[index][4] += dx * dz * massone; + rgt[index][5] += dy * dz * massone; } - if (nchunk) - MPI_Allreduce(&rgt[0][0],&rgtall[0][0],nchunk*6,MPI_DOUBLE,MPI_SUM,world); + if (nchunk) MPI_Allreduce(&rgt[0][0], &rgtall[0][0], nchunk * 6, MPI_DOUBLE, MPI_SUM, world); for (i = 0; i < nchunk; i++) { if (masstotal[i] > 0.0) { - for (j = 0; j < 6; j++) - rgtall[i][j] = rgtall[i][j]/masstotal[i]; + for (j = 0; j < 6; j++) rgtall[i][j] = rgtall[i][j] / masstotal[i]; } } } - /* ---------------------------------------------------------------------- calculate per-chunk COM, used by both scalar and tensor ------------------------------------------------------------------------- */ @@ -209,18 +183,8 @@ void ComputeGyrationChunk::com_chunk() double massone; double unwrap[3]; - // 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(); - if (tensor) size_array_rows = nchunk; - else size_vector = nchunk; - // zero local per-chunk values for (int i = 0; i < nchunk; i++) { @@ -240,19 +204,21 @@ void ComputeGyrationChunk::com_chunk() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - domain->unmap(x[i],image[i],unwrap); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + domain->unmap(x[i], image[i], unwrap); massproc[index] += massone; 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(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(massproc, masstotal, nchunk, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(&com[0][0], &comall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); for (int i = 0; i < nchunk; i++) { if (masstotal[i] > 0.0) { @@ -263,68 +229,13 @@ void ComputeGyrationChunk::com_chunk() } } -/* ---------------------------------------------------------------------- - 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 ComputeGyrationChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputeGyrationChunk::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 ComputeGyrationChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeGyrationChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeGyrationChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ void ComputeGyrationChunk::allocate() { + ComputeChunk::allocate(); memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -334,17 +245,17 @@ void ComputeGyrationChunk::allocate() memory->destroy(rgt); memory->destroy(rgtall); maxchunk = nchunk; - memory->create(massproc,maxchunk,"gyration/chunk:massproc"); - memory->create(masstotal,maxchunk,"gyration/chunk:masstotal"); - memory->create(com,maxchunk,3,"gyration/chunk:com"); - memory->create(comall,maxchunk,3,"gyration/chunk:comall"); + memory->create(massproc, maxchunk, "gyration/chunk:massproc"); + memory->create(masstotal, maxchunk, "gyration/chunk:masstotal"); + memory->create(com, maxchunk, 3, "gyration/chunk:com"); + memory->create(comall, maxchunk, 3, "gyration/chunk:comall"); if (tensor) { - memory->create(rgt,maxchunk,6,"gyration/chunk:rgt"); - memory->create(rgtall,maxchunk,6,"gyration/chunk:rgtall"); + memory->create(rgt, maxchunk, 6, "gyration/chunk:rgt"); + memory->create(rgtall, maxchunk, 6, "gyration/chunk:rgtall"); array = rgtall; } else { - memory->create(rg,maxchunk,"gyration/chunk:rg"); - memory->create(rgall,maxchunk,"gyration/chunk:rgall"); + memory->create(rg, maxchunk, "gyration/chunk:rg"); + memory->create(rgall, maxchunk, "gyration/chunk:rgall"); vector = rgall; } } @@ -355,9 +266,12 @@ void ComputeGyrationChunk::allocate() double ComputeGyrationChunk::memory_usage() { - double bytes = (bigint) maxchunk * 2 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); - if (tensor) bytes += (double) maxchunk * 2*6 * sizeof(double); - else bytes += (double) maxchunk * 2 * sizeof(double); + double bytes = ComputeChunk::memory_usage(); + bytes += (bigint) maxchunk * 2 * sizeof(double); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); + if (tensor) + bytes += (double) maxchunk * 2 * 6 * sizeof(double); + else + bytes += (double) maxchunk * 2 * sizeof(double); return bytes; } diff --git a/src/compute_gyration_chunk.h b/src/compute_gyration_chunk.h index c455f5a100..640d350a05 100644 --- a/src/compute_gyration_chunk.h +++ b/src/compute_gyration_chunk.h @@ -20,31 +20,21 @@ ComputeStyle(gyration/chunk,ComputeGyrationChunk); #ifndef LMP_COMPUTE_GYRATION_CHUNK_H #define LMP_COMPUTE_GYRATION_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { -class ComputeGyrationChunk : public Compute { +class ComputeGyrationChunk : public ComputeChunk { public: ComputeGyrationChunk(class LAMMPS *, int, char **); ~ComputeGyrationChunk() override; - void init() override; + void compute_vector() override; void compute_array() override; - void lock_enable() override; - void lock_disable() override; - int lock_length() override; - void lock(class Fix *, bigint, bigint) override; - void unlock(class Fix *) override; - double memory_usage() override; private: - int nchunk, maxchunk; - char *idchunk; - class ComputeChunkAtom *cchunk; - int tensor; double *massproc, *masstotal; @@ -53,7 +43,7 @@ class ComputeGyrationChunk : public Compute { double **rgt, **rgtall; void com_chunk(); - void allocate(); + void allocate() override; }; } // namespace LAMMPS_NS diff --git a/src/compute_inertia_chunk.cpp b/src/compute_inertia_chunk.cpp index 1e8f331fe1..0a0507c0cf 100644 --- a/src/compute_inertia_chunk.cpp +++ b/src/compute_inertia_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,25 +13,21 @@ #include "compute_inertia_chunk.h" -#include #include "atom.h" -#include "update.h" -#include "modify.h" #include "compute_chunk_atom.h" #include "domain.h" -#include "memory.h" #include "error.h" +#include "memory.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeInertiaChunk::ComputeInertiaChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - idchunk(nullptr), massproc(nullptr), masstotal(nullptr), com(nullptr), comall(nullptr), - inertia(nullptr), inertiaall(nullptr) + ComputeChunk(lmp, narg, arg), massproc(nullptr), masstotal(nullptr), com(nullptr), + comall(nullptr), inertia(nullptr), inertiaall(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute inertia/chunk command"); + if (narg != 4) error->all(FLERR, "Illegal compute inertia/chunk command"); array_flag = 1; size_array_cols = 6; @@ -40,24 +35,14 @@ ComputeInertiaChunk::ComputeInertiaChunk(LAMMPS *lmp, int narg, char **arg) : size_array_rows_variable = 1; extarray = 0; - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); - ComputeInertiaChunk::init(); - - // chunk-based data - - nchunk = 1; - maxchunk = 0; - allocate(); + ComputeInertiaChunk::allocate(); } /* ---------------------------------------------------------------------- */ ComputeInertiaChunk::~ComputeInertiaChunk() { - delete [] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -68,38 +53,15 @@ ComputeInertiaChunk::~ComputeInertiaChunk() /* ---------------------------------------------------------------------- */ -void ComputeInertiaChunk::init() -{ - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for " - "compute inertia/chunk"); - cchunk = dynamic_cast(modify->compute[icompute]); - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute inertia/chunk does not use chunk/atom compute"); -} - -/* ---------------------------------------------------------------------- */ - void ComputeInertiaChunk::compute_array() { - int i,j,index; - double dx,dy,dz,massone; + int i, j, index; + double dx, dy, dz, 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(); + ComputeChunk::compute_array(); int *ichunk = cchunk->ichunk; - if (nchunk > maxchunk) allocate(); - size_array_rows = nchunk; - // zero local per-chunk values for (i = 0; i < nchunk; i++) { @@ -120,19 +82,21 @@ void ComputeInertiaChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - domain->unmap(x[i],image[i],unwrap); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + domain->unmap(x[i], image[i], unwrap); massproc[index] += massone; 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(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(massproc, masstotal, 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) { @@ -146,83 +110,26 @@ void ComputeInertiaChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - domain->unmap(x[i],image[i],unwrap); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + domain->unmap(x[i], image[i], unwrap); dx = unwrap[0] - comall[index][0]; dy = unwrap[1] - comall[index][1]; dz = unwrap[2] - comall[index][2]; - inertia[index][0] += massone * (dy*dy + dz*dz); - inertia[index][1] += massone * (dx*dx + dz*dz); - inertia[index][2] += massone * (dx*dx + dy*dy); - inertia[index][3] -= massone * dx*dy; - inertia[index][4] -= massone * dy*dz; - inertia[index][5] -= massone * dx*dz; + inertia[index][0] += massone * (dy * dy + dz * dz); + inertia[index][1] += massone * (dx * dx + dz * dz); + inertia[index][2] += massone * (dx * dx + dy * dy); + inertia[index][3] -= massone * dx * dy; + inertia[index][4] -= massone * dy * dz; + inertia[index][5] -= massone * dx * dz; } - MPI_Allreduce(&inertia[0][0],&inertiaall[0][0],6*nchunk, - MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&inertia[0][0], &inertiaall[0][0], 6 * nchunk, MPI_DOUBLE, MPI_SUM, world); } - -/* ---------------------------------------------------------------------- - 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 ComputeInertiaChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputeInertiaChunk::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 ComputeInertiaChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeInertiaChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeInertiaChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - - /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ @@ -236,12 +143,12 @@ void ComputeInertiaChunk::allocate() memory->destroy(inertia); memory->destroy(inertiaall); maxchunk = nchunk; - memory->create(massproc,maxchunk,"inertia/chunk:massproc"); - memory->create(masstotal,maxchunk,"inertia/chunk:masstotal"); - memory->create(com,maxchunk,3,"inertia/chunk:com"); - memory->create(comall,maxchunk,3,"inertia/chunk:comall"); - memory->create(inertia,maxchunk,6,"inertia/chunk:inertia"); - memory->create(inertiaall,maxchunk,6,"inertia/chunk:inertiaall"); + memory->create(massproc, maxchunk, "inertia/chunk:massproc"); + memory->create(masstotal, maxchunk, "inertia/chunk:masstotal"); + memory->create(com, maxchunk, 3, "inertia/chunk:com"); + memory->create(comall, maxchunk, 3, "inertia/chunk:comall"); + memory->create(inertia, maxchunk, 6, "inertia/chunk:inertia"); + memory->create(inertiaall, maxchunk, 6, "inertia/chunk:inertiaall"); array = inertiaall; } @@ -251,8 +158,9 @@ void ComputeInertiaChunk::allocate() double ComputeInertiaChunk::memory_usage() { - double bytes = (bigint) maxchunk * 2 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); - bytes += (double) maxchunk * 2*6 * sizeof(double); + double bytes = ComputeChunk::memory_usage(); + bytes += (bigint) maxchunk * 2 * sizeof(double); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); + bytes += (double) maxchunk * 2 * 6 * sizeof(double); return bytes; } diff --git a/src/compute_inertia_chunk.h b/src/compute_inertia_chunk.h index 6ec318833a..17acd4ac40 100644 --- a/src/compute_inertia_chunk.h +++ b/src/compute_inertia_chunk.h @@ -20,38 +20,25 @@ ComputeStyle(inertia/chunk,ComputeInertiaChunk); #ifndef LMP_COMPUTE_INERTIA_CHUNK_H #define LMP_COMPUTE_INERTIA_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { -class ComputeInertiaChunk : public Compute { +class ComputeInertiaChunk : public ComputeChunk { public: ComputeInertiaChunk(class LAMMPS *, int, char **); ~ComputeInertiaChunk() override; - void init() override; void compute_array() override; - void lock_enable() override; - void lock_disable() override; - int lock_length() override; - void lock(class Fix *, bigint, bigint) override; - void unlock(class Fix *) override; - double memory_usage() override; private: - int nchunk, maxchunk; - char *idchunk; - class ComputeChunkAtom *cchunk; - double *massproc, *masstotal; double **com, **comall; double **inertia, **inertiaall; - void allocate(); + void allocate() override; }; - } // namespace LAMMPS_NS - #endif #endif diff --git a/src/compute_msd_chunk.cpp b/src/compute_msd_chunk.cpp index 6e23ac1a86..07234ecfdb 100644 --- a/src/compute_msd_chunk.cpp +++ b/src/compute_msd_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -22,19 +21,16 @@ #include "group.h" #include "memory.h" #include "modify.h" -#include "update.h" - -#include using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeMSDChunk::ComputeMSDChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), idchunk(nullptr), id_fix(nullptr), massproc(nullptr), - masstotal(nullptr), com(nullptr), comall(nullptr), msd(nullptr) + ComputeChunk(lmp, narg, arg), id_fix(nullptr), massproc(nullptr), masstotal(nullptr), + com(nullptr), comall(nullptr), msd(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute msd/chunk command"); + if (narg != 4) error->all(FLERR, "Illegal compute msd/chunk command"); array_flag = 1; size_array_cols = 4; @@ -42,11 +38,6 @@ ComputeMSDChunk::ComputeMSDChunk(LAMMPS *lmp, int narg, char **arg) : size_array_rows_variable = 1; extarray = 0; - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); - - firstflag = 1; ComputeMSDChunk::init(); // create a new fix STORE style for reference positions @@ -58,7 +49,7 @@ ComputeMSDChunk::ComputeMSDChunk(LAMMPS *lmp, int narg, char **arg) : id_fix = utils::strdup(std::string(id) + "_COMPUTE_STORE"); fix = dynamic_cast( - modify->add_fix(fmt::format("{} {} STORE/GLOBAL 1 1", id_fix,group->names[igroup]))); + modify->add_fix(fmt::format("{} {} STORE/GLOBAL 1 1", id_fix, group->names[igroup]))); } /* ---------------------------------------------------------------------- */ @@ -70,7 +61,6 @@ ComputeMSDChunk::~ComputeMSDChunk() if (modify->nfix) modify->delete_fix(id_fix); delete[] id_fix; - delete[] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -82,19 +72,14 @@ ComputeMSDChunk::~ComputeMSDChunk() void ComputeMSDChunk::init() { - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for compute msd/chunk"); - cchunk = dynamic_cast(modify->compute[icompute]); - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute msd/chunk does not use chunk/atom compute"); + ComputeChunk::init(); // set fix which stores reference atom coords // if firstflag, will be created in setup() if (!firstflag) { fix = dynamic_cast(modify->get_fix_by_id(id_fix)); - if (!fix) error->all(FLERR,"Could not find compute msd/chunk fix with ID {}", id_fix); + if (!fix) error->all(FLERR, "Could not find compute msd/chunk fix with ID {}", id_fix); } } @@ -113,7 +98,7 @@ void ComputeMSDChunk::setup() // otherwise reset its size now and initialize to current COM if (fix->nrow == nchunk && fix->ncol == 3) return; - fix->reset_global(nchunk,3); + fix->reset_global(nchunk, 3); double **cominit = fix->astore; for (int i = 0; i < nchunk; i++) { @@ -132,25 +117,15 @@ void ComputeMSDChunk::compute_array() 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 - - int n = cchunk->setup_chunks(); - cchunk->compute_ichunk(); + int oldnchunk = nchunk; + ComputeChunk::compute_array(); int *ichunk = cchunk->ichunk; // first time call, allocate per-chunk arrays // thereafter, require nchunk remain the same - if (firstflag) { - nchunk = n; - allocate(); - size_array_rows = nchunk; - } else if (n != nchunk) - error->all(FLERR,"Compute msd/chunk nchunk is not static"); + if (!firstflag && (oldnchunk != nchunk)) + error->all(FLERR, "Compute msd/chunk nchunk is not static"); // zero local per-chunk values @@ -171,19 +146,21 @@ void ComputeMSDChunk::compute_array() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - domain->unmap(x[i],image[i],unwrap); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + domain->unmap(x[i], image[i], unwrap); massproc[index] += massone; 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(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(massproc, masstotal, nchunk, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(&com[0][0], &comall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); for (int i = 0; i < nchunk; i++) { if (masstotal[i] > 0.0) { @@ -198,87 +175,32 @@ void ComputeMSDChunk::compute_array() if (firstflag) return; - double dx,dy,dz; + double dx, dy, dz; double **cominit = fix->astore; for (int i = 0; i < nchunk; i++) { dx = comall[i][0] - cominit[i][0]; dy = comall[i][1] - cominit[i][1]; dz = comall[i][2] - cominit[i][2]; - msd[i][0] = dx*dx; - msd[i][1] = dy*dy; - msd[i][2] = dz*dz; - msd[i][3] = dx*dx + dy*dy + dz*dz; + msd[i][0] = dx * dx; + msd[i][1] = dy * dy; + msd[i][2] = dz * dz; + msd[i][3] = dx * dx + dy * dy + dz * dz; } } -/* ---------------------------------------------------------------------- - 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 ComputeMSDChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputeMSDChunk::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 ComputeMSDChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeMSDChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeMSDChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- one-time allocate of per-chunk arrays ------------------------------------------------------------------------- */ void ComputeMSDChunk::allocate() { - memory->create(massproc,nchunk,"msd/chunk:massproc"); - memory->create(masstotal,nchunk,"msd/chunk:masstotal"); - memory->create(com,nchunk,3,"msd/chunk:com"); - memory->create(comall,nchunk,3,"msd/chunk:comall"); - memory->create(msd,nchunk,4,"msd/chunk:msd"); + ComputeChunk::allocate(); + memory->create(massproc, nchunk, "msd/chunk:massproc"); + memory->create(masstotal, nchunk, "msd/chunk:masstotal"); + memory->create(com, nchunk, 3, "msd/chunk:com"); + memory->create(comall, nchunk, 3, "msd/chunk:comall"); + memory->create(msd, nchunk, 4, "msd/chunk:msd"); array = msd; } @@ -288,8 +210,9 @@ void ComputeMSDChunk::allocate() double ComputeMSDChunk::memory_usage() { - double bytes = (bigint) nchunk * 2 * sizeof(double); - bytes += (double) nchunk * 2*3 * sizeof(double); + double bytes = ComputeChunk::memory_usage(); + bytes += (bigint) nchunk * 2 * sizeof(double); + bytes += (double) nchunk * 2 * 3 * sizeof(double); bytes += (double) nchunk * 4 * sizeof(double); return bytes; } diff --git a/src/compute_msd_chunk.h b/src/compute_msd_chunk.h index 7357f440a1..aba0c25fcb 100644 --- a/src/compute_msd_chunk.h +++ b/src/compute_msd_chunk.h @@ -20,11 +20,11 @@ ComputeStyle(msd/chunk,ComputeMSDChunk); #ifndef LMP_COMPUTE_MSD_CHUNK_H #define LMP_COMPUTE_MSD_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { -class ComputeMSDChunk : public Compute { +class ComputeMSDChunk : public ComputeChunk { public: ComputeMSDChunk(class LAMMPS *, int, char **); ~ComputeMSDChunk() override; @@ -32,27 +32,17 @@ class ComputeMSDChunk : public Compute { void setup() override; void compute_array() override; - void lock_enable() override; - void lock_disable() override; - int lock_length() override; - void lock(class Fix *, bigint, bigint) override; - void unlock(class Fix *) override; - double memory_usage() override; private: - int nchunk; - char *idchunk; - class ComputeChunkAtom *cchunk; char *id_fix; class FixStoreGlobal *fix; - int firstflag; double *massproc, *masstotal; double **com, **comall; double **msd; - void allocate(); + void allocate() override; }; } // namespace LAMMPS_NS #endif diff --git a/src/compute_omega_chunk.cpp b/src/compute_omega_chunk.cpp index 159c865188..3c345ab7a4 100644 --- a/src/compute_omega_chunk.cpp +++ b/src/compute_omega_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,16 +13,13 @@ #include "compute_omega_chunk.h" -#include #include "atom.h" -#include "update.h" -#include "modify.h" #include "compute_chunk_atom.h" #include "domain.h" -#include "math_extra.h" -#include "math_eigen.h" -#include "memory.h" #include "error.h" +#include "math_eigen.h" +#include "math_extra.h" +#include "memory.h" using namespace LAMMPS_NS; @@ -32,11 +28,11 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeOmegaChunk::ComputeOmegaChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - idchunk(nullptr),massproc(nullptr),masstotal(nullptr),com(nullptr),comall(nullptr), - inertia(nullptr),inertiaall(nullptr),angmom(nullptr),angmomall(nullptr),omega(nullptr) + ComputeChunk(lmp, narg, arg), massproc(nullptr), masstotal(nullptr), com(nullptr), + comall(nullptr), inertia(nullptr), inertiaall(nullptr), angmom(nullptr), angmomall(nullptr), + omega(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute omega/chunk command"); + if (narg != 4) error->all(FLERR, "Illegal compute omega/chunk command"); array_flag = 1; size_array_cols = 3; @@ -44,24 +40,14 @@ ComputeOmegaChunk::ComputeOmegaChunk(LAMMPS *lmp, int narg, char **arg) : size_array_rows_variable = 1; extarray = 0; - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); - ComputeOmegaChunk::init(); - - // chunk-based data - - nchunk = 1; - maxchunk = 0; - allocate(); + ComputeOmegaChunk::allocate(); } /* ---------------------------------------------------------------------- */ ComputeOmegaChunk::~ComputeOmegaChunk() { - delete [] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -75,33 +61,13 @@ ComputeOmegaChunk::~ComputeOmegaChunk() /* ---------------------------------------------------------------------- */ -void ComputeOmegaChunk::init() -{ - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for " - "compute omega/chunk"); - cchunk = dynamic_cast(modify->compute[icompute]); - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute omega/chunk does not use chunk/atom compute"); -} - -/* ---------------------------------------------------------------------- */ - void ComputeOmegaChunk::compute_array() { - int i,j,m,index; - double dx,dy,dz,massone; + int i, j, m, index; + double dx, dy, dz, 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(); + ComputeChunk::compute_array(); int *ichunk = cchunk->ichunk; if (nchunk > maxchunk) allocate(); @@ -129,19 +95,21 @@ void ComputeOmegaChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - domain->unmap(x[i],image[i],unwrap); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + domain->unmap(x[i], image[i], unwrap); massproc[index] += massone; 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(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(massproc, masstotal, 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) { @@ -155,24 +123,25 @@ void ComputeOmegaChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - domain->unmap(x[i],image[i],unwrap); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + domain->unmap(x[i], image[i], unwrap); dx = unwrap[0] - comall[index][0]; dy = unwrap[1] - comall[index][1]; dz = unwrap[2] - comall[index][2]; - inertia[index][0] += massone * (dy*dy + dz*dz); - inertia[index][1] += massone * (dx*dx + dz*dz); - inertia[index][2] += massone * (dx*dx + dy*dy); - inertia[index][3] -= massone * dx*dy; - inertia[index][4] -= massone * dy*dz; - inertia[index][5] -= massone * dx*dz; + inertia[index][0] += massone * (dy * dy + dz * dz); + inertia[index][1] += massone * (dx * dx + dz * dz); + inertia[index][2] += massone * (dx * dx + dy * dy); + inertia[index][3] -= massone * dx * dy; + inertia[index][4] -= massone * dy * dz; + inertia[index][5] -= massone * dx * dz; } - MPI_Allreduce(&inertia[0][0],&inertiaall[0][0],6*nchunk, - MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&inertia[0][0], &inertiaall[0][0], 6 * nchunk, MPI_DOUBLE, MPI_SUM, world); // compute angmom for each chunk @@ -180,37 +149,38 @@ void ComputeOmegaChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - domain->unmap(x[i],image[i],unwrap); + domain->unmap(x[i], image[i], unwrap); dx = unwrap[0] - comall[index][0]; dy = unwrap[1] - comall[index][1]; dz = unwrap[2] - comall[index][2]; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - angmom[index][0] += massone * (dy*v[i][2] - dz*v[i][1]); - angmom[index][1] += massone * (dz*v[i][0] - dx*v[i][2]); - angmom[index][2] += massone * (dx*v[i][1] - dy*v[i][0]); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + angmom[index][0] += massone * (dy * v[i][2] - dz * v[i][1]); + angmom[index][1] += massone * (dz * v[i][0] - dx * v[i][2]); + angmom[index][2] += massone * (dx * v[i][1] - dy * v[i][0]); } - MPI_Allreduce(&angmom[0][0],&angmomall[0][0],3*nchunk, - MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&angmom[0][0], &angmomall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); // compute omega for each chunk - double determinant,invdeterminant; - double idiag[3],ex[3],ey[3],ez[3],cross[3]; - double ione[3][3],inverse[3][3],evectors[3][3]; - double *iall,*mall; + double determinant, invdeterminant; + double idiag[3], ex[3], ey[3], ez[3], cross[3]; + double ione[3][3], inverse[3][3], evectors[3][3]; + double *iall, *mall; for (m = 0; m < nchunk; m++) { // determinant = triple product of rows of inertia matrix iall = &inertiaall[m][0]; - determinant = iall[0] * (iall[1]*iall[2] - iall[4]*iall[4]) + - iall[3] * (iall[4]*iall[5] - iall[3]*iall[2]) + - iall[5] * (iall[3]*iall[4] - iall[1]*iall[5]); + determinant = iall[0] * (iall[1] * iall[2] - iall[4] * iall[4]) + + iall[3] * (iall[4] * iall[5] - iall[3] * iall[2]) + + iall[5] * (iall[3] * iall[4] - iall[1] * iall[5]); ione[0][0] = iall[0]; ione[1][1] = iall[1]; @@ -223,39 +193,34 @@ void ComputeOmegaChunk::compute_array() // use L = Iw, inverting I to solve for w if (determinant > EPSILON) { - inverse[0][0] = ione[1][1]*ione[2][2] - ione[1][2]*ione[2][1]; - inverse[0][1] = -(ione[0][1]*ione[2][2] - ione[0][2]*ione[2][1]); - inverse[0][2] = ione[0][1]*ione[1][2] - ione[0][2]*ione[1][1]; + inverse[0][0] = ione[1][1] * ione[2][2] - ione[1][2] * ione[2][1]; + inverse[0][1] = -(ione[0][1] * ione[2][2] - ione[0][2] * ione[2][1]); + inverse[0][2] = ione[0][1] * ione[1][2] - ione[0][2] * ione[1][1]; - inverse[1][0] = -(ione[1][0]*ione[2][2] - ione[1][2]*ione[2][0]); - inverse[1][1] = ione[0][0]*ione[2][2] - ione[0][2]*ione[2][0]; - inverse[1][2] = -(ione[0][0]*ione[1][2] - ione[0][2]*ione[1][0]); + inverse[1][0] = -(ione[1][0] * ione[2][2] - ione[1][2] * ione[2][0]); + inverse[1][1] = ione[0][0] * ione[2][2] - ione[0][2] * ione[2][0]; + inverse[1][2] = -(ione[0][0] * ione[1][2] - ione[0][2] * ione[1][0]); - inverse[2][0] = ione[1][0]*ione[2][1] - ione[1][1]*ione[2][0]; - inverse[2][1] = -(ione[0][0]*ione[2][1] - ione[0][1]*ione[2][0]); - inverse[2][2] = ione[0][0]*ione[1][1] - ione[0][1]*ione[1][0]; + inverse[2][0] = ione[1][0] * ione[2][1] - ione[1][1] * ione[2][0]; + inverse[2][1] = -(ione[0][0] * ione[2][1] - ione[0][1] * ione[2][0]); + inverse[2][2] = ione[0][0] * ione[1][1] - ione[0][1] * ione[1][0]; - invdeterminant = 1.0/determinant; + invdeterminant = 1.0 / determinant; for (i = 0; i < 3; i++) - for (j = 0; j < 3; j++) - inverse[i][j] *= invdeterminant; + for (j = 0; j < 3; j++) inverse[i][j] *= invdeterminant; mall = &angmomall[m][0]; - omega[m][0] = inverse[0][0]*mall[0] + inverse[0][1]*mall[1] + - inverse[0][2]*mall[2]; - omega[m][1] = inverse[1][0]*mall[0] + inverse[1][1]*mall[1] + - inverse[1][2]*mall[2]; - omega[m][2] = inverse[2][0]*mall[0] + inverse[2][1]*mall[1] + - inverse[2][2]*mall[2]; + omega[m][0] = inverse[0][0] * mall[0] + inverse[0][1] * mall[1] + inverse[0][2] * mall[2]; + omega[m][1] = inverse[1][0] * mall[0] + inverse[1][1] * mall[1] + inverse[1][2] * mall[2]; + omega[m][2] = inverse[2][0] * mall[0] + inverse[2][1] * mall[1] + inverse[2][2] * mall[2]; - // handle each (nearly) singular I matrix - // due to 2-atom chunk or linear molecule - // use jacobi3() and angmom_to_omega() to calculate valid omega + // handle each (nearly) singular I matrix + // due to 2-atom chunk or linear molecule + // use jacobi3() and angmom_to_omega() to calculate valid omega } else { - int ierror = MathEigen::jacobi3(ione,idiag,evectors); - if (ierror) error->all(FLERR, - "Insufficient Jacobi rotations for omega/chunk"); + int ierror = MathEigen::jacobi3(ione, idiag, evectors); + if (ierror) error->all(FLERR, "Insufficient Jacobi rotations for omega/chunk"); ex[0] = evectors[0][0]; ex[1] = evectors[1][0]; @@ -270,88 +235,32 @@ void ComputeOmegaChunk::compute_array() // enforce 3 evectors as a right-handed coordinate system // flip 3rd vector if needed - MathExtra::cross3(ex,ey,cross); - if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez); + MathExtra::cross3(ex, ey, cross); + if (MathExtra::dot3(cross, ez) < 0.0) MathExtra::negate3(ez); // if any principal moment < scaled EPSILON, set to 0.0 double max; - max = MAX(idiag[0],idiag[1]); - max = MAX(max,idiag[2]); + max = MAX(idiag[0], idiag[1]); + max = MAX(max, idiag[2]); - if (idiag[0] < EPSILON*max) idiag[0] = 0.0; - if (idiag[1] < EPSILON*max) idiag[1] = 0.0; - if (idiag[2] < EPSILON*max) idiag[2] = 0.0; + if (idiag[0] < EPSILON * max) idiag[0] = 0.0; + if (idiag[1] < EPSILON * max) idiag[1] = 0.0; + if (idiag[2] < EPSILON * max) idiag[2] = 0.0; // calculate omega using diagonalized inertia matrix - MathExtra::angmom_to_omega(&angmomall[m][0],ex,ey,ez,idiag,&omega[m][0]); + MathExtra::angmom_to_omega(&angmomall[m][0], ex, ey, ez, idiag, &omega[m][0]); } } } - -/* ---------------------------------------------------------------------- - 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 ComputeOmegaChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputeOmegaChunk::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 ComputeOmegaChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeOmegaChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeOmegaChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ void ComputeOmegaChunk::allocate() { + ComputeChunk::allocate(); memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -362,15 +271,15 @@ void ComputeOmegaChunk::allocate() memory->destroy(angmomall); memory->destroy(omega); maxchunk = nchunk; - memory->create(massproc,maxchunk,"omega/chunk:massproc"); - memory->create(masstotal,maxchunk,"omega/chunk:masstotal"); - memory->create(com,maxchunk,3,"omega/chunk:com"); - memory->create(comall,maxchunk,3,"omega/chunk:comall"); - memory->create(inertia,maxchunk,6,"omega/chunk:inertia"); - memory->create(inertiaall,maxchunk,6,"omega/chunk:inertiaall"); - memory->create(angmom,maxchunk,3,"omega/chunk:angmom"); - memory->create(angmomall,maxchunk,3,"omega/chunk:angmomall"); - memory->create(omega,maxchunk,3,"omega/chunk:omega"); + memory->create(massproc, maxchunk, "omega/chunk:massproc"); + memory->create(masstotal, maxchunk, "omega/chunk:masstotal"); + memory->create(com, maxchunk, 3, "omega/chunk:com"); + memory->create(comall, maxchunk, 3, "omega/chunk:comall"); + memory->create(inertia, maxchunk, 6, "omega/chunk:inertia"); + memory->create(inertiaall, maxchunk, 6, "omega/chunk:inertiaall"); + memory->create(angmom, maxchunk, 3, "omega/chunk:angmom"); + memory->create(angmomall, maxchunk, 3, "omega/chunk:angmomall"); + memory->create(omega, maxchunk, 3, "omega/chunk:omega"); array = omega; } @@ -380,10 +289,11 @@ void ComputeOmegaChunk::allocate() double ComputeOmegaChunk::memory_usage() { - double bytes = (bigint) maxchunk * 2 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); - bytes += (double) maxchunk * 2*6 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); + double bytes = ComputeChunk::memory_usage(); + bytes += (bigint) maxchunk * 2 * sizeof(double); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); + bytes += (double) maxchunk * 2 * 6 * sizeof(double); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); bytes += (double) maxchunk * 3 * sizeof(double); return bytes; } diff --git a/src/compute_omega_chunk.h b/src/compute_omega_chunk.h index cbd158c525..6265be2010 100644 --- a/src/compute_omega_chunk.h +++ b/src/compute_omega_chunk.h @@ -20,40 +20,28 @@ ComputeStyle(omega/chunk,ComputeOmegaChunk); #ifndef LMP_COMPUTE_OMEGA_CHUNK_H #define LMP_COMPUTE_OMEGA_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { -class ComputeOmegaChunk : public Compute { +class ComputeOmegaChunk : public ComputeChunk { public: ComputeOmegaChunk(class LAMMPS *, int, char **); ~ComputeOmegaChunk() override; - void init() override; - void compute_array() override; - void lock_enable() override; - void lock_disable() override; - int lock_length() override; - void lock(class Fix *, bigint, bigint) override; - void unlock(class Fix *) override; + void compute_array() override; double memory_usage() override; private: - int nchunk, maxchunk; - char *idchunk; - class ComputeChunkAtom *cchunk; - double *massproc, *masstotal; double **com, **comall; double **inertia, **inertiaall; double **angmom, **angmomall; double **omega; - void allocate(); + void allocate() override; }; - } // namespace LAMMPS_NS - #endif #endif diff --git a/src/compute_property_chunk.cpp b/src/compute_property_chunk.cpp index 1fa8cf422b..190a75d1bd 100644 --- a/src/compute_property_chunk.cpp +++ b/src/compute_property_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,27 +13,21 @@ #include "compute_property_chunk.h" -#include #include "atom.h" -#include "update.h" -#include "modify.h" #include "compute_chunk_atom.h" -#include "memory.h" #include "error.h" +#include "memory.h" + +#include using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputePropertyChunk::ComputePropertyChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - idchunk(nullptr), count_one(nullptr), count_all(nullptr) + ComputeChunk(lmp, narg, arg), count_one(nullptr), count_all(nullptr) { - if (narg < 5) error->all(FLERR,"Illegal compute property/chunk command"); - - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); + if (narg < 5) utils::missing_cmd_args(FLERR, "compute property/chunk", error); ComputePropertyChunk::init(); @@ -46,40 +39,32 @@ ComputePropertyChunk::ComputePropertyChunk(LAMMPS *lmp, int narg, char **arg) : int i; for (int iarg = 4; iarg < narg; iarg++) { - i = iarg-4; + i = iarg - 4; - if (strcmp(arg[iarg],"count") == 0) { + if (strcmp(arg[iarg], "count") == 0) { pack_choice[i] = &ComputePropertyChunk::pack_count; countflag = 1; - } else if (strcmp(arg[iarg],"id") == 0) { + } else if (strcmp(arg[iarg], "id") == 0) { if (!cchunk->compress) - error->all(FLERR,"Compute chunk/atom stores no IDs for " - "compute property/chunk"); + error->all(FLERR, "Compute chunk/atom stores no IDs for compute property/chunk"); pack_choice[i] = &ComputePropertyChunk::pack_id; - } else if (strcmp(arg[iarg],"coord1") == 0) { + } else if (strcmp(arg[iarg], "coord1") == 0) { if (cchunk->ncoord < 1) - error->all(FLERR,"Compute chunk/atom stores no coord1 for " - "compute property/chunk"); + error->all(FLERR, "Compute chunk/atom stores no coord1 for compute property/chunk"); pack_choice[i] = &ComputePropertyChunk::pack_coord1; - } else if (strcmp(arg[iarg],"coord2") == 0) { + } else if (strcmp(arg[iarg], "coord2") == 0) { if (cchunk->ncoord < 2) - error->all(FLERR,"Compute chunk/atom stores no coord2 for " - "compute property/chunk"); + error->all(FLERR, "Compute chunk/atom stores no coord2 for compute property/chunk"); pack_choice[i] = &ComputePropertyChunk::pack_coord2; - } else if (strcmp(arg[iarg],"coord3") == 0) { + } else if (strcmp(arg[iarg], "coord3") == 0) { if (cchunk->ncoord < 3) - error->all(FLERR,"Compute chunk/atom stores no coord3 for " - "compute property/chunk"); + error->all(FLERR, "Compute chunk/atom stores no coord3 for compute property/chunk"); pack_choice[i] = &ComputePropertyChunk::pack_coord3; - } else error->all(FLERR, - "Invalid keyword in compute property/chunk command"); + } else + error->all(FLERR, "Unkown keyword {} in compute property/chunk command", arg[iarg]); } - // initialization - - nchunk = 1; - maxchunk = 0; - allocate(); + ComputePropertyChunk::allocate(); if (nvalues == 1) { vector_flag = 1; @@ -99,8 +84,7 @@ ComputePropertyChunk::ComputePropertyChunk(LAMMPS *lmp, int narg, char **arg) : ComputePropertyChunk::~ComputePropertyChunk() { - delete [] idchunk; - delete [] pack_choice; + delete[] pack_choice; memory->destroy(vector); memory->destroy(array); memory->destroy(count_one); @@ -109,36 +93,11 @@ ComputePropertyChunk::~ComputePropertyChunk() /* ---------------------------------------------------------------------- */ -void ComputePropertyChunk::init() -{ - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for " - "compute property/chunk"); - cchunk = dynamic_cast(modify->compute[icompute]); - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute property/chunk does not use chunk/atom compute"); -} - -/* ---------------------------------------------------------------------- */ - void ComputePropertyChunk::compute_vector() { - invoked_vector = update->ntimestep; + ComputeChunk::compute_vector(); - // compute chunk/atom assigns atoms to chunk IDs - // if need count, extract ichunk index vector from compute - // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms - - nchunk = cchunk->setup_chunks(); - if (nchunk > maxchunk) allocate(); - if (nvalues == 1) size_vector = nchunk; - else size_array_rows = nchunk; - - if (countflag) { - cchunk->compute_ichunk(); - ichunk = cchunk->ichunk; - } + if (countflag) ichunk = cchunk->ichunk; // fill vector @@ -150,101 +109,34 @@ void ComputePropertyChunk::compute_vector() void ComputePropertyChunk::compute_array() { - invoked_array = update->ntimestep; + ComputeChunk::compute_array(); - // compute chunk/atom assigns atoms to chunk IDs - // if need count, extract ichunk index vector from compute - // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms - - nchunk = cchunk->setup_chunks(); - if (nchunk > maxchunk) allocate(); - if (nvalues == 1) size_vector = nchunk; - else size_array_rows = nchunk; - - if (countflag) { - cchunk->compute_ichunk(); - ichunk = cchunk->ichunk; - } + if (countflag) ichunk = cchunk->ichunk; // fill array if (array) buf = &array[0][0]; - for (int n = 0; n < nvalues; n++) - (this->*pack_choice[n])(n); + for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n); } - -/* ---------------------------------------------------------------------- - 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 ComputePropertyChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputePropertyChunk::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 ComputePropertyChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputePropertyChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputePropertyChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ void ComputePropertyChunk::allocate() { + ComputeChunk::allocate(); memory->destroy(vector); memory->destroy(array); memory->destroy(count_one); memory->destroy(count_all); maxchunk = nchunk; - if (nvalues == 1) memory->create(vector,maxchunk,"property/chunk:vector"); - else memory->create(array,maxchunk,nvalues,"property/chunk:array"); + if (nvalues == 1) + memory->create(vector, maxchunk, "property/chunk:vector"); + else + memory->create(array, maxchunk, nvalues, "property/chunk:array"); if (countflag) { - memory->create(count_one,maxchunk,"property/chunk:count_one"); - memory->create(count_all,maxchunk,"property/chunk:count_all"); + memory->create(count_one, maxchunk, "property/chunk:count_one"); + memory->create(count_all, maxchunk, "property/chunk:count_all"); } } @@ -254,7 +146,8 @@ void ComputePropertyChunk::allocate() double ComputePropertyChunk::memory_usage() { - double bytes = (bigint) nchunk * nvalues * sizeof(double); + double bytes = ComputeChunk::memory_usage(); + bytes += (bigint) nchunk * nvalues * sizeof(double); if (countflag) bytes += (double) nchunk * 2 * sizeof(int); return bytes; } @@ -277,13 +170,13 @@ void ComputePropertyChunk::pack_count(int n) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; count_one[index]++; } } - MPI_Allreduce(count_one,count_all,nchunk,MPI_INT,MPI_SUM,world); + MPI_Allreduce(count_one, count_all, nchunk, MPI_INT, MPI_SUM, world); for (int m = 0; m < nchunk; m++) { buf[n] = count_all[m]; diff --git a/src/compute_property_chunk.h b/src/compute_property_chunk.h index 7339d6479a..3bc130e443 100644 --- a/src/compute_property_chunk.h +++ b/src/compute_property_chunk.h @@ -20,37 +20,27 @@ ComputeStyle(property/chunk,ComputePropertyChunk); #ifndef LMP_COMPUTE_CHUNK_MOLECULE_H #define LMP_COMPUTE_CHUNK_MOLECULE_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { -class ComputePropertyChunk : public Compute { +class ComputePropertyChunk : public ComputeChunk { public: ComputePropertyChunk(class LAMMPS *, int, char **); ~ComputePropertyChunk() override; - void init() override; + void compute_vector() override; void compute_array() override; - void lock_enable() override; - void lock_disable() override; - int lock_length() override; - void lock(class Fix *, bigint, bigint) override; - void unlock(class Fix *) override; - double memory_usage() override; private: - int nchunk, maxchunk; - char *idchunk; - class ComputeChunkAtom *cchunk; int *ichunk; - int nvalues, countflag; double *buf; int *count_one, *count_all; - void allocate(); + void allocate() override; typedef void (ComputePropertyChunk::*FnPtrPack)(int); FnPtrPack *pack_choice; // ptrs to pack functions @@ -61,8 +51,6 @@ class ComputePropertyChunk : public Compute { void pack_coord2(int); void pack_coord3(int); }; - } // namespace LAMMPS_NS - #endif #endif diff --git a/src/compute_reduce_chunk.cpp b/src/compute_reduce_chunk.cpp index 992ab91bea..51781eac7b 100644 --- a/src/compute_reduce_chunk.cpp +++ b/src/compute_reduce_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -30,29 +29,30 @@ using namespace LAMMPS_NS; -enum{ SUM, MINN, MAXX }; +enum { SUM, MINN, MAXX }; #define BIG 1.0e20 /* ---------------------------------------------------------------------- */ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), idchunk(nullptr), vlocal(nullptr), vglobal(nullptr), - alocal(nullptr), aglobal(nullptr), varatom(nullptr), cchunk(nullptr), ichunk(nullptr) + ComputeChunk(lmp, narg, arg), vlocal(nullptr), vglobal(nullptr), alocal(nullptr), + aglobal(nullptr), varatom(nullptr), ichunk(nullptr) { - if (narg < 6) utils::missing_cmd_args(FLERR,"compute reduce/chunk", error); + if (narg < 6) utils::missing_cmd_args(FLERR, "compute reduce/chunk", error); - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); - init_chunk(); + ComputeChunk::init(); // mode - if (strcmp(arg[4],"sum") == 0) mode = SUM; - else if (strcmp(arg[4],"min") == 0) mode = MINN; - else if (strcmp(arg[4],"max") == 0) mode = MAXX; - else error->all(FLERR,"Unknown compute reduce/chunk mode: {}", arg[4]); + if (strcmp(arg[4], "sum") == 0) + mode = SUM; + else if (strcmp(arg[4], "min") == 0) + mode = MINN; + else if (strcmp(arg[4], "max") == 0) + mode = MAXX; + else + error->all(FLERR, "Unknown compute reduce/chunk mode: {}", arg[4]); int iarg = 5; @@ -60,7 +60,7 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) : int expand = 0; char **earg; - int nargnew = utils::expand_args(FLERR,narg-iarg,&arg[iarg],1,earg,lmp); + int nargnew = utils::expand_args(FLERR, narg - iarg, &arg[iarg], 1, earg, lmp); if (earg != &arg[iarg]) expand = 1; arg = earg; @@ -78,7 +78,7 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) : val.val.c = nullptr; if ((val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE) || (argi.get_dim() > 1)) - error->all(FLERR,"Illegal compute reduce/chunk argument: {}", arg[iarg]); + error->all(FLERR, "Illegal compute reduce/chunk argument: {}", arg[iarg]); values.push_back(val); } @@ -86,7 +86,7 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) : // if wildcard expansion occurred, free earg memory from expand_args() if (expand) { - for (int i = 0; i < nargnew; i++) delete [] earg[i]; + for (int i = 0; i < nargnew; i++) delete[] earg[i]; memory->sfree(earg); } @@ -96,15 +96,15 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) : if (val.which == ArgInfo::COMPUTE) { val.val.c = modify->get_compute_by_id(val.id); if (!val.val.c) - error->all(FLERR,"Compute ID {} for compute reduce/chunk does not exist", val.id); + error->all(FLERR, "Compute ID {} for compute reduce/chunk does not exist", val.id); if (!val.val.c->peratom_flag) - error->all(FLERR,"Compute reduce/chunk compute {} does not calculate per-atom values", + error->all(FLERR, "Compute reduce/chunk compute {} does not calculate per-atom values", val.id); if ((val.argindex == 0) && (val.val.c->size_peratom_cols != 0)) - error->all(FLERR,"Compute reduce/chunk compute {} does not calculate a per-atom vector", + error->all(FLERR, "Compute reduce/chunk compute {} does not calculate a per-atom vector", val.id); if (val.argindex && (val.val.c->size_peratom_cols == 0)) - error->all(FLERR,"Compute reduce/chunk compute {} does not calculate a per-atom array", + error->all(FLERR, "Compute reduce/chunk compute {} does not calculate a per-atom array", val.id); if (val.argindex && (val.argindex > val.val.c->size_peratom_cols)) error->all(FLERR, "Compute reduce/chunk compute array {} is accessed out-of-range", val.id); @@ -112,22 +112,24 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) : } else if (val.which == ArgInfo::FIX) { val.val.f = modify->get_fix_by_id(val.id); if (!val.val.f) - error->all(FLERR,"Fix ID {} for compute reduce/chunk does not exist", val.id); + error->all(FLERR, "Fix ID {} for compute reduce/chunk does not exist", val.id); if (!val.val.f->peratom_flag) - error->all(FLERR,"Compute reduce/chunk fix {} does not calculate per-atom values", val.id); + error->all(FLERR, "Compute reduce/chunk fix {} does not calculate per-atom values", val.id); if ((val.argindex == 0) && (val.val.f->size_peratom_cols != 0)) - error->all(FLERR,"Compute reduce/chunk fix {} does not calculate a per-atom vector", val.id); + error->all(FLERR, "Compute reduce/chunk fix {} does not calculate a per-atom vector", + val.id); if (val.argindex && (val.val.f->size_peratom_cols == 0)) - error->all(FLERR,"Compute reduce/chunk fix {} does not calculate a per-atom array", val.id); + error->all(FLERR, "Compute reduce/chunk fix {} does not calculate a per-atom array", + val.id); if (val.argindex && (val.argindex > val.val.f->size_peratom_cols)) - error->all(FLERR,"Compute reduce/chunk fix {} array is accessed out-of-range", val.id); + error->all(FLERR, "Compute reduce/chunk fix {} array is accessed out-of-range", val.id); } else if (val.which == ArgInfo::VARIABLE) { val.val.v = input->variable->find(val.id.c_str()); if (val.val.v < 0) - error->all(FLERR,"Variable name {} for compute reduce/chunk does not exist", val.id); + error->all(FLERR, "Variable name {} for compute reduce/chunk does not exist", val.id); if (input->variable->atomstyle(val.val.v) == 0) - error->all(FLERR,"Compute reduce/chunk variable is not atom-style variable"); + error->all(FLERR, "Compute reduce/chunk variable is not atom-style variable"); } } @@ -146,11 +148,13 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) : // setup - if (mode == SUM) initvalue = 0.0; - else if (mode == MINN) initvalue = BIG; - else if (mode == MAXX) initvalue = -BIG; + if (mode == SUM) + initvalue = 0.0; + else if (mode == MINN) + initvalue = BIG; + else if (mode == MAXX) + initvalue = -BIG; - maxchunk = 0; vlocal = vglobal = nullptr; alocal = aglobal = nullptr; @@ -162,8 +166,6 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) : ComputeReduceChunk::~ComputeReduceChunk() { - delete[] idchunk; - memory->destroy(vlocal); memory->destroy(vglobal); memory->destroy(alocal); @@ -176,7 +178,7 @@ ComputeReduceChunk::~ComputeReduceChunk() void ComputeReduceChunk::init() { - init_chunk(); + ComputeChunk::init(); // set indices of all computes,fixes,variables @@ -184,109 +186,84 @@ void ComputeReduceChunk::init() if (val.which == ArgInfo::COMPUTE) { val.val.c = modify->get_compute_by_id(val.id); if (!val.val.c) - error->all(FLERR,"Compute ID {} for compute reduce/chunk does not exist", val.id); + error->all(FLERR, "Compute ID {} for compute reduce/chunk does not exist", val.id); } else if (val.which == ArgInfo::FIX) { val.val.f = modify->get_fix_by_id(val.id); if (!val.val.f) - error->all(FLERR,"Fix ID {} for compute reduce/chunk does not exist", val.id); + error->all(FLERR, "Fix ID {} for compute reduce/chunk does not exist", val.id); } else if (val.which == ArgInfo::VARIABLE) { val.val.v = input->variable->find(val.id.c_str()); if (val.val.v < 0) - error->all(FLERR,"Variable name {} for compute reduce/chunk does not exist", val.id); + error->all(FLERR, "Variable name {} for compute reduce/chunk does not exist", val.id); } } } /* ---------------------------------------------------------------------- */ -void ComputeReduceChunk::init_chunk() -{ - cchunk = dynamic_cast(modify->get_compute_by_id(idchunk)); - if (!cchunk) - error->all(FLERR,"Compute chunk/atom {} does not exist or is incorrect style for " - "compute reduce/chunk", idchunk); -} - -/* ---------------------------------------------------------------------- */ - void ComputeReduceChunk::compute_vector() { - invoked_vector = 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(); + ComputeChunk::compute_vector(); ichunk = cchunk->ichunk; if (!nchunk) return; - size_vector = nchunk; - if (nchunk > maxchunk) { memory->destroy(vlocal); memory->destroy(vglobal); maxchunk = nchunk; - memory->create(vlocal,maxchunk,"reduce/chunk:vlocal"); - memory->create(vglobal,maxchunk,"reduce/chunk:vglobal"); + memory->create(vlocal, maxchunk, "reduce/chunk:vlocal"); + memory->create(vglobal, maxchunk, "reduce/chunk:vglobal"); vector = vglobal; } // perform local reduction of single peratom value - compute_one(0,vlocal,1); + compute_one(0, vlocal, 1); // reduce the per-chunk values across all procs if (mode == SUM) - MPI_Allreduce(vlocal,vglobal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(vlocal, vglobal, nchunk, MPI_DOUBLE, MPI_SUM, world); else if (mode == MINN) - MPI_Allreduce(vlocal,vglobal,nchunk,MPI_DOUBLE,MPI_MIN,world); + MPI_Allreduce(vlocal, vglobal, nchunk, MPI_DOUBLE, MPI_MIN, world); else if (mode == MAXX) - MPI_Allreduce(vlocal,vglobal,nchunk,MPI_DOUBLE,MPI_MAX,world); + MPI_Allreduce(vlocal, vglobal, nchunk, MPI_DOUBLE, MPI_MAX, world); } /* ---------------------------------------------------------------------- */ void ComputeReduceChunk::compute_array() { - 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(); + ComputeChunk::compute_array(); ichunk = cchunk->ichunk; if (!nchunk) return; - size_array_rows = nchunk; - if (nchunk > maxchunk) { memory->destroy(alocal); memory->destroy(aglobal); maxchunk = nchunk; - memory->create(alocal,maxchunk,values.size(),"reduce/chunk:alocal"); - memory->create(aglobal,maxchunk,values.size(),"reduce/chunk:aglobal"); + memory->create(alocal, maxchunk, values.size(), "reduce/chunk:alocal"); + memory->create(aglobal, maxchunk, values.size(), "reduce/chunk:aglobal"); array = aglobal; } // perform local reduction of all peratom values - for (std::size_t m = 0; m < values.size(); m++) compute_one(m,&alocal[0][m],values.size()); + for (std::size_t m = 0; m < values.size(); m++) compute_one(m, &alocal[0][m], values.size()); // reduce the per-chunk values across all procs if (mode == SUM) - MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*values.size(),MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&alocal[0][0], &aglobal[0][0], nchunk * values.size(), MPI_DOUBLE, MPI_SUM, + world); else if (mode == MINN) - MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*values.size(),MPI_DOUBLE,MPI_MIN,world); + MPI_Allreduce(&alocal[0][0], &aglobal[0][0], nchunk * values.size(), MPI_DOUBLE, MPI_MIN, + world); else if (mode == MAXX) - MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*values.size(),MPI_DOUBLE,MPI_MAX,world); + MPI_Allreduce(&alocal[0][0], &aglobal[0][0], nchunk * values.size(), MPI_DOUBLE, MPI_MAX, + world); } /* ---------------------------------------------------------------------- */ @@ -295,7 +272,7 @@ void ComputeReduceChunk::compute_one(int m, double *vchunk, int nstride) { // initialize per-chunk values in accumulation vector - for (std::size_t i = 0; i < values.size()*nchunk; i += nstride) vchunk[i] = initvalue; + for (std::size_t i = 0; i < values.size() * nchunk; i += nstride) vchunk[i] = initvalue; // loop over my atoms // use peratom input and chunk ID of each atom to update vector @@ -321,61 +298,61 @@ void ComputeReduceChunk::compute_one(int m, double *vchunk, int nstride) double *vcompute = val.val.c->vector_atom; for (int i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - combine(vchunk[index*nstride],vcompute[i]); + combine(vchunk[index * nstride], vcompute[i]); } } else { double **acompute = val.val.c->array_atom; int argindexm1 = val.argindex - 1; for (int i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - combine(vchunk[index*nstride],acompute[i][argindexm1]); + combine(vchunk[index * nstride], acompute[i][argindexm1]); } } - // access fix fields, check if fix frequency is a match + // access fix fields, check if fix frequency is a match } else if (val.which == ArgInfo::FIX) { if (update->ntimestep % val.val.f->peratom_freq) - error->all(FLERR,"Fix used in compute reduce/chunk not computed at compatible time"); + error->all(FLERR, "Fix used in compute reduce/chunk not computed at compatible time"); if (val.argindex == 0) { double *vfix = val.val.f->vector_atom; for (int i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - combine(vchunk[index*nstride],vfix[i]); + combine(vchunk[index * nstride], vfix[i]); } } else { double **afix = val.val.f->array_atom; int argindexm1 = val.argindex - 1; for (int i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - combine(vchunk[index*nstride],afix[i][argindexm1]); + combine(vchunk[index * nstride], afix[i][argindexm1]); } } - // evaluate atom-style variable + // evaluate atom-style variable } else if (val.which == ArgInfo::VARIABLE) { if (atom->nmax > maxatom) { memory->destroy(varatom); maxatom = atom->nmax; - memory->create(varatom,maxatom,"reduce/chunk:varatom"); + memory->create(varatom, maxatom, "reduce/chunk:varatom"); } - input->variable->compute_atom(val.val.v,igroup,varatom,1,0); + input->variable->compute_atom(val.val.v, igroup, varatom, 1, 0); for (int i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - combine(vchunk[index*nstride],varatom[i]); + combine(vchunk[index * nstride], varatom[i]); } } } @@ -386,7 +363,8 @@ void ComputeReduceChunk::compute_one(int m, double *vchunk, int nstride) void ComputeReduceChunk::combine(double &one, double two) { - if (mode == SUM) one += two; + if (mode == SUM) + one += two; else if (mode == MINN) { if (two < one) one = two; } else if (mode == MAXX) { @@ -394,67 +372,16 @@ void ComputeReduceChunk::combine(double &one, double two) } } -/* ---------------------------------------------------------------------- - 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 ComputeReduceChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, if it still exists -------------------------------------------------------------------------- */ - -void ComputeReduceChunk::lock_disable() -{ - cchunk = dynamic_cast(modify->get_compute_by_id(idchunk)); - if (cchunk) cchunk->lockcount--; -} - -/* ---------------------------------------------------------------------- - calculate and return # of chunks = length of vector/array -------------------------------------------------------------------------- */ - -int ComputeReduceChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeReduceChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeReduceChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- memory usage of local data ------------------------------------------------------------------------- */ double ComputeReduceChunk::memory_usage() { - double bytes = (bigint) maxatom * sizeof(double); - if (values.size() == 1) bytes += (double) maxchunk * 2 * sizeof(double); - else bytes += (double) maxchunk * values.size() * 2 * sizeof(double); + double bytes = (double) maxatom * sizeof(double) + ComputeChunk::memory_usage(); + if (values.size() == 1) + bytes += (double) maxchunk * 2 * sizeof(double); + else + bytes += (double) maxchunk * values.size() * 2 * sizeof(double); return bytes; } diff --git a/src/compute_reduce_chunk.h b/src/compute_reduce_chunk.h index 505e2b4a21..4055956d2d 100644 --- a/src/compute_reduce_chunk.h +++ b/src/compute_reduce_chunk.h @@ -20,11 +20,11 @@ ComputeStyle(reduce/chunk,ComputeReduceChunk); #ifndef LMP_COMPUTE_REDUCE_CHUNK_H #define LMP_COMPUTE_REDUCE_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { -class ComputeReduceChunk : public Compute { +class ComputeReduceChunk : public ComputeChunk { public: ComputeReduceChunk(class LAMMPS *, int, char **); ~ComputeReduceChunk() override; @@ -32,12 +32,6 @@ class ComputeReduceChunk : public Compute { void compute_vector() override; void compute_array() override; - void lock_enable() override; - void lock_disable() override; - int lock_length() override; - void lock(class Fix *, bigint, bigint) override; - void unlock(class Fix *) override; - double memory_usage() override; private: @@ -53,19 +47,14 @@ class ComputeReduceChunk : public Compute { }; std::vector values; - char *idchunk; - - int mode, nchunk; - int maxchunk, maxatom; + int mode, maxatom; double initvalue; double *vlocal, *vglobal; double **alocal, **aglobal; double *varatom; - class ComputeChunkAtom *cchunk; int *ichunk; - void init_chunk(); void compute_one(int, double *, int); void combine(double &, double); }; diff --git a/src/compute_temp_chunk.cpp b/src/compute_temp_chunk.cpp index e27dcb49f2..72c2471b87 100644 --- a/src/compute_temp_chunk.cpp +++ b/src/compute_temp_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -27,42 +26,40 @@ using namespace LAMMPS_NS; -enum{TEMP,KECOM,INTERNAL}; +enum { TEMP, KECOM, INTERNAL }; /* ---------------------------------------------------------------------- */ ComputeTempChunk::ComputeTempChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - which(nullptr), idchunk(nullptr), id_bias(nullptr), sum(nullptr), sumall(nullptr), count(nullptr), - countall(nullptr), massproc(nullptr), masstotal(nullptr), vcm(nullptr), vcmall(nullptr) + ComputeChunk(lmp, narg, arg), which(nullptr), id_bias(nullptr), sum(nullptr), sumall(nullptr), + count(nullptr), countall(nullptr), massproc(nullptr), masstotal(nullptr), vcm(nullptr), + vcmall(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute temp/chunk command"); - scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); - biasflag = 0; ComputeTempChunk::init(); // optional per-chunk values - nvalues = narg-4; + nvalues = narg - 4; which = new int[nvalues]; nvalues = 0; int iarg = 4; while (iarg < narg) { - if (strcmp(arg[iarg],"temp") == 0) which[nvalues] = TEMP; - else if (strcmp(arg[iarg],"kecom") == 0) which[nvalues] = KECOM; - else if (strcmp(arg[iarg],"internal") == 0) which[nvalues] = INTERNAL; - else break; + if (strcmp(arg[iarg], "temp") == 0) + which[nvalues] = TEMP; + else if (strcmp(arg[iarg], "kecom") == 0) + which[nvalues] = KECOM; + else if (strcmp(arg[iarg], "internal") == 0) + which[nvalues] = INTERNAL; + else + break; iarg++; nvalues++; } @@ -76,56 +73,49 @@ ComputeTempChunk::ComputeTempChunk(LAMMPS *lmp, int narg, char **arg) : cdof = 0.0; while (iarg < narg) { - if (strcmp(arg[iarg],"com") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/chunk command"); - comflag = utils::logical(FLERR,arg[iarg+1],false,lmp); + if (strcmp(arg[iarg], "com") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute temp/chunk command"); + comflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bias") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/chunk command"); + } else if (strcmp(arg[iarg], "bias") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute temp/chunk command"); biasflag = 1; - id_bias = utils::strdup(arg[iarg+1]); + id_bias = utils::strdup(arg[iarg + 1]); iarg += 2; - } else if (strcmp(arg[iarg],"adof") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/chunk command"); - adof = utils::numeric(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "adof") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute temp/chunk command"); + adof = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"cdof") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/chunk command"); - cdof = utils::numeric(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "cdof") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute temp/chunk command"); + cdof = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else error->all(FLERR,"Illegal compute temp/chunk command"); + } else + error->all(FLERR, "Illegal compute temp/chunk command"); } // error check on bias compute if (biasflag) { - int i = modify->find_compute(id_bias); - if (i < 0) - error->all(FLERR,"Could not find compute ID for temperature bias"); - tbias = modify->compute[i]; - if (tbias->tempflag == 0) - error->all(FLERR,"Bias compute does not calculate temperature"); - if (tbias->tempbias == 0) - error->all(FLERR,"Bias compute does not calculate a velocity bias"); + tbias = modify->get_compute_by_id(id_bias); + if (!tbias) error->all(FLERR, "Could not find compute {} for temperature bias", id_bias); + + if (tbias->tempflag == 0) error->all(FLERR, "Bias compute does not calculate temperature"); + if (tbias->tempbias == 0) error->all(FLERR, "Bias compute does not calculate a velocity bias"); } // this compute only calculates a bias, if comflag is set // won't be two biases since comflag and biasflag cannot both be set if (comflag && biasflag) - error->all(FLERR,"Cannot use both com and bias with compute temp/chunk"); + error->all(FLERR, "Cannot use both com and bias with compute temp/chunk"); if (comflag) tempbias = 1; // vector data vector = new double[size_vector]; - // chunk-based data - - nchunk = 1; - maxchunk = 0; - - if (nvalues) { + if (nvalues) { array_flag = 1; size_array_cols = nvalues; size_array_rows = 0; @@ -133,7 +123,7 @@ ComputeTempChunk::ComputeTempChunk(LAMMPS *lmp, int narg, char **arg) : extarray = 0; } - allocate(); + ComputeTempChunk::allocate(); comstep = -1; } @@ -141,10 +131,9 @@ ComputeTempChunk::ComputeTempChunk(LAMMPS *lmp, int narg, char **arg) : ComputeTempChunk::~ComputeTempChunk() { - delete [] idchunk; - delete [] which; - delete [] id_bias; - delete [] vector; + delete[] which; + delete[] id_bias; + delete[] vector; memory->destroy(sum); memory->destroy(sumall); memory->destroy(count); @@ -160,19 +149,11 @@ ComputeTempChunk::~ComputeTempChunk() void ComputeTempChunk::init() { - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for " - "compute temp/chunk"); - cchunk = dynamic_cast(modify->compute[icompute]); - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute temp/chunk does not use chunk/atom compute"); + ComputeChunk::init(); if (biasflag) { - int i = modify->find_compute(id_bias); - if (i < 0) - error->all(FLERR,"Could not find compute ID for temperature bias"); - tbias = modify->compute[i]; + tbias = modify->get_compute_by_id(id_bias); + if (!tbias) error->all(FLERR, "Could not find compute ID {} for temperature bias", id_bias); } } @@ -180,7 +161,7 @@ void ComputeTempChunk::init() double ComputeTempChunk::compute_scalar() { - int i,index; + int i, index; invoked_scalar = update->ntimestep; @@ -224,47 +205,45 @@ double ComputeTempChunk::compute_scalar() if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - rmass[i]; + t += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]) * rmass[i]; mycount++; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - mass[type[i]]; + t += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]) * mass[type[i]]; mycount++; - } + } } } else { - double vx,vy,vz; + double vx, vy, vz; if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; - t += (vx*vx + vy*vy + vz*vz) * rmass[i]; + t += (vx * vx + vy * vy + vz * vz) * rmass[i]; mycount++; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; - t += (vx*vx + vy*vy + vz*vz) * mass[type[i]]; + t += (vx * vx + vy * vy + vz * vz) * mass[type[i]]; mycount++; - } + } } } @@ -274,16 +253,15 @@ double ComputeTempChunk::compute_scalar() // final temperature - MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&t, &scalar, 1, MPI_DOUBLE, MPI_SUM, world); double rcount = mycount; double allcount; - MPI_Allreduce(&rcount,&allcount,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&rcount, &allcount, 1, MPI_DOUBLE, MPI_SUM, world); - double dof = nchunk*cdof + adof*allcount; + double dof = nchunk * cdof + adof * allcount; double tfactor = 0.0; if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); - if (dof < 0.0 && allcount > 0.0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); + if (dof < 0.0 && allcount > 0.0) error->all(FLERR, "Temperature compute degrees of freedom < 0"); scalar *= tfactor; return scalar; } @@ -292,22 +270,11 @@ double ComputeTempChunk::compute_scalar() void ComputeTempChunk::compute_vector() { - int i,index; + int i, index; - invoked_vector = update->ntimestep; - - // calculate chunk assignments, - // since only atoms in chunks contribute to global temperature - // 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(); + ComputeChunk::compute_vector(); int *ichunk = cchunk->ichunk; - if (nchunk > maxchunk) allocate(); - // remove velocity bias if (biasflag) { @@ -329,40 +296,44 @@ void ComputeTempChunk::compute_vector() int *mask = atom->mask; int nlocal = atom->nlocal; - double massone,t[6]; + double massone, t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; if (!comflag) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + t[0] += massone * v[i][0] * v[i][0]; + t[1] += massone * v[i][1] * v[i][1]; + t[2] += massone * v[i][2] * v[i][2]; + t[3] += massone * v[i][0] * v[i][1]; + t[4] += massone * v[i][0] * v[i][2]; + t[5] += massone * v[i][1] * v[i][2]; } } else { - double vx,vy,vz; + double vx, vy, vz; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; - t[0] += massone * vx*vx; - t[1] += massone * vy*vy; - t[2] += massone * vz*vz; - t[3] += massone * vx*vy; - t[4] += massone * vx*vz; - t[5] += massone * vy*vz; + t[0] += massone * vx * vx; + t[1] += massone * vy * vy; + t[2] += massone * vz * vz; + t[3] += massone * vx * vy; + t[4] += massone * vx * vz; + t[5] += massone * vy * vz; } } @@ -372,7 +343,7 @@ void ComputeTempChunk::compute_vector() // final KE - MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(t, vector, 6, MPI_DOUBLE, MPI_SUM, world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } @@ -380,17 +351,7 @@ void ComputeTempChunk::compute_vector() void ComputeTempChunk::compute_array() { - 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(); - - if (nchunk > maxchunk) allocate(); - size_array_rows = nchunk; + ComputeChunk::compute_array(); // remove velocity bias @@ -409,9 +370,12 @@ void ComputeTempChunk::compute_array() // compute each value for (int i = 0; i < nvalues; i++) { - if (which[i] == TEMP) temperature(i); - else if (which[i] == KECOM) kecom(i); - else if (which[i] == INTERNAL) internal(i); + if (which[i] == TEMP) + temperature(i); + else if (which[i] == KECOM) + kecom(i); + else if (which[i] == INTERNAL) + internal(i); } // restore velocity bias @@ -425,7 +389,7 @@ void ComputeTempChunk::compute_array() void ComputeTempChunk::vcm_compute() { - int i,index; + int i, index; double massone; // avoid re-computing VCM more than once per step @@ -448,18 +412,20 @@ void ComputeTempChunk::vcm_compute() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; vcm[index][0] += v[i][0] * massone; vcm[index][1] += v[i][1] * massone; vcm[index][2] += v[i][2] * massone; massproc[index] += massone; } - MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); - MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&vcm[0][0], &vcmall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(massproc, masstotal, nchunk, MPI_DOUBLE, MPI_SUM, world); for (i = 0; i < nchunk; i++) { if (masstotal[i] > 0.0) { @@ -478,7 +444,7 @@ void ComputeTempChunk::vcm_compute() void ComputeTempChunk::temperature(int icol) { - int i,index; + int i, index; int *ichunk = cchunk->ichunk; // zero local per-chunk values @@ -501,65 +467,65 @@ void ComputeTempChunk::temperature(int icol) if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - sum[index] += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - rmass[i]; + sum[index] += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]) * rmass[i]; count[index]++; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - sum[index] += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - mass[type[i]]; + sum[index] += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]) * mass[type[i]]; count[index]++; - } + } } } else { - double vx,vy,vz; + double vx, vy, vz; if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; - sum[index] += (vx*vx + vy*vy + vz*vz) * rmass[i]; + sum[index] += (vx * vx + vy * vy + vz * vz) * rmass[i]; count[index]++; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; - sum[index] += (vx*vx + vy*vy + vz*vz) * mass[type[i]]; + sum[index] += (vx * vx + vy * vy + vz * vz) * mass[type[i]]; count[index]++; - } + } } } // sum across procs - MPI_Allreduce(sum,sumall,nchunk,MPI_DOUBLE,MPI_SUM,world); - MPI_Allreduce(count,countall,nchunk,MPI_INT,MPI_SUM,world); + MPI_Allreduce(sum, sumall, nchunk, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(count, countall, nchunk, MPI_INT, MPI_SUM, world); // normalize temperatures by per-chunk DOF - double dof,tfactor; + double dof, tfactor; double mvv2e = force->mvv2e; double boltz = force->boltz; for (i = 0; i < nchunk; i++) { - dof = cdof + adof*countall[i]; - if (dof > 0.0) tfactor = mvv2e / (dof * boltz); - else tfactor = 0.0; + dof = cdof + adof * countall[i]; + if (dof > 0.0) + tfactor = mvv2e / (dof * boltz); + else + tfactor = 0.0; array[i][icol] = tfactor * sumall[i]; } } @@ -585,37 +551,35 @@ void ComputeTempChunk::kecom(int icol) int *type = atom->type; int nlocal = atom->nlocal; - double vx,vy,vz; + double vx, vy, vz; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; vx = vcmall[index][0]; vy = vcmall[index][1]; vz = vcmall[index][2]; - sum[index] += (vx*vx + vy*vy + vz*vz) * rmass[i]; + sum[index] += (vx * vx + vy * vy + vz * vz) * rmass[i]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; vx = vcmall[index][0]; vy = vcmall[index][1]; vz = vcmall[index][2]; - sum[index] += (vx*vx + vy*vy + vz*vz) * mass[type[i]]; + sum[index] += (vx * vx + vy * vy + vz * vz) * mass[type[i]]; } } // sum across procs - MPI_Allreduce(sum,sumall,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(sum, sumall, nchunk, MPI_DOUBLE, MPI_SUM, world); double mvv2e = force->mvv2e; - for (int i = 0; i < nchunk; i++) - array[i][icol] = 0.5 * mvv2e * sumall[i]; - + for (int i = 0; i < nchunk; i++) array[i][icol] = 0.5 * mvv2e * sumall[i]; } /* ---------------------------------------------------------------------- @@ -641,36 +605,35 @@ void ComputeTempChunk::internal(int icol) int *type = atom->type; int nlocal = atom->nlocal; - double vx,vy,vz; + double vx, vy, vz; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; - sum[index] += (vx*vx + vy*vy + vz*vz) * rmass[i]; + sum[index] += (vx * vx + vy * vy + vz * vz) * rmass[i]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; vx = v[i][0] - vcmall[index][0]; vy = v[i][1] - vcmall[index][1]; vz = v[i][2] - vcmall[index][2]; - sum[index] += (vx*vx + vy*vy + vz*vz) * mass[type[i]]; + sum[index] += (vx * vx + vy * vy + vz * vz) * mass[type[i]]; } } // sum across procs - MPI_Allreduce(sum,sumall,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(sum, sumall, nchunk, MPI_DOUBLE, MPI_SUM, world); double mvv2e = force->mvv2e; - for (int i = 0; i < nchunk; i++) - array[i][icol] = 0.5 * mvv2e * sumall[i]; + for (int i = 0; i < nchunk; i++) array[i][icol] = 0.5 * mvv2e * sumall[i]; } /* ---------------------------------------------------------------------- @@ -683,7 +646,7 @@ void ComputeTempChunk::internal(int icol) void ComputeTempChunk::remove_bias(int i, double *v) { - int index = cchunk->ichunk[i]-1; + int index = cchunk->ichunk[i] - 1; if (index < 0) return; v[0] -= vcmall[index][0]; v[1] -= vcmall[index][1]; @@ -705,7 +668,7 @@ void ComputeTempChunk::remove_bias_all() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; v[i][0] -= vcmall[index][0]; v[i][1] -= vcmall[index][1]; @@ -720,7 +683,7 @@ void ComputeTempChunk::remove_bias_all() void ComputeTempChunk::restore_bias(int i, double *v) { - int index = cchunk->ichunk[i]-1; + int index = cchunk->ichunk[i] - 1; if (index < 0) return; v[0] += vcmall[index][0]; v[1] += vcmall[index][1]; @@ -743,7 +706,7 @@ void ComputeTempChunk::restore_bias_all() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; v[i][0] += vcmall[index][0]; v[i][1] += vcmall[index][1]; @@ -751,89 +714,34 @@ void ComputeTempChunk::restore_bias_all() } } -/* ---------------------------------------------------------------------- - 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 ComputeTempChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputeTempChunk::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 ComputeTempChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeTempChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeTempChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ void ComputeTempChunk::allocate() { + ComputeChunk::allocate(); memory->destroy(sum); memory->destroy(sumall); memory->destroy(count); memory->destroy(countall); memory->destroy(array); maxchunk = nchunk; - memory->create(sum,maxchunk,"temp/chunk:sum"); - memory->create(sumall,maxchunk,"temp/chunk:sumall"); - memory->create(count,maxchunk,"temp/chunk:count"); - memory->create(countall,maxchunk,"temp/chunk:countall"); - memory->create(array,maxchunk,nvalues,"temp/chunk:array"); + memory->create(sum, maxchunk, "temp/chunk:sum"); + memory->create(sumall, maxchunk, "temp/chunk:sumall"); + memory->create(count, maxchunk, "temp/chunk:count"); + memory->create(countall, maxchunk, "temp/chunk:countall"); + memory->create(array, maxchunk, nvalues, "temp/chunk:array"); if (comflag || nvalues) { memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(vcm); memory->destroy(vcmall); - memory->create(massproc,maxchunk,"vcm/chunk:massproc"); - memory->create(masstotal,maxchunk,"vcm/chunk:masstotal"); - memory->create(vcm,maxchunk,3,"vcm/chunk:vcm"); - memory->create(vcmall,maxchunk,3,"vcm/chunk:vcmall"); + memory->create(massproc, maxchunk, "vcm/chunk:massproc"); + memory->create(masstotal, maxchunk, "vcm/chunk:masstotal"); + memory->create(vcm, maxchunk, 3, "vcm/chunk:vcm"); + memory->create(vcmall, maxchunk, 3, "vcm/chunk:vcmall"); } } @@ -843,12 +751,12 @@ void ComputeTempChunk::allocate() double ComputeTempChunk::memory_usage() { - double bytes = (bigint) maxchunk * 2 * sizeof(double); + double bytes = (double) maxchunk * 2 * sizeof(double) + ComputeChunk::memory_usage(); bytes += (double) maxchunk * 2 * sizeof(int); bytes += (double) maxchunk * nvalues * sizeof(double); if (comflag || nvalues) { bytes += (double) maxchunk * 2 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); } return bytes; } diff --git a/src/compute_temp_chunk.h b/src/compute_temp_chunk.h index fc9b6d78fd..30e211b175 100644 --- a/src/compute_temp_chunk.h +++ b/src/compute_temp_chunk.h @@ -20,11 +20,11 @@ ComputeStyle(temp/chunk,ComputeTempChunk); #ifndef LMP_COMPUTE_TEMP_CHUNK_H #define LMP_COMPUTE_TEMP_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { -class ComputeTempChunk : public Compute { +class ComputeTempChunk : public ComputeChunk { public: ComputeTempChunk(class LAMMPS *, int, char **); ~ComputeTempChunk() override; @@ -38,20 +38,12 @@ class ComputeTempChunk : public Compute { void restore_bias(int, double *) override; void restore_bias_all() override; - void lock_enable() override; - void lock_disable() override; - int lock_length() override; - void lock(class Fix *, bigint, bigint) override; - void unlock(class Fix *) override; - double memory_usage() override; private: - int nchunk, maxchunk, comflag, biasflag; + int comflag, biasflag; int nvalues; int *which; - char *idchunk; - class ComputeChunkAtom *cchunk; double adof, cdof; char *id_bias; class Compute *tbias; // ptr to additional bias compute @@ -66,10 +58,8 @@ class ComputeTempChunk : public Compute { void temperature(int); void kecom(int); void internal(int); - void allocate(); + void allocate() override; }; - } // namespace LAMMPS_NS - #endif #endif diff --git a/src/compute_torque_chunk.cpp b/src/compute_torque_chunk.cpp index 36aca5dd5f..522ba5a03e 100644 --- a/src/compute_torque_chunk.cpp +++ b/src/compute_torque_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,24 +13,21 @@ #include "compute_torque_chunk.h" -#include #include "atom.h" -#include "update.h" -#include "modify.h" #include "compute_chunk_atom.h" #include "domain.h" -#include "memory.h" #include "error.h" +#include "memory.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTorqueChunk::ComputeTorqueChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - idchunk(nullptr), massproc(nullptr), masstotal(nullptr), com(nullptr), comall(nullptr), torque(nullptr), torqueall(nullptr) + ComputeChunk(lmp, narg, arg), massproc(nullptr), masstotal(nullptr), com(nullptr), + comall(nullptr), torque(nullptr), torqueall(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute torque/chunk command"); + if (narg != 4) error->all(FLERR, "Illegal compute torque/chunk command"); array_flag = 1; size_array_cols = 3; @@ -39,24 +35,14 @@ ComputeTorqueChunk::ComputeTorqueChunk(LAMMPS *lmp, int narg, char **arg) : size_array_rows_variable = 1; extarray = 0; - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); - ComputeTorqueChunk::init(); - - // chunk-based data - - nchunk = 1; - maxchunk = 0; - allocate(); + ComputeTorqueChunk::allocate(); } /* ---------------------------------------------------------------------- */ ComputeTorqueChunk::~ComputeTorqueChunk() { - delete [] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -67,38 +53,15 @@ ComputeTorqueChunk::~ComputeTorqueChunk() /* ---------------------------------------------------------------------- */ -void ComputeTorqueChunk::init() -{ - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for " - "compute torque/chunk"); - cchunk = dynamic_cast(modify->compute[icompute]); - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute torque/chunk does not use chunk/atom compute"); -} - -/* ---------------------------------------------------------------------- */ - void ComputeTorqueChunk::compute_array() { - int i,index; - double dx,dy,dz,massone; + int i, index; + double dx, dy, dz, 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(); + ComputeChunk::compute_array(); int *ichunk = cchunk->ichunk; - if (nchunk > maxchunk) allocate(); - size_array_rows = nchunk; - // zero local per-chunk values for (i = 0; i < nchunk; i++) { @@ -119,19 +82,21 @@ void ComputeTorqueChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - domain->unmap(x[i],image[i],unwrap); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + domain->unmap(x[i], image[i], unwrap); massproc[index] += massone; 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(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(massproc, masstotal, 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) { @@ -147,83 +112,26 @@ void ComputeTorqueChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - domain->unmap(x[i],image[i],unwrap); + domain->unmap(x[i], image[i], unwrap); dx = unwrap[0] - comall[index][0]; dy = unwrap[1] - comall[index][1]; dz = unwrap[2] - comall[index][2]; - torque[index][0] += dy*f[i][2] - dz*f[i][1]; - torque[index][1] += dz*f[i][0] - dx*f[i][2]; - torque[index][2] += dx*f[i][1] - dy*f[i][0]; + torque[index][0] += dy * f[i][2] - dz * f[i][1]; + torque[index][1] += dz * f[i][0] - dx * f[i][2]; + torque[index][2] += dx * f[i][1] - dy * f[i][0]; } - MPI_Allreduce(&torque[0][0],&torqueall[0][0],3*nchunk, - MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&torque[0][0], &torqueall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); } - -/* ---------------------------------------------------------------------- - 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 ComputeTorqueChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputeTorqueChunk::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 ComputeTorqueChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeTorqueChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeTorqueChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ void ComputeTorqueChunk::allocate() { + ComputeChunk::allocate(); memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -231,12 +139,12 @@ void ComputeTorqueChunk::allocate() memory->destroy(torque); memory->destroy(torqueall); maxchunk = nchunk; - memory->create(massproc,maxchunk,"torque/chunk:massproc"); - memory->create(masstotal,maxchunk,"torque/chunk:masstotal"); - memory->create(com,maxchunk,3,"torque/chunk:com"); - memory->create(comall,maxchunk,3,"torque/chunk:comall"); - memory->create(torque,maxchunk,3,"torque/chunk:torque"); - memory->create(torqueall,maxchunk,3,"torque/chunk:torqueall"); + memory->create(massproc, maxchunk, "torque/chunk:massproc"); + memory->create(masstotal, maxchunk, "torque/chunk:masstotal"); + memory->create(com, maxchunk, 3, "torque/chunk:com"); + memory->create(comall, maxchunk, 3, "torque/chunk:comall"); + memory->create(torque, maxchunk, 3, "torque/chunk:torque"); + memory->create(torqueall, maxchunk, 3, "torque/chunk:torqueall"); array = torqueall; } @@ -246,8 +154,8 @@ void ComputeTorqueChunk::allocate() double ComputeTorqueChunk::memory_usage() { - double bytes = (bigint) maxchunk * 2 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); + double bytes = (double) maxchunk * 2 * sizeof(double) + ComputeChunk::memory_usage(); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); return bytes; } diff --git a/src/compute_torque_chunk.h b/src/compute_torque_chunk.h index c810f7dc03..27752d6843 100644 --- a/src/compute_torque_chunk.h +++ b/src/compute_torque_chunk.h @@ -20,35 +20,24 @@ ComputeStyle(torque/chunk,ComputeTorqueChunk); #ifndef LMP_COMPUTE_TORQUE_CHUNK_H #define LMP_COMPUTE_TORQUE_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { -class ComputeTorqueChunk : public Compute { +class ComputeTorqueChunk : public ComputeChunk { public: ComputeTorqueChunk(class LAMMPS *, int, char **); ~ComputeTorqueChunk() override; - void init() override; void compute_array() override; - void lock_enable() override; - void lock_disable() override; - int lock_length() override; - void lock(class Fix *, bigint, bigint) override; - void unlock(class Fix *) override; - double memory_usage() override; private: - int nchunk, maxchunk; - char *idchunk; - class ComputeChunkAtom *cchunk; - double *massproc, *masstotal; double **com, **comall; double **torque, **torqueall; - void allocate(); + void allocate() override; }; } // namespace LAMMPS_NS diff --git a/src/compute_vcm_chunk.cpp b/src/compute_vcm_chunk.cpp index c1d22bf5e6..07b504f2ae 100644 --- a/src/compute_vcm_chunk.cpp +++ b/src/compute_vcm_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,25 +13,22 @@ #include "compute_vcm_chunk.h" -#include #include "atom.h" -#include "update.h" -#include "modify.h" #include "compute_chunk_atom.h" -#include "memory.h" #include "error.h" +#include "memory.h" using namespace LAMMPS_NS; -enum{ONCE,NFREQ,EVERY}; +enum { ONCE, NFREQ, EVERY }; /* ---------------------------------------------------------------------- */ ComputeVCMChunk::ComputeVCMChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - idchunk(nullptr), massproc(nullptr), masstotal(nullptr), vcm(nullptr), vcmall(nullptr) + ComputeChunk(lmp, narg, arg), massproc(nullptr), masstotal(nullptr), vcm(nullptr), + vcmall(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute vcm/chunk command"); + if (narg != 4) error->all(FLERR, "Illegal compute vcm/chunk command"); array_flag = 1; size_array_cols = 3; @@ -40,26 +36,14 @@ ComputeVCMChunk::ComputeVCMChunk(LAMMPS *lmp, int narg, char **arg) : size_array_rows_variable = 1; extarray = 0; - // ID of compute chunk/atom - - idchunk = utils::strdup(arg[3]); - ComputeVCMChunk::init(); - - // chunk-based data - - nchunk = 1; - maxchunk = 0; - allocate(); - - firstflag = massneed = 1; + ComputeVCMChunk::allocate(); } /* ---------------------------------------------------------------------- */ ComputeVCMChunk::~ComputeVCMChunk() { - delete [] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(vcm); @@ -68,18 +52,6 @@ ComputeVCMChunk::~ComputeVCMChunk() /* ---------------------------------------------------------------------- */ -void ComputeVCMChunk::init() -{ - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for compute vcm/chunk"); - cchunk = dynamic_cast(modify->compute[icompute]); - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute vcm/chunk does not use chunk/atom compute"); -} - -/* ---------------------------------------------------------------------- */ - void ComputeVCMChunk::setup() { // one-time calculation of per-chunk mass @@ -98,23 +70,12 @@ void ComputeVCMChunk::compute_array() int index; double massone; - 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(); + ComputeChunk::compute_array(); int *ichunk = cchunk->ichunk; - if (nchunk > maxchunk) allocate(); - size_array_rows = nchunk; - // zero local per-chunk values - for (int i = 0; i < nchunk; i++) - vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0; + for (int i = 0; i < nchunk; i++) vcm[i][0] = vcm[i][1] = vcm[i][2] = 0.0; if (massneed) for (int i = 0; i < nchunk; i++) massproc[i] = 0.0; @@ -129,100 +90,47 @@ void ComputeVCMChunk::compute_array() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; vcm[index][0] += v[i][0] * massone; vcm[index][1] += v[i][1] * massone; vcm[index][2] += v[i][2] * massone; if (massneed) massproc[index] += massone; } - MPI_Allreduce(&vcm[0][0],&vcmall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); - if (massneed) - MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&vcm[0][0], &vcmall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); + if (massneed) MPI_Allreduce(massproc, masstotal, nchunk, MPI_DOUBLE, MPI_SUM, world); for (int i = 0; i < nchunk; i++) { if (masstotal[i] > 0.0) { vcmall[i][0] /= masstotal[i]; vcmall[i][1] /= masstotal[i]; vcmall[i][2] /= masstotal[i]; - } else vcmall[i][0] = vcmall[i][1] = vcmall[i][2] = 0.0; + } else + vcmall[i][0] = vcmall[i][1] = vcmall[i][2] = 0.0; } } -/* ---------------------------------------------------------------------- - 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 ComputeVCMChunk::lock_enable() -{ - cchunk->lockcount++; -} - -/* ---------------------------------------------------------------------- - decrement lock counter in compute chunk/atom, it if still exists -------------------------------------------------------------------------- */ - -void ComputeVCMChunk::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 ComputeVCMChunk::lock_length() -{ - nchunk = cchunk->setup_chunks(); - return nchunk; -} - -/* ---------------------------------------------------------------------- - set the lock from startstep to stopstep -------------------------------------------------------------------------- */ - -void ComputeVCMChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) -{ - cchunk->lock(fixptr,startstep,stopstep); -} - -/* ---------------------------------------------------------------------- - unset the lock -------------------------------------------------------------------------- */ - -void ComputeVCMChunk::unlock(Fix *fixptr) -{ - cchunk->unlock(fixptr); -} - /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ void ComputeVCMChunk::allocate() { + ComputeChunk::allocate(); memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(vcm); memory->destroy(vcmall); maxchunk = nchunk; - memory->create(massproc,maxchunk,"vcm/chunk:massproc"); - memory->create(masstotal,maxchunk,"vcm/chunk:masstotal"); - memory->create(vcm,maxchunk,3,"vcm/chunk:vcm"); - memory->create(vcmall,maxchunk,3,"vcm/chunk:vcmall"); + memory->create(massproc, maxchunk, "vcm/chunk:massproc"); + memory->create(masstotal, maxchunk, "vcm/chunk:masstotal"); + memory->create(vcm, maxchunk, 3, "vcm/chunk:vcm"); + memory->create(vcmall, maxchunk, 3, "vcm/chunk:vcmall"); array = vcmall; } @@ -232,7 +140,7 @@ void ComputeVCMChunk::allocate() double ComputeVCMChunk::memory_usage() { - double bytes = (bigint) maxchunk * 2 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); + double bytes = (double) maxchunk * 2 * sizeof(double) + ComputeChunk::memory_usage(); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); return bytes; } diff --git a/src/compute_vcm_chunk.h b/src/compute_vcm_chunk.h index 418ae63c3a..af94a235d9 100644 --- a/src/compute_vcm_chunk.h +++ b/src/compute_vcm_chunk.h @@ -20,39 +20,26 @@ ComputeStyle(vcm/chunk,ComputeVCMChunk); #ifndef LMP_COMPUTE_VCM_CHUNK_H #define LMP_COMPUTE_VCM_CHUNK_H -#include "compute.h" +#include "compute_chunk.h" namespace LAMMPS_NS { -class ComputeVCMChunk : public Compute { +class ComputeVCMChunk : public ComputeChunk { public: ComputeVCMChunk(class LAMMPS *, int, char **); ~ComputeVCMChunk() override; - void init() override; + void setup() override; void compute_array() override; - void lock_enable() override; - void lock_disable() override; - int lock_length() override; - void lock(class Fix *, bigint, bigint) override; - void unlock(class Fix *) override; - double memory_usage() override; private: - int nchunk, maxchunk; - int firstflag, massneed; - char *idchunk; - class ComputeChunkAtom *cchunk; - double *massproc, *masstotal; double **vcm, **vcmall; - void allocate(); + void allocate() override; }; - } // namespace LAMMPS_NS - #endif #endif