tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

enable and apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-06-17 13:47:23 -04:00
parent d9646ee537
commit 9fd6bde0ed
8 changed files with 337 additions and 378 deletions
Download Patch File Download Diff File Expand all files Collapse all files

200
src/ML-SNAP/compute_sna_grid.cpp
View File

@ -11,22 +11,22 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "compute_grid.h"
#include "compute_sna_grid.h"
#include "sna.h"
#include "atom.h"
#include "update.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "neighbor.h"
#include "pair.h"
#include "domain.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "sna.h"
#include "tokenizer.h"
#include "update.h"
#include <cmath>
#include <cstring>
@ -34,8 +34,7 @@
using namespace LAMMPS_NS;
ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) :
ComputeGrid(lmp, narg, arg), cutsq(nullptr),
radelem(nullptr), wjelem(nullptr)
ComputeGrid(lmp, narg, arg), cutsq(nullptr), radelem(nullptr), wjelem(nullptr)
{
double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
@ -48,9 +47,9 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) :
narg -= nargbase;
int ntypes = atom->ntypes;
int nargmin = 6+2*ntypes;
int nargmin = 6 + 2 * ntypes;
if (narg < nargmin) error->all(FLERR,"Illegal compute sna/grid command");
if (narg < nargmin) error->all(FLERR, "Illegal compute sna/grid command");
// default values
@ -66,30 +65,28 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) :
// process required arguments
memory->create(radelem,ntypes+1,"sna/grid:radelem"); // offset by 1 to match up with types
memory->create(wjelem,ntypes+1,"sna/grid:wjelem");
memory->create(radelem, ntypes + 1, "sna/grid:radelem"); // offset by 1 to match up with types
memory->create(wjelem, ntypes + 1, "sna/grid:wjelem");
rcutfac = atof(arg[3]);
rfac0 = atof(arg[4]);
twojmax = atoi(arg[5]);
for(int i = 0; i < ntypes; i++)
radelem[i+1] = atof(arg[6+i]);
for(int i = 0; i < ntypes; i++)
wjelem[i+1] = atof(arg[6+ntypes+i]);
for (int i = 0; i < ntypes; i++) radelem[i + 1] = atof(arg[6 + i]);
for (int i = 0; i < ntypes; i++) wjelem[i + 1] = atof(arg[6 + ntypes + i]);
// construct cutsq
double cut;
cutmax = 0.0;
memory->create(cutsq,ntypes+1,ntypes+1,"sna/grid:cutsq");
for(int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/grid:cutsq");
for (int i = 1; i <= ntypes; i++) {
cut = 2.0 * radelem[i] * rcutfac;
if (cut > cutmax) cutmax = cut;
cutsq[i][i] = cut*cut;
for(int j = i+1; j <= ntypes; j++) {
cut = (radelem[i]+radelem[j])*rcutfac;
cutsq[i][j] = cutsq[j][i] = cut*cut;
cutsq[i][i] = cut * cut;
for (int j = i + 1; j <= ntypes; j++) {
cut = (radelem[i] + radelem[j]) * rcutfac;
cutsq[i][j] = cutsq[j][i] = cut * cut;
}
}
@ -103,90 +100,81 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) :
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid command");
rmin0 = atof(arg[iarg+1]);
if (strcmp(arg[iarg], "rmin0") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid command");
rmin0 = atof(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"switchflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid command");
switchflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid command");
switchflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"bzeroflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid command");
bzeroflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "bzeroflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid command");
bzeroflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"quadraticflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid command");
quadraticflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "quadraticflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid command");
quadraticflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid command");
} else if (strcmp(arg[iarg], "chem") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid command");
chemflag = 1;
memory->create(map,ntypes+1,"compute_sna_grid:map");
nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
memory->create(map, ntypes + 1, "compute_sna_grid:map");
nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
for (int i = 0; i < ntypes; i++) {
int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp);
if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute sna/grid command");
map[i+1] = jelem;
int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp);
if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute sna/grid command");
map[i + 1] = jelem;
}
iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2 + ntypes;
} else if (strcmp(arg[iarg], "bnormflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid command");
bnormflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid command");
wselfallflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "wselfallflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid command");
wselfallflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid command");
switchinnerflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchinnerflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid command");
switchinnerflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"sinner") == 0) {
} else if (strcmp(arg[iarg], "sinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute sna/grid command");
memory->create(sinnerelem,ntypes+1,"snap:sinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute sna/grid command");
memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem");
for (int i = 0; i < ntypes; i++)
sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
sinnerflag = 1;
iarg += ntypes;
} else if (strcmp(arg[iarg],"dinner") == 0) {
} else if (strcmp(arg[iarg], "dinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute sna/grid command");
memory->create(dinnerelem,ntypes+1,"snap:dinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute sna/grid command");
memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem");
for (int i = 0; i < ntypes; i++)
dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
dinnerflag = 1;
iarg += ntypes;
} else error->all(FLERR,"Illegal compute sna/grid command");
} else
error->all(FLERR, "Illegal compute sna/grid command");
}
if (switchinnerflag && !(sinnerflag && dinnerflag))
error->all(FLERR,"Illegal compute sna/grid command: switchinnerflag = 1, missing sinner/dinner keyword");
error->all(
FLERR,
"Illegal compute sna/grid command: switchinnerflag = 1, missing sinner/dinner keyword");
if (!switchinnerflag && (sinnerflag || dinnerflag))
error->all(FLERR,"Illegal compute sna/grid command: switchinnerflag = 0, unexpected sinner/dinner keyword");
error->all(
FLERR,
"Illegal compute sna/grid command: switchinnerflag = 0, unexpected sinner/dinner keyword");
snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag,
chemflag, bnormflag, wselfallflag,
nelements, switchinnerflag);
snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag,
wselfallflag, nelements, switchinnerflag);
ncoeff = snaptr->ncoeff;
nvalues = ncoeff;
if (quadraticflag) nvalues += (ncoeff*(ncoeff+1))/2;
if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2;
size_array_cols = size_array_cols_base + nvalues;
array_flag = 1;
}
@ -209,9 +197,8 @@ void ComputeSNAGrid::init()
{
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"sna/grid") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute sna/grid");
if (strcmp(modify->compute[i]->style, "sna/grid") == 0) count++;
if (count > 1 && comm->me == 0) error->warning(FLERR, "More than one compute sna/grid");
snaptr->init();
}
@ -230,9 +217,9 @@ void ComputeSNAGrid::compute_array()
// compute sna for each gridpoint
double** const x = atom->x;
const int* const mask = atom->mask;
int * const type = atom->type;
double **const x = atom->x;
const int *const mask = atom->mask;
int *const type = atom->type;
const int ntotal = atom->nlocal + atom->nghost;
// insure rij, inside, and typej are of size jnum
@ -243,7 +230,7 @@ void ComputeSNAGrid::compute_array()
for (int iy = nylo; iy <= nyhi; iy++)
for (int ix = nxlo; ix <= nxhi; ix++) {
double xgrid[3];
const int igrid = iz*(nx*ny) + iy*nx + ix;
const int igrid = iz * (nx * ny) + iy * nx + ix;
grid2x(igrid, xgrid);
const double xtmp = xgrid[0];
const double ytmp = xgrid[1];
@ -253,8 +240,7 @@ void ComputeSNAGrid::compute_array()
const int itype = 1;
int ielem = 0;
if (chemflag)
ielem = map[itype];
if (chemflag) ielem = map[itype];
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
@ -270,11 +256,10 @@ void ComputeSNAGrid::compute_array()
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx*delx + dely*dely + delz*delz;
const double rsq = delx * delx + dely * dely + delz * delz;
int jtype = type[j];
int jelem = 0;
if (chemflag)
jelem = map[jtype];
if (chemflag) jelem = map[jtype];
if (rsq < cutsq[jtype][jtype] && rsq > 1e-20) {
snaptr->rij[ninside][0] = delx;
@ -282,7 +267,7 @@ void ComputeSNAGrid::compute_array()
snaptr->rij[ninside][2] = delz;
snaptr->inside[ninside] = j;
snaptr->wj[ninside] = wjelem[jtype];
snaptr->rcutij[ninside] = 2.0*radelem[jtype]*rcutfac;
snaptr->rcutij[ninside] = 2.0 * radelem[jtype] * rcutfac;
if (switchinnerflag) {
snaptr->sinnerij[ninside] = sinnerelem[jelem];
snaptr->dinnerij[ninside] = dinnerelem[jelem];
@ -299,7 +284,7 @@ void ComputeSNAGrid::compute_array()
// linear contributions
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
gridlocal[size_array_cols_base+icoeff][iz][iy][ix] = snaptr->blist[icoeff];
gridlocal[size_array_cols_base + icoeff][iz][iy][ix] = snaptr->blist[icoeff];
// quadratic contributions
@ -307,28 +292,28 @@ void ComputeSNAGrid::compute_array()
int ncount = ncoeff;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bveci = snaptr->blist[icoeff];
gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = 0.5*bveci*bveci;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++)
gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = bveci*snaptr->blist[jcoeff];
gridlocal[size_array_cols_base + ncount++][iz][iy][ix] = 0.5 * bveci * bveci;
for (int jcoeff = icoeff + 1; jcoeff < ncoeff; jcoeff++)
gridlocal[size_array_cols_base + ncount++][iz][iy][ix] =
bveci * snaptr->blist[jcoeff];
}
}
}
memset(&grid[0][0],0,size_array_rows*size_array_cols*sizeof(double));
memset(&grid[0][0], 0, size_array_rows * size_array_cols * sizeof(double));
for (int iz = nzlo; iz <= nzhi; iz++)
for (int iy = nylo; iy <= nyhi; iy++)
for (int ix = nxlo; ix <= nxhi; ix++) {
const int igrid = iz*(nx*ny) + iy*nx + ix;
const int igrid = iz * (nx * ny) + iy * nx + ix;
for (int j = 0; j < nvalues; j++)
grid[igrid][size_array_cols_base + j] = gridlocal[size_array_cols_base + j][iz][iy][ix];
}
MPI_Allreduce(&grid[0][0],&gridall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&grid[0][0], &gridall[0][0], size_array_rows * size_array_cols, MPI_DOUBLE, MPI_SUM,
world);
assign_coords_all();
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
@ -336,9 +321,8 @@ void ComputeSNAGrid::compute_array()
double ComputeSNAGrid::memory_usage()
{
double nbytes = snaptr->memory_usage(); // SNA object
int n = atom->ntypes+1;
nbytes += (double)n*sizeof(int); // map
int n = atom->ntypes + 1;
nbytes += (double) n * sizeof(int); // map
return nbytes;
}

