ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Introduce ComputeChunk class with shared functionality of all /chunk computes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-03-18 05:55:03 -04:00
parent fce1f8e0af
commit 1ccb0f8d8d
29 changed files with 958 additions and 2050 deletions
Download Patch File Download Diff File Expand all files Collapse all files

190
src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.cpp
View File

@ -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 <cmath>
#include <cstring>
@ -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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(modify->compute[icompute]);
cchunk->lockcount--;
}
}
/* ----------------------------------------------------------------------
calculate and return # of chunks = length of vector/array
------------------------------------------------------------------------- */
int ComputeDipoleTIP4PChunk::lock_length()
{
nchunk = cchunk->setup_chunks();
return nchunk;
}
/* ----------------------------------------------------------------------
set the lock from startstep to stopstep
------------------------------------------------------------------------- */
void ComputeDipoleTIP4PChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep)
{
cchunk->lock(fixptr,startstep,stopstep);
}
/* ----------------------------------------------------------------------
unset the lock
------------------------------------------------------------------------- */
void ComputeDipoleTIP4PChunk::unlock(Fix *fixptr)
{
cchunk->unlock(fixptr);
}
/* ----------------------------------------------------------------------
free and reallocate per-chunk arrays
------------------------------------------------------------------------- */
void ComputeDipoleTIP4PChunk::allocate()
{
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];

18
src/EXTRA-COMPUTE/compute_dipole_tip4p_chunk.h
View File

@ -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

55
src/EXTRA-COMPUTE/compute_gyration_shape_chunk.cpp
View File

@ -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;
}

98
src/compute_angmom_chunk.cpp
View File

@ -18,17 +18,13 @@
#include "domain.h"
#include "error.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include <cstring>
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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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;

20
src/compute_angmom_chunk.h
View File

@ -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

156
src/compute_chunk.cpp Normal file
View File

@ -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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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);
}

49
src/compute_chunk.h Normal file
View File

@ -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

141
src/compute_com_chunk.cpp
View File

@ -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 <cstring>
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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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;
}

22
src/compute_com_chunk.h
View File

@ -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

184
src/compute_dipole_chunk.cpp
View File

@ -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 <cmath>
#include <cstring>
@ -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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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;
}

18
src/compute_dipole_chunk.h
View File

@ -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

212
src/compute_gyration_chunk.cpp
View File

@ -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 <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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;
}

18
src/compute_gyration_chunk.h
View File

@ -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

170
src/compute_inertia_chunk.cpp
View File

@ -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 <cstring>
#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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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;
}

19
src/compute_inertia_chunk.h
View File

@ -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

143
src/compute_msd_chunk.cpp
View File

@ -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 <cstring>
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<FixStoreGlobal *>(
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<ComputeChunkAtom *>(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<FixStoreGlobal *>(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<ComputeChunkAtom *>(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;
}

16
src/compute_msd_chunk.h
View File

@ -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

272
src/compute_omega_chunk.cpp
View File

@ -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 <cstring>
#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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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;
}

20
src/compute_omega_chunk.h
View File

@ -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

177
src/compute_property_chunk.cpp
View File

@ -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 <cstring>
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "compute_chunk_atom.h"
#include "memory.h"
#include "error.h"
#include "memory.h"
#include <cstring>
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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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];

20
src/compute_property_chunk.h
View File

@ -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

231
src/compute_reduce_chunk.cpp
View File

@ -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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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;
}

17
src/compute_reduce_chunk.h
View File

@ -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<value_t> 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);
};

390
src/compute_temp_chunk.cpp
View File

@ -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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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;
}

18
src/compute_temp_chunk.h
View File

@ -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

156
src/compute_torque_chunk.cpp
View File

@ -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 <cstring>
#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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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;
}

17
src/compute_torque_chunk.h
View File

@ -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

140
src/compute_vcm_chunk.cpp
View File

@ -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 <cstring>
#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<ComputeChunkAtom *>(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<ComputeChunkAtom *>(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;
}

21
src/compute_vcm_chunk.h
View File

@ -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