2
src/ML-SNAP/compute_sna_grid.h
View File

@ -50,7 +50,7 @@ class ComputeSNAGrid : public ComputeGrid {
int quadraticflag;
};
}
} // namespace LAMMPS_NS
#endif
#endif

193
src/ML-SNAP/compute_sna_grid_local.cpp
View File

@ -11,22 +11,22 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "compute_grid_local.h"
#include "compute_sna_grid_local.h"
#include "sna.h"
#include "atom.h"
#include "update.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "neighbor.h"
#include "pair.h"
#include "domain.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "sna.h"
#include "tokenizer.h"
#include "update.h"
#include <cmath>
#include <cstring>
@ -34,8 +34,7 @@
using namespace LAMMPS_NS;
ComputeSNAGridLocal::ComputeSNAGridLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeGridLocal(lmp, narg, arg), cutsq(nullptr),
radelem(nullptr), wjelem(nullptr)
ComputeGridLocal(lmp, narg, arg), cutsq(nullptr), radelem(nullptr), wjelem(nullptr)
{
double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
@ -48,9 +47,9 @@ ComputeSNAGridLocal::ComputeSNAGridLocal(LAMMPS *lmp, int narg, char **arg) :
narg -= nargbase;
int ntypes = atom->ntypes;
int nargmin = 6+2*ntypes;
int nargmin = 6 + 2 * ntypes;
if (narg < nargmin) error->all(FLERR,"Illegal compute sna/grid/local command");
if (narg < nargmin) error->all(FLERR, "Illegal compute sna/grid/local command");
// default values
@ -71,30 +70,29 @@ ComputeSNAGridLocal::ComputeSNAGridLocal(LAMMPS *lmp, int narg, char **arg) :
// process required arguments
memory->create(radelem,ntypes+1,"sna/grid/local:radelem"); // offset by 1 to match up with types
memory->create(wjelem,ntypes+1,"sna/grid/local:wjelem");
memory->create(radelem, ntypes + 1,
"sna/grid/local:radelem"); // offset by 1 to match up with types
memory->create(wjelem, ntypes + 1, "sna/grid/local:wjelem");
rcutfac = atof(arg[3]);
rfac0 = atof(arg[4]);
twojmax = atoi(arg[5]);
for(int i = 0; i < ntypes; i++)
radelem[i+1] = atof(arg[6+i]);
for(int i = 0; i < ntypes; i++)
wjelem[i+1] = atof(arg[6+ntypes+i]);
for (int i = 0; i < ntypes; i++) radelem[i + 1] = atof(arg[6 + i]);
for (int i = 0; i < ntypes; i++) wjelem[i + 1] = atof(arg[6 + ntypes + i]);
// construct cutsq
double cut;
cutmax = 0.0;
memory->create(cutsq,ntypes+1,ntypes+1,"sna/grid/local:cutsq");
for(int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/grid/local:cutsq");
for (int i = 1; i <= ntypes; i++) {
cut = 2.0 * radelem[i] * rcutfac;
if (cut > cutmax) cutmax = cut;
cutsq[i][i] = cut*cut;
for(int j = i+1; j <= ntypes; j++) {
cut = (radelem[i]+radelem[j])*rcutfac;
cutsq[i][j] = cutsq[j][i] = cut*cut;
cutsq[i][i] = cut * cut;
for (int j = i + 1; j <= ntypes; j++) {
cut = (radelem[i] + radelem[j]) * rcutfac;
cutsq[i][j] = cutsq[j][i] = cut * cut;
}
}
@ -103,90 +101,82 @@ ComputeSNAGridLocal::ComputeSNAGridLocal(LAMMPS *lmp, int narg, char **arg) :
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid/local command");
rmin0 = atof(arg[iarg+1]);
if (strcmp(arg[iarg], "rmin0") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid/local command");
rmin0 = atof(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"switchflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid/local command");
switchflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid/local command");
switchflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"bzeroflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid/local command");
bzeroflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "bzeroflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid/local command");
bzeroflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"quadraticflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid/local command");
quadraticflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "quadraticflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid/local command");
quadraticflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid/local command");
} else if (strcmp(arg[iarg], "chem") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid/local command");
chemflag = 1;
memory->create(map,ntypes+1,"compute_sna_grid_local:map");
nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
memory->create(map, ntypes + 1, "compute_sna_grid_local:map");
nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
for (int i = 0; i < ntypes; i++) {
int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp);
int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp);
if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute sna/grid/local command");
map[i+1] = jelem;
error->all(FLERR, "Illegal compute sna/grid/local command");
map[i + 1] = jelem;
}
iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid/local command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2 + ntypes;
} else if (strcmp(arg[iarg], "bnormflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid/local command");
bnormflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid/local command");
wselfallflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "wselfallflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid/local command");
wselfallflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/grid/local command");
switchinnerflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchinnerflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute sna/grid/local command");
switchinnerflag = atoi(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg],"sinner") == 0) {
} else if (strcmp(arg[iarg], "sinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute sna/grid/local command");
memory->create(sinnerelem,ntypes+1,"snap:sinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute sna/grid/local command");
memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem");
for (int i = 0; i < ntypes; i++)
sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
sinnerflag = 1;
iarg += ntypes;
} else if (strcmp(arg[iarg],"dinner") == 0) {
} else if (strcmp(arg[iarg], "dinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute sna/grid/local command");
memory->create(dinnerelem,ntypes+1,"snap:dinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute sna/grid/local command");
memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem");
for (int i = 0; i < ntypes; i++)
dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
dinnerflag = 1;
iarg += ntypes;
} else error->all(FLERR,"Illegal compute sna/grid/local command");
} else
error->all(FLERR, "Illegal compute sna/grid/local command");
}
if (switchinnerflag && !(sinnerflag && dinnerflag))
error->all(FLERR,"Illegal compute sna/grid/local command: switchinnerflag = 1, missing sinner/dinner keyword");
error->all(FLERR,
"Illegal compute sna/grid/local command: switchinnerflag = 1, missing sinner/dinner "
"keyword");
if (!switchinnerflag && (sinnerflag || dinnerflag))
error->all(FLERR,"Illegal compute sna/grid/local command: switchinnerflag = 0, unexpected sinner/dinner keyword");
error->all(FLERR,
"Illegal compute sna/grid/local command: switchinnerflag = 0, unexpected "
"sinner/dinner keyword");
snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag,
chemflag, bnormflag, wselfallflag,
nelements, switchinnerflag);
snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag,
wselfallflag, nelements, switchinnerflag);
ncoeff = snaptr->ncoeff;
nvalues = ncoeff;
if (quadraticflag) nvalues += (ncoeff*(ncoeff+1))/2;
if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2;
size_local_cols = size_local_cols_base + nvalues;
}
@ -208,9 +198,8 @@ void ComputeSNAGridLocal::init()
{
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"sna/grid/local") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute sna/grid/local");
if (strcmp(modify->compute[i]->style, "sna/grid/local") == 0) count++;
if (count > 1 && comm->me == 0) error->warning(FLERR, "More than one compute sna/grid/local");
snaptr->init();
}
@ -229,9 +218,9 @@ void ComputeSNAGridLocal::compute_local()
// compute sna for each gridpoint
double** const x = atom->x;
const int* const mask = atom->mask;
int * const type = atom->type;
double **const x = atom->x;
const int *const mask = atom->mask;
int *const type = atom->type;
const int ntotal = atom->nlocal + atom->nghost;
// insure rij, inside, and typej are of size jnum
@ -252,8 +241,7 @@ void ComputeSNAGridLocal::compute_local()
const int itype = 1;
int ielem = 0;
if (chemflag)
ielem = map[itype];
if (chemflag) ielem = map[itype];
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
@ -269,23 +257,23 @@ void ComputeSNAGridLocal::compute_local()
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx*delx + dely*dely + delz*delz;
const double rsq = delx * delx + dely * dely + delz * delz;
int jtype = type[j];
int jelem = 0;
if (chemflag)
jelem = map[jtype];
if (chemflag) jelem = map[jtype];
if (rsq < cutsq[jtype][jtype] && rsq > 1e-20) {
snaptr->rij[ninside][0] = delx;
snaptr->rij[ninside][1] = dely;
snaptr->rij[ninside][2] = delz;
snaptr->inside[ninside] = j;
snaptr->wj[ninside] = wjelem[jtype];
snaptr->rcutij[ninside] = 2.0*radelem[jtype]*rcutfac;
snaptr->rcutij[ninside] = 2.0 * radelem[jtype] * rcutfac;
if (switchinnerflag) {
snaptr->sinnerij[ninside] = sinnerelem[jelem];
snaptr->dinnerij[ninside] = dinnerelem[jelem];
}
if (chemflag) snaptr->element[ninside] = jelem; // element index for multi-element snap
if (chemflag)
snaptr->element[ninside] = jelem; // element index for multi-element snap
ninside++;
}
}
@ -297,7 +285,7 @@ void ComputeSNAGridLocal::compute_local()
// linear contributions
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
alocal[igrid][size_local_cols_base+icoeff] = snaptr->blist[icoeff];
alocal[igrid][size_local_cols_base + icoeff] = snaptr->blist[icoeff];
// quadratic contributions
@ -305,17 +293,15 @@ void ComputeSNAGridLocal::compute_local()
int ncount = ncoeff;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bveci = snaptr->blist[icoeff];
alocal[igrid][size_local_cols_base+ncount++] = 0.5*bveci*bveci;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++)
alocal[igrid][size_local_cols_base+ncount++] =
bveci*snaptr->blist[jcoeff];
alocal[igrid][size_local_cols_base + ncount++] = 0.5 * bveci * bveci;
for (int jcoeff = icoeff + 1; jcoeff < ncoeff; jcoeff++)
alocal[igrid][size_local_cols_base + ncount++] = bveci * snaptr->blist[jcoeff];
}
}
igrid++;
}
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
@ -323,9 +309,8 @@ void ComputeSNAGridLocal::compute_local()
double ComputeSNAGridLocal::memory_usage()
{
double nbytes = snaptr->memory_usage(); // SNA object
int n = atom->ntypes+1;
nbytes += (double)n*sizeof(int); // map
int n = atom->ntypes + 1;
nbytes += (double) n * sizeof(int); // map
return nbytes;
}

2
src/ML-SNAP/compute_sna_grid_local.h
View File

@ -50,7 +50,7 @@ class ComputeSNAGridLocal : public ComputeGridLocal {
int quadraticflag;
};
}
} // namespace LAMMPS_NS
#endif
#endif

121
src/compute_grid.cpp
View File

@ -12,25 +12,26 @@
------------------------------------------------------------------------- */
#include "compute_grid.h"
#include <mpi.h>
#include <cstring>
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "comm.h"
#include "modify.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), grid(nullptr), gridall(nullptr), gridlocal(nullptr)
Compute(lmp, narg, arg), grid(nullptr), gridall(nullptr), gridlocal(nullptr)
{
if (narg < 6) error->all(FLERR,"Illegal compute grid command");
if (narg < 6) error->all(FLERR, "Illegal compute grid command");
array_flag = 1;
size_array_cols = 0;
@ -39,19 +40,19 @@ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) :
int iarg0 = 3;
int iarg = iarg0;
if (strcmp(arg[iarg],"grid") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal compute grid command");
nx = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
ny = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
nz = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
if (nx <= 0 || ny <= 0 || nz <= 0)
error->all(FLERR,"All grid dimensions must be positive");
if (strcmp(arg[iarg], "grid") == 0) {
if (iarg + 4 > narg) error->all(FLERR, "Illegal compute grid command");
nx = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
ny = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
nz = utils::inumeric(FLERR, arg[iarg + 3], false, lmp);
if (nx <= 0 || ny <= 0 || nz <= 0) error->all(FLERR, "All grid dimensions must be positive");
iarg += 4;
} else error->all(FLERR,"Illegal compute grid command");
} else
error->all(FLERR, "Illegal compute grid command");
nargbase = iarg - iarg0;
size_array_rows = nx*ny*nz;
size_array_rows = nx * ny * nz;
size_array_cols_base = 3;
gridlocal_allocated = 0;
}
@ -65,9 +66,7 @@ ComputeGrid::~ComputeGrid()
/* ---------------------------------------------------------------------- */
void ComputeGrid::init()
{
}
void ComputeGrid::init() {}
/* ---------------------------------------------------------------------- */
@ -85,18 +84,17 @@ void ComputeGrid::setup()
void ComputeGrid::grid2x(int igrid, double *x)
{
int iz = igrid / (nx*ny);
igrid -= iz * (nx*ny);
int iz = igrid / (nx * ny);
igrid -= iz * (nx * ny);
int iy = igrid / nx;
igrid -= iy * nx;
int ix = igrid;
x[0] = ix*delx;
x[1] = iy*dely;
x[2] = iz*delz;
x[0] = ix * delx;
x[1] = iy * dely;
x[2] = iz * delz;
if (triclinic) domain->lamda2x(x, x);
}
/* ----------------------------------------------------------------------
@ -107,7 +105,7 @@ void ComputeGrid::assign_coords_all()
{
double x[3];
for (int igrid = 0; igrid < size_array_rows; igrid++) {
grid2x(igrid,x);
grid2x(igrid, x);
gridall[igrid][0] = x[0];
gridall[igrid][1] = x[1];
gridall[igrid][2] = x[2];
@ -122,18 +120,17 @@ void ComputeGrid::allocate()
{
// allocate arrays
printf("In allocate() %d %d \n", size_array_rows,size_array_cols);
memory->create(grid,size_array_rows,size_array_cols,"grid:grid");
memory->create(gridall,size_array_rows,size_array_cols,"grid:gridall");
printf("In allocate() %d %d \n", size_array_rows, size_array_cols);
memory->create(grid, size_array_rows, size_array_cols, "grid:grid");
memory->create(gridall, size_array_rows, size_array_cols, "grid:gridall");
if (nxlo <= nxhi && nylo <= nyhi && nzlo <= nzhi) {
gridlocal_allocated = 1;
memory->create4d_offset(gridlocal,size_array_cols,nzlo,nzhi,nylo,nyhi,
nxlo,nxhi,"grid:gridlocal");
memory->create4d_offset(gridlocal, size_array_cols, nzlo, nzhi, nylo, nyhi, nxlo, nxhi,
"grid:gridlocal");
}
array = gridall;
}
/* ----------------------------------------------------------------------
free arrays
------------------------------------------------------------------------- */
@ -144,19 +141,18 @@ void ComputeGrid::deallocate()
memory->destroy(gridall);
if (gridlocal_allocated) {
gridlocal_allocated = 0;
memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo);
memory->destroy4d_offset(gridlocal, nzlo, nylo, nxlo);
}
array = nullptr;
}
/* ----------------------------------------------------------------------
set global grid
------------------------------------------------------------------------- */
void ComputeGrid::set_grid_global()
{
// calculate grid layout
// calculate grid layout
triclinic = domain->triclinic;
@ -176,13 +172,13 @@ void ComputeGrid::set_grid_global()
double yprd = prd[1];
double zprd = prd[2];
delxinv = nx/xprd;
delyinv = ny/yprd;
delzinv = nz/zprd;
delxinv = nx / xprd;
delyinv = ny / yprd;
delzinv = nz / zprd;
delx = 1.0/delxinv;
dely = 1.0/delyinv;
delz = 1.0/delzinv;
delx = 1.0 / delxinv;
dely = 1.0 / delyinv;
delz = 1.0 / delzinv;
}
/* ----------------------------------------------------------------------
@ -202,15 +198,15 @@ void ComputeGrid::set_grid_local()
// the 2 equality if tests insure a consistent decision
// as to which proc owns it
double xfraclo,xfrachi,yfraclo,yfrachi,zfraclo,zfrachi;
double xfraclo, xfrachi, yfraclo, yfrachi, zfraclo, zfrachi;
if (comm->layout != Comm::LAYOUT_TILED) {
xfraclo = comm->xsplit[comm->myloc[0]];
xfrachi = comm->xsplit[comm->myloc[0]+1];
xfrachi = comm->xsplit[comm->myloc[0] + 1];
yfraclo = comm->ysplit[comm->myloc[1]];
yfrachi = comm->ysplit[comm->myloc[1]+1];
yfrachi = comm->ysplit[comm->myloc[1] + 1];
zfraclo = comm->zsplit[comm->myloc[2]];
zfrachi = comm->zsplit[comm->myloc[2]+1];
zfrachi = comm->zsplit[comm->myloc[2] + 1];
} else {
xfraclo = comm->mysplit[0][0];
xfrachi = comm->mysplit[0][1];
@ -220,23 +216,22 @@ void ComputeGrid::set_grid_local()
zfrachi = comm->mysplit[2][1];
}
nxlo = static_cast<int> (xfraclo * nx);
if (1.0*nxlo != xfraclo*nx) nxlo++;
nxhi = static_cast<int> (xfrachi * nx);
if (1.0*nxhi == xfrachi*nx) nxhi--;
nxlo = static_cast<int>(xfraclo * nx);
if (1.0 * nxlo != xfraclo * nx) nxlo++;
nxhi = static_cast<int>(xfrachi * nx);
if (1.0 * nxhi == xfrachi * nx) nxhi--;
nylo = static_cast<int> (yfraclo * ny);
if (1.0*nylo != yfraclo*ny) nylo++;
nyhi = static_cast<int> (yfrachi * ny);
if (1.0*nyhi == yfrachi*ny) nyhi--;
nylo = static_cast<int>(yfraclo * ny);
if (1.0 * nylo != yfraclo * ny) nylo++;
nyhi = static_cast<int>(yfrachi * ny);
if (1.0 * nyhi == yfrachi * ny) nyhi--;
nzlo = static_cast<int> (zfraclo * nz);
if (1.0*nzlo != zfraclo*nz) nzlo++;
nzhi = static_cast<int> (zfrachi * nz);
if (1.0*nzhi == zfrachi*nz) nzhi--;
nzlo = static_cast<int>(zfraclo * nz);
if (1.0 * nzlo != zfraclo * nz) nzlo++;
nzhi = static_cast<int>(zfrachi * nz);
if (1.0 * nzhi == zfrachi * nz) nzhi--;
ngridlocal = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1);
}
/* ----------------------------------------------------------------------
@ -245,10 +240,8 @@ void ComputeGrid::set_grid_local()
double ComputeGrid::memory_usage()
{
double nbytes = size_array_rows*size_array_cols *
sizeof(double); // grid
nbytes += size_array_rows*size_array_cols *
sizeof(double); // gridall
nbytes += size_array_cols*ngridlocal*sizeof(double); // gridlocal
double nbytes = size_array_rows * size_array_cols * sizeof(double); // grid
nbytes += size_array_rows * size_array_cols * sizeof(double); // gridall
nbytes += size_array_cols * ngridlocal * sizeof(double); // gridlocal
return nbytes;
}

47
src/compute_grid.h
View File

@ -20,7 +20,6 @@ namespace LAMMPS_NS {
class ComputeGrid : public Compute {
public:
ComputeGrid(class LAMMPS *, int, char **);
~ComputeGrid() override;
void init() override;
@ -30,31 +29,31 @@ class ComputeGrid : public Compute {
double memory_usage() override;
protected:
int nx, ny, nz; // global grid dimensions
int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive
int ngridlocal; // number of local grid points
int nvalues; // number of values per grid point
double **grid; // global grid
double **gridall; // global grid summed over procs
double ****gridlocal; // local grid
int triclinic; // triclinic flag
double *boxlo, *prd; // box info (units real/ortho or reduced/tri)
double *sublo, *subhi; // subdomain info (units real/ortho or reduced/tri)
double delxinv,delyinv,delzinv; // inverse grid spacing
double delx,dely,delz; // grid spacing
int nargbase; // number of base class args
double cutmax; // largest cutoff distance
int size_array_cols_base; // number of columns used for coords, etc.
int gridlocal_allocated; // shows if gridlocal allocated
int nx, ny, nz; // global grid dimensions
int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive
int ngridlocal; // number of local grid points
int nvalues; // number of values per grid point
double **grid; // global grid
double **gridall; // global grid summed over procs
double ****gridlocal; // local grid
int triclinic; // triclinic flag
double *boxlo, *prd; // box info (units real/ortho or reduced/tri)
double *sublo, *subhi; // subdomain info (units real/ortho or reduced/tri)
double delxinv, delyinv, delzinv; // inverse grid spacing
double delx, dely, delz; // grid spacing
int nargbase; // number of base class args
double cutmax; // largest cutoff distance
int size_array_cols_base; // number of columns used for coords, etc.
int gridlocal_allocated; // shows if gridlocal allocated
void allocate(); // create arrays
void deallocate(); // free arrays
void grid2x(int, double*); // convert grid point to coord
void assign_coords_all(); // assign coords for global grid
void set_grid_global(); // set global grid
void set_grid_local(); // set bounds for local grid
void allocate(); // create arrays
void deallocate(); // free arrays
void grid2x(int, double *); // convert grid point to coord
void assign_coords_all(); // assign coords for global grid
void set_grid_global(); // set global grid
void set_grid_local(); // set bounds for local grid
};
}
} // namespace LAMMPS_NS
#endif

107
src/compute_grid_local.cpp
View File

@ -12,29 +12,30 @@
------------------------------------------------------------------------- */
#include "compute_grid_local.h"
#include <mpi.h>
#include <cstring>
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "comm.h"
#include "modify.h"
#include "update.h"
#include <cstring>
// For the subdomain test below; grid-points and subdomain boundaries
// sometimes differ by minimal amounts (in the order of 2e-17).
#define EPSILON 1.0e-10
static constexpr double EPSILON = 1.0e-10;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeGridLocal::ComputeGridLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), alocal(nullptr)
Compute(lmp, narg, arg), alocal(nullptr)
{
if (narg < 6) error->all(FLERR,"Illegal compute grid/local command");
if (narg < 6) error->all(FLERR, "Illegal compute grid/local command");
local_flag = 1;
size_local_cols = 0;
@ -43,15 +44,16 @@ ComputeGridLocal::ComputeGridLocal(LAMMPS *lmp, int narg, char **arg) :
int iarg0 = 3;
int iarg = iarg0;
if (strcmp(arg[iarg],"grid") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal compute grid/local command");
nx = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
ny = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
nz = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
if (strcmp(arg[iarg], "grid") == 0) {
if (iarg + 4 > narg) error->all(FLERR, "Illegal compute grid/local command");
nx = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
ny = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
nz = utils::inumeric(FLERR, arg[iarg + 3], false, lmp);
if (nx <= 0 || ny <= 0 || nz <= 0)
error->all(FLERR,"All grid/local dimensions must be positive");
error->all(FLERR, "All grid/local dimensions must be positive");
iarg += 4;
} else error->all(FLERR,"Illegal compute grid/local command");
} else
error->all(FLERR, "Illegal compute grid/local command");
nargbase = iarg - iarg0;
@ -68,9 +70,7 @@ ComputeGridLocal::~ComputeGridLocal()
/* ---------------------------------------------------------------------- */
void ComputeGridLocal::init()
{
}
void ComputeGridLocal::init() {}
/* ---------------------------------------------------------------------- */
@ -89,9 +89,9 @@ void ComputeGridLocal::setup()
void ComputeGridLocal::grid2x(int ix, int iy, int iz, double *x)
{
x[0] = ix*delx;
x[1] = iy*dely;
x[2] = iz*delz;
x[0] = ix * delx;
x[1] = iy * dely;
x[2] = iz * delz;
if (triclinic) domain->lamda2x(x, x);
}
@ -103,12 +103,11 @@ void ComputeGridLocal::grid2x(int ix, int iy, int iz, double *x)
void ComputeGridLocal::grid2lamda(int ix, int iy, int iz, double *x)
{
x[0] = ix*delx;
x[1] = iy*dely;
x[2] = iz*delz;
x[0] = ix * delx;
x[1] = iy * dely;
x[2] = iz * delz;
}
/* ----------------------------------------------------------------------
create arrays
------------------------------------------------------------------------- */
@ -141,7 +140,7 @@ void ComputeGridLocal::deallocate()
void ComputeGridLocal::set_grid_global()
{
// calculate grid layout
// calculate grid layout
triclinic = domain->triclinic;
@ -161,13 +160,13 @@ void ComputeGridLocal::set_grid_global()
double yprd = prd[1];
double zprd = prd[2];
delxinv = nx/xprd;
delyinv = ny/yprd;
delzinv = nz/zprd;
delxinv = nx / xprd;
delyinv = ny / yprd;
delzinv = nz / zprd;
delx = 1.0/delxinv;
dely = 1.0/delyinv;
delz = 1.0/delzinv;
delx = 1.0 / delxinv;
dely = 1.0 / delyinv;
delz = 1.0 / delzinv;
}
/* ----------------------------------------------------------------------
@ -187,15 +186,15 @@ void ComputeGridLocal::set_grid_local()
// the 2 equality if tests insure a consistent decision
// as to which proc owns it
double xfraclo,xfrachi,yfraclo,yfrachi,zfraclo,zfrachi;
double xfraclo, xfrachi, yfraclo, yfrachi, zfraclo, zfrachi;
if (comm->layout != Comm::LAYOUT_TILED) {
xfraclo = comm->xsplit[comm->myloc[0]];
xfrachi = comm->xsplit[comm->myloc[0]+1];
xfrachi = comm->xsplit[comm->myloc[0] + 1];
yfraclo = comm->ysplit[comm->myloc[1]];
yfrachi = comm->ysplit[comm->myloc[1]+1];
yfrachi = comm->ysplit[comm->myloc[1] + 1];
zfraclo = comm->zsplit[comm->myloc[2]];
zfrachi = comm->zsplit[comm->myloc[2]+1];
zfrachi = comm->zsplit[comm->myloc[2] + 1];
} else {
xfraclo = comm->mysplit[0][0];
xfrachi = comm->mysplit[0][1];
@ -205,20 +204,20 @@ void ComputeGridLocal::set_grid_local()
zfrachi = comm->mysplit[2][1];
}
nxlo = static_cast<int> (xfraclo * nx);
if (1.0*nxlo != xfraclo*nx) nxlo++;
nxhi = static_cast<int> (xfrachi * nx);
if (1.0*nxhi == xfrachi*nx) nxhi--;
nxlo = static_cast<int>(xfraclo * nx);
if (1.0 * nxlo != xfraclo * nx) nxlo++;
nxhi = static_cast<int>(xfrachi * nx);
if (1.0 * nxhi == xfrachi * nx) nxhi--;
nylo = static_cast<int> (yfraclo * ny);
if (1.0*nylo != yfraclo*ny) nylo++;
nyhi = static_cast<int> (yfrachi * ny);
if (1.0*nyhi == yfrachi*ny) nyhi--;
nylo = static_cast<int>(yfraclo * ny);
if (1.0 * nylo != yfraclo * ny) nylo++;
nyhi = static_cast<int>(yfrachi * ny);
if (1.0 * nyhi == yfrachi * ny) nyhi--;
nzlo = static_cast<int> (zfraclo * nz);
if (1.0*nzlo != zfraclo*nz) nzlo++;
nzhi = static_cast<int> (zfrachi * nz);
if (1.0*nzhi == zfrachi*nz) nzhi--;
nzlo = static_cast<int>(zfraclo * nz);
if (1.0 * nzlo != zfraclo * nz) nzlo++;
nzhi = static_cast<int>(zfrachi * nz);
if (1.0 * nzhi == zfrachi * nz) nzhi--;
size_local_rows = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1);
}
@ -248,10 +247,10 @@ void ComputeGridLocal::assign_coords()
// ensure gridpoint is not strictly outside subdomain
if ((sublo[0]-xgrid[0]) > EPSILON || (xgrid[0]-subhi[0]) > EPSILON ||
(sublo[1]-xgrid[1]) > EPSILON || (xgrid[1]-subhi[1]) > EPSILON ||
(sublo[2]-xgrid[2]) > EPSILON || (xgrid[2]-subhi[2]) > EPSILON)
error->one(FLERR,"Invalid gridpoint position in compute grid/local");
if ((sublo[0] - xgrid[0]) > EPSILON || (xgrid[0] - subhi[0]) > EPSILON ||
(sublo[1] - xgrid[1]) > EPSILON || (xgrid[1] - subhi[1]) > EPSILON ||
(sublo[2] - xgrid[2]) > EPSILON || (xgrid[2] - subhi[2]) > EPSILON)
error->one(FLERR, "Invalid gridpoint position in compute grid/local");
// convert lamda to x, y, z, after sudomain check
@ -270,6 +269,6 @@ void ComputeGridLocal::assign_coords()
double ComputeGridLocal::memory_usage()
{
int nbytes = size_local_rows * size_local_cols * sizeof(double); // gridlocal
int nbytes = size_local_rows * size_local_cols * sizeof(double); // gridlocal
return nbytes;
}

43
src/compute_grid_local.h
View File

@ -20,7 +20,6 @@ namespace LAMMPS_NS {
class ComputeGridLocal : public Compute {
public:
ComputeGridLocal(class LAMMPS *, int, char **);
~ComputeGridLocal() override;
void init() override;
@ -30,29 +29,29 @@ class ComputeGridLocal : public Compute {
double memory_usage() override;
protected:
int nx, ny, nz; // global grid dimensions
int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive
int nvalues; // number of values per grid point
double **alocal; // pointer to Compute::array_local
int triclinic; // triclinic flag
double *boxlo, *prd; // box info (units real/ortho or reduced/tri)
double *sublo, *subhi; // subdomain info (units real/ortho or reduced/tri)
double delxinv,delyinv,delzinv; // inverse grid spacing
double delx,dely,delz; // grid spacing
int nargbase; // number of base class args
double cutmax; // largest cutoff distance
int size_local_cols_base; // number of columns used for coords, etc.
int gridlocal_allocated; // shows if gridlocal allocated
int nx, ny, nz; // global grid dimensions
int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive
int nvalues; // number of values per grid point
double **alocal; // pointer to Compute::array_local
int triclinic; // triclinic flag
double *boxlo, *prd; // box info (units real/ortho or reduced/tri)
double *sublo, *subhi; // subdomain info (units real/ortho or reduced/tri)
double delxinv, delyinv, delzinv; // inverse grid spacing
double delx, dely, delz; // grid spacing
int nargbase; // number of base class args
double cutmax; // largest cutoff distance
int size_local_cols_base; // number of columns used for coords, etc.
int gridlocal_allocated; // shows if gridlocal allocated
void allocate(); // create arrays
void deallocate(); // free arrays
void grid2x(int, int, int, double*); // convert global indices to coordinates
void grid2lamda(int, int, int, double*); // convert global indices to lamda coordinates
void set_grid_global(); // set global grid
void set_grid_local(); // set bounds for local grid
void assign_coords(); // assign coords for grid
void allocate(); // create arrays
void deallocate(); // free arrays
void grid2x(int, int, int, double *); // convert global indices to coordinates
void grid2lamda(int, int, int, double *); // convert global indices to lamda coordinates
void set_grid_global(); // set global grid
void set_grid_local(); // set bounds for local grid
void assign_coords(); // assign coords for grid
};
}
} // namespace LAMMPS_NS
#endif