From ea9c1002fed3261274eeac0769b165a2f27fcea2 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 11 Oct 2019 17:51:19 -0600 Subject: [PATCH 01/19] Created placeholders for ComputeGrid and ComputeSNAGrid classes --- src/compute_grid.cpp | 186 ++++++++++++++++++++++++++ src/compute_grid.h | 69 ++++++++++ src/compute_sna_grid.cpp | 277 +++++++++++++++++++++++++++++++++++++++ src/compute_sna_grid.h | 75 +++++++++++ 4 files changed, 607 insertions(+) create mode 100644 src/compute_grid.cpp create mode 100644 src/compute_grid.h create mode 100644 src/compute_sna_grid.cpp create mode 100644 src/compute_sna_grid.h diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp new file mode 100644 index 0000000000..c67d72160b --- /dev/null +++ b/src/compute_grid.cpp @@ -0,0 +1,186 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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_grid.h" +#include +#include +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{ONCE,NFREQ,EVERY}; + +/* ---------------------------------------------------------------------- */ + +ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), + idchunk(NULL), masstotal(NULL), massproc(NULL), com(NULL), comall(NULL) +{ + if (narg != 4) error->all(FLERR,"Illegal compute com/chunk command"); + + array_flag = 1; + size_array_cols = 3; + size_array_rows = 0; + size_array_rows_variable = 1; + extarray = 0; + + // ID of compute chunk/atom + + int n = strlen(arg[3]) + 1; + idchunk = new char[n]; + strcpy(idchunk,arg[3]); + + init(); + + // chunk-based data + + nchunk = 1; + maxchunk = 0; + allocate(); + + firstflag = massneed = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeGrid::~ComputeGrid() +{ + delete [] idchunk; + memory->destroy(massproc); + memory->destroy(masstotal); + memory->destroy(com); + memory->destroy(comall); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGrid::init() +{ + int icompute = modify->find_compute(idchunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for compute com/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Compute com/chunk does not use chunk/atom compute"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGrid::setup() +{ + // one-time calculation of per-chunk mass + // done in setup, so that ComputeChunkAtom::setup() is already called + + if (firstflag && cchunk->idsflag == ONCE) { + compute_array(); + firstflag = massneed = 0; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGrid::compute_array() +{ + int index; + double massone; + double unwrap[3]; + + invoked_array = update->ntimestep; + + // compute chunk/atom assigns atoms to chunk IDs + // extract ichunk index vector from compute + // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms + + nchunk = cchunk->setup_chunks(); + cchunk->compute_ichunk(); + int *ichunk = cchunk->ichunk; + + if (nchunk > maxchunk) allocate(); + size_array_rows = nchunk; + + // zero local per-chunk values + + for (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; + + // compute COM for each chunk + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + imageint *image = atom->image; + double *mass = atom->mass; + double *rmass = atom->rmass; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + 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); + 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); + + 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; + } +} + +/* ---------------------------------------------------------------------- + free and reallocate per-chunk arrays +------------------------------------------------------------------------- */ + +void ComputeGrid::allocate() +{ + memory->destroy(massproc); + memory->destroy(masstotal); + 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"); + array = comall; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeGrid::memory_usage() +{ + double bytes = (bigint) maxchunk * 2 * sizeof(double); + bytes += (bigint) maxchunk * 2*3 * sizeof(double); + return bytes; +} diff --git a/src/compute_grid.h b/src/compute_grid.h new file mode 100644 index 0000000000..cece7050de --- /dev/null +++ b/src/compute_grid.h @@ -0,0 +1,69 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(grid,ComputeGrid) + +#else + +#ifndef LMP_COMPUTE_GRID_H +#define LMP_COMPUTE_GRID_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeGrid : public Compute { + public: + char *idchunk; // fields accessed by other classes + double *masstotal; + + ComputeGrid(class LAMMPS *, int, char **); + ~ComputeGrid(); + void init(); + void setup(); + void compute_array(); + + void lock_enable(); + void lock_disable(); + int lock_length(); + void lock(class Fix *, bigint, bigint); + void unlock(class Fix *); + + double memory_usage(); + + private: + int nchunk,maxchunk; + int firstflag,massneed; + + double *massproc; + double **com,**comall; + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +*/ diff --git a/src/compute_sna_grid.cpp b/src/compute_sna_grid.cpp new file mode 100644 index 0000000000..58f7e9fc38 --- /dev/null +++ b/src/compute_sna_grid.cpp @@ -0,0 +1,277 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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_sna_grid.h" +#include +#include +#include "sna.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "force.h" +#include "pair.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), cutsq(NULL), list(NULL), sna(NULL), + radelem(NULL), wjelem(NULL) +{ + double rmin0, rfac0; + int twojmax, switchflag, bzeroflag; + radelem = NULL; + wjelem = NULL; + + int ntypes = atom->ntypes; + int nargmin = 6+2*ntypes; + + if (narg < nargmin) error->all(FLERR,"Illegal compute sna/grid command"); + + // default values + + rmin0 = 0.0; + switchflag = 1; + bzeroflag = 1; + quadraticflag = 0; + + // offset by 1 to match up with types + + memory->create(radelem,ntypes+1,"sna/grid:radelem"); + 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]); + + // 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; + 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; + } + } + + // process optional args + + 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]); + 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]); + 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]); + 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]); + iarg += 2; + } else error->all(FLERR,"Illegal compute sna/grid command"); + } + + snaptr = new SNA(lmp,rfac0,twojmax, + rmin0,switchflag,bzeroflag); + + ncoeff = snaptr->ncoeff; + size_peratom_cols = ncoeff; + if (quadraticflag) size_peratom_cols += (ncoeff*(ncoeff+1))/2; + peratom_flag = 1; + + nmax = 0; + sna = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeSNAGrid::~ComputeSNAGrid() +{ + memory->destroy(sna); + memory->destroy(radelem); + memory->destroy(wjelem); + memory->destroy(cutsq); + delete snaptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGrid::init() +{ + if (force->pair == NULL) + error->all(FLERR,"Compute sna/grid requires a pair style be defined"); + + if (cutmax > force->pair->cutforce) + error->all(FLERR,"Compute sna/grid cutoff is longer than pairwise cutoff"); + + // need an occasional full neighbor list + + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + neighbor->requests[irequest]->occasional = 1; + + 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"); + snaptr->init(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGrid::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGrid::compute_pergrid() +{ + invoked_peratom = update->ntimestep; + + // grow sna array if necessary + + if (atom->nmax > nmax) { + memory->destroy(sna); + nmax = atom->nmax; + memory->create(sna,nmax,size_peratom_cols,"sna/grid:sna"); + array_atom = sna; + } + + // invoke full neighbor list (will copy or build if necessary) + + neighbor->build_one(list); + + const int inum = list->inum; + const int* const ilist = list->ilist; + const int* const numneigh = list->numneigh; + int** const firstneigh = list->firstneigh; + int * const type = atom->type; + + // compute sna for each atom in group + // use full neighbor list to count atoms less than cutoff + + double** const x = atom->x; + const int* const mask = atom->mask; + + for (int ii = 0; ii < inum; ii++) { + const int i = ilist[ii]; + if (mask[i] & groupbit) { + + const double xtmp = x[i][0]; + const double ytmp = x[i][1]; + const double ztmp = x[i][2]; + const int itype = type[i]; + const double radi = radelem[itype]; + const int* const jlist = firstneigh[i]; + const int jnum = numneigh[i]; + + // insure rij, inside, and typej are of size jnum + + snaptr->grow_rij(jnum); + + // rij[][3] = displacements between atom I and those neighbors + // inside = indices of neighbors of I within cutoff + // typej = types of neighbors of I within cutoff + + int ninside = 0; + for (int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + j &= NEIGHMASK; + + 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; + int jtype = type[j]; + if (rsq < cutsq[itype][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] = (radi+radelem[jtype])*rcutfac; + ninside++; + } + } + + snaptr->compute_ui(ninside); + snaptr->compute_zi(); + snaptr->compute_bi(); + for (int icoeff = 0; icoeff < ncoeff; icoeff++) + sna[i][icoeff] = snaptr->blist[icoeff]; + if (quadraticflag) { + int ncount = ncoeff; + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + double bi = snaptr->blist[icoeff]; + + // diagonal element of quadratic matrix + + sna[i][ncount++] = 0.5*bi*bi; + + // upper-triangular elements of quadratic matrix + + for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) + sna[i][ncount++] = bi*snaptr->blist[jcoeff]; + } + } + } else { + for (int icoeff = 0; icoeff < size_peratom_cols; icoeff++) + sna[i][icoeff] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + memory usage +------------------------------------------------------------------------- */ + +double ComputeSNAGrid::memory_usage() +{ + double bytes = nmax*size_peratom_cols * sizeof(double); // sna + bytes += snaptr->memory_usage(); // SNA object + + return bytes; +} + diff --git a/src/compute_sna_grid.h b/src/compute_sna_grid.h new file mode 100644 index 0000000000..f85d9679cc --- /dev/null +++ b/src/compute_sna_grid.h @@ -0,0 +1,75 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(sna/grid,ComputeSNAGrid) + +#else + +#ifndef LMP_COMPUTE_SNA_GRID_H +#define LMP_COMPUTE_SNA_GRID_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeSNAGrid : public Compute { + public: + ComputeSNAGrid(class LAMMPS *, int, char **); + ~ComputeSNAGrid(); + void init(); + void init_list(int, class NeighList *); + void compute_grid(); + double memory_usage(); + + private: + int nmax; + int ncoeff; + double **cutsq; + class NeighList *list; + double **sna; + double rcutfac; + double *radelem; + double *wjelem; + class SNA* snaptr; + double cutmax; + int quadraticflag; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Compute sna/grid requires a pair style be defined + +Self-explanatory. + +E: Compute sna/grid cutoff is longer than pairwise cutoff + +Self-explanatory. + +W: More than one compute sna/grid + +Self-explanatory. + +*/ From 762ecf8f0e5bb8fb146bb7ef42bd80e9a4e7f9e1 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Sat, 19 Oct 2019 17:03:19 -0600 Subject: [PATCH 02/19] Completed serial version with PBC, but incorrect --- src/compute_grid.cpp | 255 +++++++++++++++++++++++---------------- src/compute_grid.h | 45 ++++--- src/compute_sna_grid.cpp | 166 +++++++++++++------------ src/compute_sna_grid.h | 8 +- 4 files changed, 266 insertions(+), 208 deletions(-) diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp index c67d72160b..d8c7eb3697 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -18,169 +18,218 @@ #include "update.h" #include "modify.h" #include "domain.h" +#include "force.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; -enum{ONCE,NFREQ,EVERY}; - /* ---------------------------------------------------------------------- */ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - idchunk(NULL), masstotal(NULL), massproc(NULL), com(NULL), comall(NULL) + Compute(lmp, narg, arg) { - if (narg != 4) error->all(FLERR,"Illegal compute com/chunk command"); + if (narg < 6) error->all(FLERR,"Illegal compute grid command"); array_flag = 1; - size_array_cols = 3; + size_array_cols = 0; size_array_rows = 0; - size_array_rows_variable = 1; - extarray = 0; + extarray = 1; - // ID of compute chunk/atom + int iarg0 = 3; + int iarg = iarg0; + if (strcmp(arg[iarg],"grid") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal compute grid command"); + nx = force->inumeric(FLERR,arg[iarg+1]); + ny = force->inumeric(FLERR,arg[iarg+2]); + nz = force->inumeric(FLERR,arg[iarg+3]); + 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"); - int n = strlen(arg[3]) + 1; - idchunk = new char[n]; - strcpy(idchunk,arg[3]); + nargbase = iarg - iarg0; - init(); - - // chunk-based data - - nchunk = 1; - maxchunk = 0; - allocate(); - - firstflag = massneed = 1; + size_array_rows = nx*ny*nz; } /* ---------------------------------------------------------------------- */ ComputeGrid::~ComputeGrid() { - delete [] idchunk; - memory->destroy(massproc); - memory->destroy(masstotal); - memory->destroy(com); - memory->destroy(comall); } /* ---------------------------------------------------------------------- */ void ComputeGrid::init() { - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - error->all(FLERR,"Chunk/atom compute does not exist for compute com/chunk"); - cchunk = (ComputeChunkAtom *) modify->compute[icompute]; - if (strcmp(cchunk->style,"chunk/atom") != 0) - error->all(FLERR,"Compute com/chunk does not use chunk/atom compute"); } /* ---------------------------------------------------------------------- */ void ComputeGrid::setup() { - // one-time calculation of per-chunk mass - // done in setup, so that ComputeChunkAtom::setup() is already called + + // calculate grid layout - if (firstflag && cchunk->idsflag == ONCE) { - compute_array(); - firstflag = massneed = 0; + triclinic = domain->triclinic; + + if (triclinic == 0) { + prd = domain->prd; + boxlo = domain->boxlo; + } else { + prd = domain->prd_lamda; + boxlo = domain->boxlo_lamda; } + + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; + + delxinv = nx/xprd; + delyinv = ny/yprd; + delzinv = nz/zprd; + + delx = 1.0/delxinv; + dely = 1.0/delyinv; + delz = 1.0/delzinv; + + // sufficient conditions for stencil bounding rcut + + // require |delz*mz|^2 <= rcut^2 + // require |dely*my|^2 <= rcut^2 + |delyz*mz_max|^2 + // require |delx*mx|^2 <= rcut^2 + |delxz*mz_max|^2 + |delxy*my_max|^2 + + double delxy = domain->xy/ny; + double delxz = domain->xz/nz; + double delyz = domain->yz/nz; + + if (!triclinic) { + mz = cutmax*delzinv + 1; + my = sqrt(cutmax*cutmax + pow(delyz*mz,2))*delyinv + 1; + mx = sqrt(cutmax*cutmax + pow(delxz*mz,2) + + pow(delxy*my,2))*delxinv + 1; + } else { + double delxinvtmp = nx/domain->xprd; + double delyinvtmp = ny/domain->yprd; + double delzinvtmp = nz/domain->zprd; + mz = cutmax*delzinvtmp + 1; + my = sqrt(cutmax*cutmax + pow(delyz*mz,2))*delyinvtmp + 1; + mx = sqrt(cutmax*cutmax + pow(delxz*mz,2) + + pow(delxy*my,2))*delxinvtmp + 1; + } + + printf("mx = %d\n",mx); + printf("my = %d\n",my); + printf("mz = %d\n",mz); + + // size global grid to accomodate periodic interactions + + nxfull = nx + 2*mx; + nyfull = ny + 2*my; + nzfull = nz + 2*mz; + nxyfull = nxfull * nyfull; + + printf("nxfull = %d\n",nxfull); + printf("nyfull = %d\n",nyfull); + printf("nzfull = %d\n",nzfull); + + x0full = boxlo[0] - mx*delx; + y0full = boxlo[1] - my*dely; + z0full = boxlo[2] - mz*delz; + + allocate(); } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + convert grid index to box coords +------------------------------------------------------------------------- */ -void ComputeGrid::compute_array() +void ComputeGrid::igridfull2x(int igrid, double *x) { - int index; - double massone; - double unwrap[3]; + int iz = igrid / nxyfull; + igrid -= iz*nxyfull; + int iy = igrid / nxfull; + igrid -= igrid*nxfull; + int ix = igrid; - invoked_array = update->ntimestep; + x[0] = x0full+ix*delx; + x[1] = y0full+iy*dely; + x[2] = z0full+iz*delz; - // 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 + if (triclinic) domain->lamda2x(x, x); - nchunk = cchunk->setup_chunks(); - cchunk->compute_ichunk(); - int *ichunk = cchunk->ichunk; +} - if (nchunk > maxchunk) allocate(); - size_array_rows = nchunk; +/* ---------------------------------------------------------------------- + gather global array from full grid +------------------------------------------------------------------------- */ - // zero local per-chunk values +void ComputeGrid::gather_global_array() +{ + int iarray; + memset(&array[0][0],0,size_array_rows*size_array_cols*sizeof(double)); - 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; + for (int igrid = 0; igrid < ngridfull; igrid++) { - // compute COM for each chunk + // inefficient, should exploit shared ix structure - double **x = atom->x; - int *mask = atom->mask; - int *type = atom->type; - imageint *image = atom->image; - double *mass = atom->mass; - double *rmass = atom->rmass; - int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - 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); - 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); - - 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; + iarray = igridfull2iarray(igrid); + for (int icol = 0; icol < size_array_cols; icol++) + array[iarray][icol] += gridfull[igrid][icol]; } } /* ---------------------------------------------------------------------- - free and reallocate per-chunk arrays + convert full grid index to compute array index + inefficient, should exploit shared ix structure +------------------------------------------------------------------------- */ + +int ComputeGrid::igridfull2iarray(int igrid) +{ + int iz = igrid / nxyfull; + igrid -= iz*nxyfull; + int iy = igrid / nxfull; + igrid -= igrid*nxfull; + int ix = igrid; + + ix -= mx; + iy -= my; + iz -= mz; + + while (ix < 0) ix += nx; + while (iy < 0) iy += ny; + while (iz < 0) iz += nz; + + while (ix >= nx) ix -= nx; + while (iy >= ny) iy -= ny; + while (iz >= nz) iz -= nz; + + int iarray = (iz * ny + iy) * nx + ix; + + return iarray; +} + +/* ---------------------------------------------------------------------- + free and reallocate arrays ------------------------------------------------------------------------- */ void ComputeGrid::allocate() { - memory->destroy(massproc); - memory->destroy(masstotal); - 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"); - array = comall; -} + ngridfull = nxfull*nyfull*nzfull; + // grow global array if necessary + + memory->destroy(array); + memory->create(array,size_array_rows,size_array_cols,"sna/grid:array"); + memory->create(gridfull,ngridfull,size_array_cols,"sna/grid:gridfull"); +} /* ---------------------------------------------------------------------- memory usage of local data ------------------------------------------------------------------------- */ double ComputeGrid::memory_usage() { - double bytes = (bigint) maxchunk * 2 * sizeof(double); - bytes += (bigint) maxchunk * 2*3 * sizeof(double); - return bytes; + int nbytes = 0; + return nbytes; } diff --git a/src/compute_grid.h b/src/compute_grid.h index cece7050de..61775186ff 100644 --- a/src/compute_grid.h +++ b/src/compute_grid.h @@ -11,12 +11,6 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#ifdef COMPUTE_CLASS - -ComputeStyle(grid,ComputeGrid) - -#else - #ifndef LMP_COMPUTE_GRID_H #define LMP_COMPUTE_GRID_H @@ -26,36 +20,39 @@ namespace LAMMPS_NS { class ComputeGrid : public Compute { public: - char *idchunk; // fields accessed by other classes - double *masstotal; ComputeGrid(class LAMMPS *, int, char **); - ~ComputeGrid(); + virtual ~ComputeGrid(); void init(); void setup(); - void compute_array(); - - void lock_enable(); - void lock_disable(); - int lock_length(); - void lock(class Fix *, bigint, bigint); - void unlock(class Fix *); + virtual void compute_array() = 0; double memory_usage(); + protected: + int nx, ny, nz; // grid dimensions + int nxfull, nyfull, nzfull; // grid dimensions with ghost points + int nxyfull; // nx_full*ny_full + int ngridfull; // number of full grid points + double **gridfull; // full grid points + int mx, my, mz; // cutmax stencil dimensions + int triclinic; // triclinic flag + double *boxlo, *prd; // box info (units real/ortho or reduced/tri) + double delxinv,delyinv,delzinv; // inverse grid spacing + double delx,dely,delz; // grid spacing + double x0full, y0full, z0full; // origin of full grid + int nargbase; // number of base class args + double cutmax; // largest cutoff distance + virtual void allocate(); + void igridfull2x(int, double*); // convert full grid point to coord + void gather_global_array(); // gather global array from full grid + int igridfull2iarray(int); // convert full grid index to compute array index + private: - int nchunk,maxchunk; - int firstflag,massneed; - - double *massproc; - double **com,**comall; - - void allocate(); }; } -#endif #endif /* ERROR/WARNING messages: diff --git a/src/compute_sna_grid.cpp b/src/compute_sna_grid.cpp index 58f7e9fc38..1f998ced9b 100644 --- a/src/compute_sna_grid.cpp +++ b/src/compute_sna_grid.cpp @@ -11,6 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +#include "compute_grid.h" #include "compute_sna_grid.h" #include #include @@ -30,7 +31,7 @@ using namespace LAMMPS_NS; ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), cutsq(NULL), list(NULL), sna(NULL), + ComputeGrid(lmp, narg, arg), cutsq(NULL), list(NULL), sna(NULL), radelem(NULL), wjelem(NULL) { double rmin0, rfac0; @@ -38,6 +39,13 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : radelem = NULL; wjelem = NULL; + // skip over arguments used by base class + // so that argument positions are identical to + // regular per-atom compute + + arg += nargbase; + narg -= nargbase; + int ntypes = atom->ntypes; int nargmin = 6+2*ntypes; @@ -58,6 +66,7 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : rcutfac = atof(arg[3]); rfac0 = atof(arg[4]); twojmax = atoi(arg[5]); + printf("rcutfac = %g rfac0 = %g twojmax = %d\n",rcutfac, rfac0, twojmax); for(int i = 0; i < ntypes; i++) radelem[i+1] = atof(arg[6+i]); @@ -73,6 +82,7 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : cut = 2.0*radelem[i]*rcutfac; if (cut > cutmax) cutmax = cut; cutsq[i][i] = cut*cut; + printf("i = %d cutsq[i][i] = %g\n",i,cutsq[i][i]); for(int j = i+1; j <= ntypes; j++) { cut = (radelem[i]+radelem[j])*rcutfac; cutsq[i][j] = cutsq[j][i] = cut*cut; @@ -84,6 +94,7 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : int iarg = nargmin; while (iarg < narg) { + printf("iarg = %d arg = %s\n",iarg, arg[iarg]); if (strcmp(arg[iarg],"rmin0") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute sna/grid command"); @@ -105,18 +116,19 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : quadraticflag = atoi(arg[iarg+1]); iarg += 2; } else error->all(FLERR,"Illegal compute sna/grid command"); + } + printf("rmin0 = %g, bzeroflag = %d, quadraticflag = %d\n", + rmin0, bzeroflag, quadraticflag); + snaptr = new SNA(lmp,rfac0,twojmax, rmin0,switchflag,bzeroflag); ncoeff = snaptr->ncoeff; - size_peratom_cols = ncoeff; - if (quadraticflag) size_peratom_cols += (ncoeff*(ncoeff+1))/2; - peratom_flag = 1; - - nmax = 0; - sna = NULL; + size_array_cols = ncoeff; + if (quadraticflag) size_array_cols += (ncoeff*(ncoeff+1))/2; + array_flag = 1; } /* ---------------------------------------------------------------------- */ @@ -166,101 +178,101 @@ void ComputeSNAGrid::init_list(int /*id*/, NeighList *ptr) /* ---------------------------------------------------------------------- */ -void ComputeSNAGrid::compute_pergrid() +void ComputeSNAGrid::compute_array() { - invoked_peratom = update->ntimestep; + invoked_array = update->ntimestep; - // grow sna array if necessary +// // invoke full neighbor list (will copy or build if necessary) - if (atom->nmax > nmax) { - memory->destroy(sna); - nmax = atom->nmax; - memory->create(sna,nmax,size_peratom_cols,"sna/grid:sna"); - array_atom = sna; - } +// neighbor->build_one(list); - // invoke full neighbor list (will copy or build if necessary) +// const int inum = list->inum; +// const int* const ilist = list->ilist; +// const int* const numneigh = list->numneigh; +// int** const firstneigh = list->firstneigh; - neighbor->build_one(list); - - const int inum = list->inum; - const int* const ilist = list->ilist; - const int* const numneigh = list->numneigh; - int** const firstneigh = list->firstneigh; int * const type = atom->type; - // compute sna for each atom in group - // use full neighbor list to count atoms less than cutoff + // compute sna for each gridpoint double** const x = atom->x; const int* const mask = atom->mask; + const int ntotal = atom->nlocal + atom->nghost; - for (int ii = 0; ii < inum; ii++) { - const int i = ilist[ii]; - if (mask[i] & groupbit) { + printf("ngridfull = %d\n",ngridfull); + for (int igrid = 0; igrid < ngridfull; igrid++) { + printf("igrid = %d\n",igrid); + double rtmp[3]; + igridfull2x(igrid, rtmp); + const double xtmp = rtmp[0]; + const double ytmp = rtmp[1]; + const double ztmp = rtmp[2]; - const double xtmp = x[i][0]; - const double ytmp = x[i][1]; - const double ztmp = x[i][2]; - const int itype = type[i]; - const double radi = radelem[itype]; - const int* const jlist = firstneigh[i]; - const int jnum = numneigh[i]; + // rij[][3] = displacements between atom I and those neighbors + // inside = indices of neighbors of I within cutoff + // typej = types of neighbors of I within cutoff - // insure rij, inside, and typej are of size jnum + int ninside = 0; + for (int j = 0; j < ntotal; j++) { - snaptr->grow_rij(jnum); + // check that j is in comute group - // rij[][3] = displacements between atom I and those neighbors - // inside = indices of neighbors of I within cutoff - // typej = types of neighbors of I within cutoff + if (!(mask[j] & groupbit)) continue; - int ninside = 0; - for (int jj = 0; jj < jnum; jj++) { - int j = jlist[jj]; - j &= NEIGHMASK; + // insure rij, inside, and typej are of size jnum - 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; - int jtype = type[j]; - if (rsq < cutsq[itype][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] = (radi+radelem[jtype])*rcutfac; - ninside++; - } + snaptr->grow_rij(ninside+1); + + 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; + int jtype = type[j]; + if (rsq < cutsq[jtype][jtype] && rsq>1e-20) { + printf("ninside = %d\n",ninside); + 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; + ninside++; } + } - snaptr->compute_ui(ninside); - snaptr->compute_zi(); - snaptr->compute_bi(); - for (int icoeff = 0; icoeff < ncoeff; icoeff++) - sna[i][icoeff] = snaptr->blist[icoeff]; - if (quadraticflag) { - int ncount = ncoeff; - for (int icoeff = 0; icoeff < ncoeff; icoeff++) { - double bi = snaptr->blist[icoeff]; + snaptr->compute_ui(ninside); + snaptr->compute_zi(); + snaptr->compute_bi(); + for (int icoeff = 0; icoeff < ncoeff; icoeff++) + sna[igrid][icoeff] = snaptr->blist[icoeff]; + printf("igrid = %d B0 = %g\n",igrid,sna[igrid][0]); + if (quadraticflag) { + int ncount = ncoeff; + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + double bi = snaptr->blist[icoeff]; - // diagonal element of quadratic matrix + // diagonal element of quadratic matrix - sna[i][ncount++] = 0.5*bi*bi; + sna[igrid][ncount++] = 0.5*bi*bi; - // upper-triangular elements of quadratic matrix + // upper-triangular elements of quadratic matrix - for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) - sna[i][ncount++] = bi*snaptr->blist[jcoeff]; - } + for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) + sna[igrid][ncount++] = bi*snaptr->blist[jcoeff]; } - } else { - for (int icoeff = 0; icoeff < size_peratom_cols; icoeff++) - sna[i][icoeff] = 0.0; } } + gather_global_array(); +} + +/* ---------------------------------------------------------------------- + allocate array in base class and then set up pointers +------------------------------------------------------------------------- */ + +void ComputeSNAGrid::allocate() +{ + ComputeGrid::allocate(); + sna = gridfull; } /* ---------------------------------------------------------------------- @@ -269,7 +281,7 @@ void ComputeSNAGrid::compute_pergrid() double ComputeSNAGrid::memory_usage() { - double bytes = nmax*size_peratom_cols * sizeof(double); // sna + double bytes = size_array_rows*size_array_cols * sizeof(double); // grid bytes += snaptr->memory_usage(); // SNA object return bytes; diff --git a/src/compute_sna_grid.h b/src/compute_sna_grid.h index f85d9679cc..2bd18a915e 100644 --- a/src/compute_sna_grid.h +++ b/src/compute_sna_grid.h @@ -20,18 +20,19 @@ ComputeStyle(sna/grid,ComputeSNAGrid) #ifndef LMP_COMPUTE_SNA_GRID_H #define LMP_COMPUTE_SNA_GRID_H -#include "compute.h" +#include "compute_grid.h" namespace LAMMPS_NS { -class ComputeSNAGrid : public Compute { +class ComputeSNAGrid : public ComputeGrid { public: ComputeSNAGrid(class LAMMPS *, int, char **); ~ComputeSNAGrid(); void init(); void init_list(int, class NeighList *); - void compute_grid(); + void compute_array(); double memory_usage(); + void allocate(); private: int nmax; @@ -43,7 +44,6 @@ class ComputeSNAGrid : public Compute { double *radelem; double *wjelem; class SNA* snaptr; - double cutmax; int quadraticflag; }; From 8374280383487f86c1f69f437501e1a77c253e60 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Wed, 23 Oct 2019 18:46:28 -0600 Subject: [PATCH 03/19] Got a first pass working for ortho and tri grids --- examples/snap/in.grid | 71 +++++++++++++++++ examples/snap/in.grid.tri | 116 ++++++++++++++++++++++++++++ src/{ => SNAP}/compute_sna_grid.cpp | 0 src/{ => SNAP}/compute_sna_grid.h | 0 src/compute_grid.cpp | 56 ++++++++++---- src/compute_grid.h | 2 + 6 files changed, 231 insertions(+), 14 deletions(-) create mode 100644 examples/snap/in.grid create mode 100644 examples/snap/in.grid.tri rename src/{ => SNAP}/compute_sna_grid.cpp (100%) rename src/{ => SNAP}/compute_sna_grid.h (100%) diff --git a/examples/snap/in.grid b/examples/snap/in.grid new file mode 100644 index 0000000000..bec57dc9ff --- /dev/null +++ b/examples/snap/in.grid @@ -0,0 +1,71 @@ +# Demonstrate bispectrum computes + +# Initialize simulation + +variable nsteps index 0 +variable nrep index 1 +variable a index 3.316 +variable ngrid index 2 + +units metal + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrep} +variable ny equal ${nrep} +variable nz equal ${nrep} + +boundary p p p + +lattice custom $a & + a1 1 0 0 & + a2 0 1 0 & + a3 0 0 1 & + basis 0 0 0 & + basis 0.5 0.5 0.5 & +# origin 0.25 0.25 0.25 + +region box block 0 ${nx} 0 ${ny} 0 ${nz} +create_box 1 box +create_atoms 1 box + +mass 1 180.88 + +# choose potential + +include Ta06A.snap + +# define grid compute and atom compute + +group snapgroup type 1 +variable twojmax equal 2 +variable rcutfac equal 4.67637 +variable rfac0 equal 0.99363 +variable rmin0 equal 0 +variable wj equal 1 +variable radelem equal 0.5 +variable bzero equal 0 +variable quad equal 0 +variable switch equal 0 + +compute b all sna/atom & + ${rcutfac} ${rfac0} ${twojmax} ${radelem} & + ${wj} rmin0 ${rmin0} bzeroflag ${bzero} & + quadraticflag ${quad} switchflag ${switch} + +compute mygrid all sna/grid grid ${ngrid} ${ngrid} ${ngrid} & + ${rcutfac} ${rfac0} ${twojmax} ${radelem} & + ${wj} rmin0 ${rmin0} bzeroflag ${bzero} & + quadraticflag ${quad} switchflag ${switch} + +# define output + +thermo_style custom step temp ke pe vol c_mygrid[1][1] c_mygrid[2][1] c_mygrid[3][1] c_mygrid[4][1] c_mygrid[5][1] c_mygrid[6][1] c_mygrid[7][1] c_mygrid[8][1] +thermo_modify norm yes + +dump mydump_b all custom 1000 dump_b id c_b[*] + +# run + +run 0 + diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri new file mode 100644 index 0000000000..c02138f66a --- /dev/null +++ b/examples/snap/in.grid.tri @@ -0,0 +1,116 @@ +# Demonstrate bispectrum computes + +# Initialize simulation + +variable nsteps index 0 +variable nrep index 1 +variable a index 3.316 +variable ngrid index 4 + +variable ngridx equal 3*${ngrid} +variable ngridy equal 2*${ngrid} +variable ngridz equal 1*${ngrid} + +units metal + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrep} +variable ny equal ${nrep} +variable nz equal ${nrep} + +boundary p p p + +lattice custom $a & + a1 1 0 0 & + a2 1 1 0 & + a3 1 1 1 & + basis 0 0 0 & + basis 0.0 0.0 0.5 & +# origin 0.25 0.25 0.25 + +box tilt large +region box prism 0 ${nx} 0 ${ny} 0 ${nz} ${ny} ${nz} ${nz} +create_box 1 box +create_atoms 1 box + +mass 1 180.88 + +# choose potential + +include Ta06A.snap + +# define grid compute and atom compute + +group snapgroup type 1 +variable twojmax equal 2 +variable rcutfac equal 4.67637 +variable rfac0 equal 0.99363 +variable rmin0 equal 0 +variable wj equal 1 +variable radelem equal 0.5 +variable bzero equal 0 +variable quad equal 0 +variable switch equal 0 + +compute b all sna/atom & + ${rcutfac} ${rfac0} ${twojmax} ${radelem} & + ${wj} rmin0 ${rmin0} bzeroflag ${bzero} & + quadraticflag ${quad} switchflag ${switch} + +compute mygrid all sna/grid grid ${ngridx} ${ngridy} ${ngridz} & + ${rcutfac} ${rfac0} ${twojmax} ${radelem} & + ${wj} rmin0 ${rmin0} bzeroflag ${bzero} & + quadraticflag ${quad} switchflag ${switch} + +# define output + +thermo_style custom step temp ke pe vol & + c_mygrid[1][1] c_mygrid[2][1] c_mygrid[3][1] c_mygrid[4][1] c_mygrid[5][1] c_mygrid[6][1] c_mygrid[7][1] c_mygrid[8][1] c_mygrid[9][1] & + c_mygrid[10][1] c_mygrid[11][1] c_mygrid[12][1] c_mygrid[13][1] c_mygrid[14][1] c_mygrid[15][1] c_mygrid[16][1] c_mygrid[17][1] c_mygrid[18][1] c_mygrid[19][1] & + c_mygrid[20][1] c_mygrid[21][1] c_mygrid[22][1] c_mygrid[23][1] c_mygrid[24][1] c_mygrid[25][1] c_mygrid[26][1] c_mygrid[27][1] c_mygrid[28][1] c_mygrid[29][1] & + c_mygrid[30][1] c_mygrid[31][1] c_mygrid[32][1] c_mygrid[33][1] c_mygrid[34][1] c_mygrid[35][1] c_mygrid[36][1] c_mygrid[37][1] c_mygrid[38][1] c_mygrid[39][1] & + c_mygrid[40][1] c_mygrid[41][1] c_mygrid[42][1] c_mygrid[43][1] c_mygrid[44][1] c_mygrid[45][1] c_mygrid[46][1] c_mygrid[47][1] c_mygrid[48][1] c_mygrid[49][1] & + c_mygrid[50][1] c_mygrid[51][1] c_mygrid[52][1] c_mygrid[53][1] c_mygrid[54][1] c_mygrid[55][1] c_mygrid[56][1] c_mygrid[57][1] c_mygrid[58][1] c_mygrid[59][1] & + c_mygrid[60][1] c_mygrid[61][1] c_mygrid[62][1] c_mygrid[63][1] c_mygrid[64][1] c_mygrid[65][1] c_mygrid[66][1] c_mygrid[67][1] c_mygrid[68][1] c_mygrid[69][1] & + c_mygrid[70][1] c_mygrid[71][1] c_mygrid[72][1] c_mygrid[73][1] c_mygrid[74][1] c_mygrid[75][1] c_mygrid[76][1] c_mygrid[77][1] c_mygrid[78][1] c_mygrid[79][1] & + c_mygrid[80][1] c_mygrid[81][1] c_mygrid[82][1] c_mygrid[83][1] c_mygrid[84][1] c_mygrid[85][1] c_mygrid[86][1] c_mygrid[87][1] c_mygrid[88][1] c_mygrid[89][1] & + c_mygrid[90][1] c_mygrid[91][1] c_mygrid[92][1] c_mygrid[93][1] c_mygrid[94][1] c_mygrid[95][1] c_mygrid[96][1] c_mygrid[97][1] c_mygrid[98][1] c_mygrid[99][1] & + c_mygrid[101][1] c_mygrid[102][1] c_mygrid[103][1] c_mygrid[104][1] c_mygrid[105][1] c_mygrid[106][1] c_mygrid[107][1] c_mygrid[108][1] c_mygrid[109][1] & + c_mygrid[110][1] c_mygrid[111][1] c_mygrid[112][1] c_mygrid[113][1] c_mygrid[114][1] c_mygrid[115][1] c_mygrid[116][1] c_mygrid[117][1] c_mygrid[118][1] c_mygrid[119][1] & + c_mygrid[120][1] c_mygrid[121][1] c_mygrid[122][1] c_mygrid[123][1] c_mygrid[124][1] c_mygrid[125][1] c_mygrid[126][1] c_mygrid[127][1] c_mygrid[128][1] c_mygrid[129][1] & + c_mygrid[130][1] c_mygrid[131][1] c_mygrid[132][1] c_mygrid[133][1] c_mygrid[134][1] c_mygrid[135][1] c_mygrid[136][1] c_mygrid[137][1] c_mygrid[138][1] c_mygrid[139][1] & + c_mygrid[140][1] c_mygrid[141][1] c_mygrid[142][1] c_mygrid[143][1] c_mygrid[144][1] c_mygrid[145][1] c_mygrid[146][1] c_mygrid[147][1] c_mygrid[148][1] c_mygrid[149][1] & + c_mygrid[150][1] c_mygrid[151][1] c_mygrid[152][1] c_mygrid[153][1] c_mygrid[154][1] c_mygrid[155][1] c_mygrid[156][1] c_mygrid[157][1] c_mygrid[158][1] c_mygrid[159][1] & + c_mygrid[160][1] c_mygrid[161][1] c_mygrid[162][1] c_mygrid[163][1] c_mygrid[164][1] c_mygrid[165][1] c_mygrid[166][1] c_mygrid[167][1] c_mygrid[168][1] c_mygrid[169][1] & + c_mygrid[170][1] c_mygrid[171][1] c_mygrid[172][1] c_mygrid[173][1] c_mygrid[174][1] c_mygrid[175][1] c_mygrid[176][1] c_mygrid[177][1] c_mygrid[178][1] c_mygrid[179][1] & + c_mygrid[180][1] c_mygrid[181][1] c_mygrid[182][1] c_mygrid[183][1] c_mygrid[184][1] c_mygrid[185][1] c_mygrid[186][1] c_mygrid[187][1] c_mygrid[188][1] c_mygrid[189][1] & + c_mygrid[190][1] c_mygrid[191][1] c_mygrid[192][1] c_mygrid[193][1] c_mygrid[194][1] c_mygrid[195][1] c_mygrid[196][1] c_mygrid[197][1] c_mygrid[198][1] c_mygrid[199][1] & + c_mygrid[201][1] c_mygrid[202][1] c_mygrid[203][1] c_mygrid[204][1] c_mygrid[205][1] c_mygrid[206][1] c_mygrid[207][1] c_mygrid[208][1] c_mygrid[209][1] & + c_mygrid[210][1] c_mygrid[211][1] c_mygrid[212][1] c_mygrid[213][1] c_mygrid[214][1] c_mygrid[215][1] c_mygrid[216][1] c_mygrid[217][1] c_mygrid[218][1] c_mygrid[219][1] & + c_mygrid[220][1] c_mygrid[221][1] c_mygrid[222][1] c_mygrid[223][1] c_mygrid[224][1] c_mygrid[225][1] c_mygrid[226][1] c_mygrid[227][1] c_mygrid[228][1] c_mygrid[229][1] & + c_mygrid[230][1] c_mygrid[231][1] c_mygrid[232][1] c_mygrid[233][1] c_mygrid[234][1] c_mygrid[235][1] c_mygrid[236][1] c_mygrid[237][1] c_mygrid[238][1] c_mygrid[239][1] & + c_mygrid[240][1] c_mygrid[241][1] c_mygrid[242][1] c_mygrid[243][1] c_mygrid[244][1] c_mygrid[245][1] c_mygrid[246][1] c_mygrid[247][1] c_mygrid[248][1] c_mygrid[249][1] & + c_mygrid[250][1] c_mygrid[251][1] c_mygrid[252][1] c_mygrid[253][1] c_mygrid[254][1] c_mygrid[255][1] c_mygrid[256][1] c_mygrid[257][1] c_mygrid[258][1] c_mygrid[259][1] & + c_mygrid[260][1] c_mygrid[261][1] c_mygrid[262][1] c_mygrid[263][1] c_mygrid[264][1] c_mygrid[265][1] c_mygrid[266][1] c_mygrid[267][1] c_mygrid[268][1] c_mygrid[269][1] & + c_mygrid[270][1] c_mygrid[271][1] c_mygrid[272][1] c_mygrid[273][1] c_mygrid[274][1] c_mygrid[275][1] c_mygrid[276][1] c_mygrid[277][1] c_mygrid[278][1] c_mygrid[279][1] & + c_mygrid[280][1] c_mygrid[281][1] c_mygrid[282][1] c_mygrid[283][1] c_mygrid[284][1] c_mygrid[285][1] c_mygrid[286][1] c_mygrid[287][1] c_mygrid[288][1] c_mygrid[289][1] & + c_mygrid[290][1] c_mygrid[291][1] c_mygrid[292][1] c_mygrid[293][1] c_mygrid[294][1] c_mygrid[295][1] c_mygrid[296][1] c_mygrid[297][1] c_mygrid[298][1] c_mygrid[299][1] & + c_mygrid[301][1] c_mygrid[302][1] c_mygrid[303][1] c_mygrid[304][1] c_mygrid[305][1] c_mygrid[306][1] c_mygrid[307][1] c_mygrid[308][1] c_mygrid[309][1] & + c_mygrid[310][1] c_mygrid[311][1] c_mygrid[312][1] c_mygrid[313][1] c_mygrid[314][1] c_mygrid[315][1] c_mygrid[316][1] c_mygrid[317][1] c_mygrid[318][1] c_mygrid[319][1] & + c_mygrid[320][1] c_mygrid[321][1] c_mygrid[322][1] c_mygrid[323][1] c_mygrid[324][1] c_mygrid[325][1] c_mygrid[326][1] c_mygrid[327][1] c_mygrid[328][1] c_mygrid[329][1] & + c_mygrid[330][1] c_mygrid[331][1] c_mygrid[332][1] c_mygrid[333][1] c_mygrid[334][1] c_mygrid[335][1] c_mygrid[336][1] c_mygrid[337][1] c_mygrid[338][1] c_mygrid[339][1] & + c_mygrid[340][1] c_mygrid[341][1] c_mygrid[342][1] c_mygrid[343][1] c_mygrid[344][1] c_mygrid[345][1] c_mygrid[346][1] c_mygrid[347][1] c_mygrid[348][1] c_mygrid[349][1] & + c_mygrid[350][1] c_mygrid[351][1] c_mygrid[352][1] c_mygrid[353][1] c_mygrid[354][1] c_mygrid[355][1] c_mygrid[356][1] c_mygrid[357][1] c_mygrid[358][1] c_mygrid[359][1] & + c_mygrid[360][1] c_mygrid[361][1] c_mygrid[362][1] c_mygrid[363][1] c_mygrid[364][1] c_mygrid[365][1] c_mygrid[366][1] c_mygrid[367][1] c_mygrid[368][1] c_mygrid[369][1] & + c_mygrid[370][1] c_mygrid[371][1] c_mygrid[372][1] c_mygrid[373][1] c_mygrid[374][1] c_mygrid[375][1] c_mygrid[376][1] c_mygrid[377][1] c_mygrid[378][1] c_mygrid[379][1] & + c_mygrid[380][1] c_mygrid[381][1] c_mygrid[382][1] c_mygrid[383][1] c_mygrid[384][1] + +thermo_modify norm yes + +dump mydump_b all custom 1000 dump_b id c_b[*] + +# run + +run 0 + diff --git a/src/compute_sna_grid.cpp b/src/SNAP/compute_sna_grid.cpp similarity index 100% rename from src/compute_sna_grid.cpp rename to src/SNAP/compute_sna_grid.cpp diff --git a/src/compute_sna_grid.h b/src/SNAP/compute_sna_grid.h similarity index 100% rename from src/compute_sna_grid.h rename to src/SNAP/compute_sna_grid.h diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp index d8c7eb3697..ec969fd01f 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -34,7 +34,7 @@ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : array_flag = 1; size_array_cols = 0; size_array_rows = 0; - extarray = 1; + extarray = 0; int iarg0 = 3; int iarg = iarg0; @@ -119,10 +119,6 @@ void ComputeGrid::setup() + pow(delxy*my,2))*delxinvtmp + 1; } - printf("mx = %d\n",mx); - printf("my = %d\n",my); - printf("mz = %d\n",mz); - // size global grid to accomodate periodic interactions nxfull = nx + 2*mx; @@ -130,10 +126,6 @@ void ComputeGrid::setup() nzfull = nz + 2*mz; nxyfull = nxfull * nyfull; - printf("nxfull = %d\n",nxfull); - printf("nyfull = %d\n",nyfull); - printf("nzfull = %d\n",nzfull); - x0full = boxlo[0] - mx*delx; y0full = boxlo[1] - my*dely; z0full = boxlo[2] - mz*delz; @@ -148,9 +140,9 @@ void ComputeGrid::setup() void ComputeGrid::igridfull2x(int igrid, double *x) { int iz = igrid / nxyfull; - igrid -= iz*nxyfull; + igrid -= iz * nxyfull; int iy = igrid / nxfull; - igrid -= igrid*nxfull; + igrid -= iy * nxfull; int ix = igrid; x[0] = x0full+ix*delx; @@ -161,14 +153,28 @@ void ComputeGrid::igridfull2x(int igrid, double *x) } +/* ---------------------------------------------------------------------- + copy local grid to global array +------------------------------------------------------------------------- */ + +void ComputeGrid::copy_local_grid() +{ + int igridfull; + for (int iarray = 0; iarray < size_array_rows; iarray++) { + igridfull = iarray2igridfull(iarray); + for (int icol = 0; icol < size_array_cols; icol++) + array[iarray][icol] = gridfull[igridfull][icol]; + } +} + /* ---------------------------------------------------------------------- gather global array from full grid ------------------------------------------------------------------------- */ void ComputeGrid::gather_global_array() { - int iarray; - memset(&array[0][0],0,size_array_rows*size_array_cols*sizeof(double)); + int iarray; + memset(&array[0][0],0,size_array_rows*size_array_cols*sizeof(double)); for (int igrid = 0; igrid < ngridfull; igrid++) { @@ -190,7 +196,7 @@ int ComputeGrid::igridfull2iarray(int igrid) int iz = igrid / nxyfull; igrid -= iz*nxyfull; int iy = igrid / nxfull; - igrid -= igrid*nxfull; + igrid -= iy*nxfull; int ix = igrid; ix -= mx; @@ -210,6 +216,28 @@ int ComputeGrid::igridfull2iarray(int igrid) return iarray; } +/* ---------------------------------------------------------------------- + convert compute array index to full grid index + inefficient, should exploit shared ix structure +------------------------------------------------------------------------- */ + +int ComputeGrid::iarray2igridfull(int iarray) +{ + int iz = iarray / (nx*ny); + iarray -= iz*(nx*ny); + int iy = iarray / nx; + iarray -= iy*nx; + int ix = iarray; + + ix += mx; + iy += my; + iz += mz; + + int igrid = (iz * nyfull + iy) * nxfull + ix; + + return igrid; +} + /* ---------------------------------------------------------------------- free and reallocate arrays ------------------------------------------------------------------------- */ diff --git a/src/compute_grid.h b/src/compute_grid.h index 61775186ff..bc757c37bc 100644 --- a/src/compute_grid.h +++ b/src/compute_grid.h @@ -46,7 +46,9 @@ class ComputeGrid : public Compute { virtual void allocate(); void igridfull2x(int, double*); // convert full grid point to coord void gather_global_array(); // gather global array from full grid + void copy_local_grid(); // copy local grid to global array int igridfull2iarray(int); // convert full grid index to compute array index + int iarray2igridfull(int); // convert compute array index to full grid index private: }; From 0fc325c98b2b4ffd8fb24cd53ce732221add6ac0 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Wed, 23 Oct 2019 18:56:21 -0600 Subject: [PATCH 04/19] Got a first pass working for ortho and tri grids --- src/SNAP/compute_sna_grid.cpp | 24 +++++++++--------------- 1 file changed, 9 insertions(+), 15 deletions(-) diff --git a/src/SNAP/compute_sna_grid.cpp b/src/SNAP/compute_sna_grid.cpp index 1f998ced9b..26bd272d9f 100644 --- a/src/SNAP/compute_sna_grid.cpp +++ b/src/SNAP/compute_sna_grid.cpp @@ -66,7 +66,6 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : rcutfac = atof(arg[3]); rfac0 = atof(arg[4]); twojmax = atoi(arg[5]); - printf("rcutfac = %g rfac0 = %g twojmax = %d\n",rcutfac, rfac0, twojmax); for(int i = 0; i < ntypes; i++) radelem[i+1] = atof(arg[6+i]); @@ -82,7 +81,6 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : cut = 2.0*radelem[i]*rcutfac; if (cut > cutmax) cutmax = cut; cutsq[i][i] = cut*cut; - printf("i = %d cutsq[i][i] = %g\n",i,cutsq[i][i]); for(int j = i+1; j <= ntypes; j++) { cut = (radelem[i]+radelem[j])*rcutfac; cutsq[i][j] = cutsq[j][i] = cut*cut; @@ -94,7 +92,6 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : int iarg = nargmin; while (iarg < narg) { - printf("iarg = %d arg = %s\n",iarg, arg[iarg]); if (strcmp(arg[iarg],"rmin0") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute sna/grid command"); @@ -119,9 +116,6 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : } - printf("rmin0 = %g, bzeroflag = %d, quadraticflag = %d\n", - rmin0, bzeroflag, quadraticflag); - snaptr = new SNA(lmp,rfac0,twojmax, rmin0,switchflag,bzeroflag); @@ -199,9 +193,12 @@ void ComputeSNAGrid::compute_array() const int* const mask = atom->mask; const int ntotal = atom->nlocal + atom->nghost; + // insure rij, inside, and typej are of size jnum + + snaptr->grow_rij(ntotal); + printf("ngridfull = %d\n",ngridfull); for (int igrid = 0; igrid < ngridfull; igrid++) { - printf("igrid = %d\n",igrid); double rtmp[3]; igridfull2x(igrid, rtmp); const double xtmp = rtmp[0]; @@ -215,21 +212,17 @@ void ComputeSNAGrid::compute_array() int ninside = 0; for (int j = 0; j < ntotal; j++) { - // check that j is in comute group + // check that j is in compute group if (!(mask[j] & groupbit)) continue; - // insure rij, inside, and typej are of size jnum - - snaptr->grow_rij(ninside+1); - 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; int jtype = type[j]; if (rsq < cutsq[jtype][jtype] && rsq>1e-20) { - printf("ninside = %d\n",ninside); + // printf("ninside = %d\n",ninside); snaptr->rij[ninside][0] = delx; snaptr->rij[ninside][1] = dely; snaptr->rij[ninside][2] = delz; @@ -245,7 +238,7 @@ void ComputeSNAGrid::compute_array() snaptr->compute_bi(); for (int icoeff = 0; icoeff < ncoeff; icoeff++) sna[igrid][icoeff] = snaptr->blist[icoeff]; - printf("igrid = %d B0 = %g\n",igrid,sna[igrid][0]); + // printf("igrid = %d %g %g %g %d B0 = %g\n",igrid,xtmp,ytmp,ztmp,ninside,sna[igrid][0]); if (quadraticflag) { int ncount = ncoeff; for (int icoeff = 0; icoeff < ncoeff; icoeff++) { @@ -262,7 +255,8 @@ void ComputeSNAGrid::compute_array() } } } - gather_global_array(); + // gather_global_array(); + copy_local_grid(); } /* ---------------------------------------------------------------------- From 5956908cfd2bbd097be8126cc4f062d5cc9cc519 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Wed, 23 Oct 2019 21:51:36 -0600 Subject: [PATCH 05/19] Added cols for coords --- examples/snap/in.grid | 6 ++- examples/snap/in.grid.tri | 78 +++++++++++++++++------------------ src/SNAP/compute_sna_grid.cpp | 6 +-- src/compute_grid.cpp | 41 +++++++++++++++++- src/compute_grid.h | 5 ++- 5 files changed, 90 insertions(+), 46 deletions(-) diff --git a/examples/snap/in.grid b/examples/snap/in.grid index bec57dc9ff..5293558fab 100644 --- a/examples/snap/in.grid +++ b/examples/snap/in.grid @@ -60,7 +60,11 @@ compute mygrid all sna/grid grid ${ngrid} ${ngrid} ${ngrid} & # define output -thermo_style custom step temp ke pe vol c_mygrid[1][1] c_mygrid[2][1] c_mygrid[3][1] c_mygrid[4][1] c_mygrid[5][1] c_mygrid[6][1] c_mygrid[7][1] c_mygrid[8][1] +thermo_style custom step temp ke pe vol & + c_mygrid[1][1] c_mygrid[2][1] c_mygrid[3][1] c_mygrid[4][1] c_mygrid[5][1] c_mygrid[6][1] c_mygrid[7][1] c_mygrid[8][1] & + c_mygrid[1][2] c_mygrid[2][2] c_mygrid[3][2] c_mygrid[4][2] c_mygrid[5][2] c_mygrid[6][2] c_mygrid[7][2] c_mygrid[8][2] & + c_mygrid[1][3] c_mygrid[2][3] c_mygrid[3][3] c_mygrid[4][3] c_mygrid[5][3] c_mygrid[6][3] c_mygrid[7][3] c_mygrid[8][3] & + c_mygrid[1][4] c_mygrid[2][4] c_mygrid[3][4] c_mygrid[4][4] c_mygrid[5][4] c_mygrid[6][4] c_mygrid[7][4] c_mygrid[8][4] thermo_modify norm yes dump mydump_b all custom 1000 dump_b id c_b[*] diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri index c02138f66a..3c471d9b66 100644 --- a/examples/snap/in.grid.tri +++ b/examples/snap/in.grid.tri @@ -66,45 +66,45 @@ compute mygrid all sna/grid grid ${ngridx} ${ngridy} ${ngridz} & # define output thermo_style custom step temp ke pe vol & - c_mygrid[1][1] c_mygrid[2][1] c_mygrid[3][1] c_mygrid[4][1] c_mygrid[5][1] c_mygrid[6][1] c_mygrid[7][1] c_mygrid[8][1] c_mygrid[9][1] & - c_mygrid[10][1] c_mygrid[11][1] c_mygrid[12][1] c_mygrid[13][1] c_mygrid[14][1] c_mygrid[15][1] c_mygrid[16][1] c_mygrid[17][1] c_mygrid[18][1] c_mygrid[19][1] & - c_mygrid[20][1] c_mygrid[21][1] c_mygrid[22][1] c_mygrid[23][1] c_mygrid[24][1] c_mygrid[25][1] c_mygrid[26][1] c_mygrid[27][1] c_mygrid[28][1] c_mygrid[29][1] & - c_mygrid[30][1] c_mygrid[31][1] c_mygrid[32][1] c_mygrid[33][1] c_mygrid[34][1] c_mygrid[35][1] c_mygrid[36][1] c_mygrid[37][1] c_mygrid[38][1] c_mygrid[39][1] & - c_mygrid[40][1] c_mygrid[41][1] c_mygrid[42][1] c_mygrid[43][1] c_mygrid[44][1] c_mygrid[45][1] c_mygrid[46][1] c_mygrid[47][1] c_mygrid[48][1] c_mygrid[49][1] & - c_mygrid[50][1] c_mygrid[51][1] c_mygrid[52][1] c_mygrid[53][1] c_mygrid[54][1] c_mygrid[55][1] c_mygrid[56][1] c_mygrid[57][1] c_mygrid[58][1] c_mygrid[59][1] & - c_mygrid[60][1] c_mygrid[61][1] c_mygrid[62][1] c_mygrid[63][1] c_mygrid[64][1] c_mygrid[65][1] c_mygrid[66][1] c_mygrid[67][1] c_mygrid[68][1] c_mygrid[69][1] & - c_mygrid[70][1] c_mygrid[71][1] c_mygrid[72][1] c_mygrid[73][1] c_mygrid[74][1] c_mygrid[75][1] c_mygrid[76][1] c_mygrid[77][1] c_mygrid[78][1] c_mygrid[79][1] & - c_mygrid[80][1] c_mygrid[81][1] c_mygrid[82][1] c_mygrid[83][1] c_mygrid[84][1] c_mygrid[85][1] c_mygrid[86][1] c_mygrid[87][1] c_mygrid[88][1] c_mygrid[89][1] & - c_mygrid[90][1] c_mygrid[91][1] c_mygrid[92][1] c_mygrid[93][1] c_mygrid[94][1] c_mygrid[95][1] c_mygrid[96][1] c_mygrid[97][1] c_mygrid[98][1] c_mygrid[99][1] & - c_mygrid[101][1] c_mygrid[102][1] c_mygrid[103][1] c_mygrid[104][1] c_mygrid[105][1] c_mygrid[106][1] c_mygrid[107][1] c_mygrid[108][1] c_mygrid[109][1] & - c_mygrid[110][1] c_mygrid[111][1] c_mygrid[112][1] c_mygrid[113][1] c_mygrid[114][1] c_mygrid[115][1] c_mygrid[116][1] c_mygrid[117][1] c_mygrid[118][1] c_mygrid[119][1] & - c_mygrid[120][1] c_mygrid[121][1] c_mygrid[122][1] c_mygrid[123][1] c_mygrid[124][1] c_mygrid[125][1] c_mygrid[126][1] c_mygrid[127][1] c_mygrid[128][1] c_mygrid[129][1] & - c_mygrid[130][1] c_mygrid[131][1] c_mygrid[132][1] c_mygrid[133][1] c_mygrid[134][1] c_mygrid[135][1] c_mygrid[136][1] c_mygrid[137][1] c_mygrid[138][1] c_mygrid[139][1] & - c_mygrid[140][1] c_mygrid[141][1] c_mygrid[142][1] c_mygrid[143][1] c_mygrid[144][1] c_mygrid[145][1] c_mygrid[146][1] c_mygrid[147][1] c_mygrid[148][1] c_mygrid[149][1] & - c_mygrid[150][1] c_mygrid[151][1] c_mygrid[152][1] c_mygrid[153][1] c_mygrid[154][1] c_mygrid[155][1] c_mygrid[156][1] c_mygrid[157][1] c_mygrid[158][1] c_mygrid[159][1] & - c_mygrid[160][1] c_mygrid[161][1] c_mygrid[162][1] c_mygrid[163][1] c_mygrid[164][1] c_mygrid[165][1] c_mygrid[166][1] c_mygrid[167][1] c_mygrid[168][1] c_mygrid[169][1] & - c_mygrid[170][1] c_mygrid[171][1] c_mygrid[172][1] c_mygrid[173][1] c_mygrid[174][1] c_mygrid[175][1] c_mygrid[176][1] c_mygrid[177][1] c_mygrid[178][1] c_mygrid[179][1] & - c_mygrid[180][1] c_mygrid[181][1] c_mygrid[182][1] c_mygrid[183][1] c_mygrid[184][1] c_mygrid[185][1] c_mygrid[186][1] c_mygrid[187][1] c_mygrid[188][1] c_mygrid[189][1] & - c_mygrid[190][1] c_mygrid[191][1] c_mygrid[192][1] c_mygrid[193][1] c_mygrid[194][1] c_mygrid[195][1] c_mygrid[196][1] c_mygrid[197][1] c_mygrid[198][1] c_mygrid[199][1] & - c_mygrid[201][1] c_mygrid[202][1] c_mygrid[203][1] c_mygrid[204][1] c_mygrid[205][1] c_mygrid[206][1] c_mygrid[207][1] c_mygrid[208][1] c_mygrid[209][1] & - c_mygrid[210][1] c_mygrid[211][1] c_mygrid[212][1] c_mygrid[213][1] c_mygrid[214][1] c_mygrid[215][1] c_mygrid[216][1] c_mygrid[217][1] c_mygrid[218][1] c_mygrid[219][1] & - c_mygrid[220][1] c_mygrid[221][1] c_mygrid[222][1] c_mygrid[223][1] c_mygrid[224][1] c_mygrid[225][1] c_mygrid[226][1] c_mygrid[227][1] c_mygrid[228][1] c_mygrid[229][1] & - c_mygrid[230][1] c_mygrid[231][1] c_mygrid[232][1] c_mygrid[233][1] c_mygrid[234][1] c_mygrid[235][1] c_mygrid[236][1] c_mygrid[237][1] c_mygrid[238][1] c_mygrid[239][1] & - c_mygrid[240][1] c_mygrid[241][1] c_mygrid[242][1] c_mygrid[243][1] c_mygrid[244][1] c_mygrid[245][1] c_mygrid[246][1] c_mygrid[247][1] c_mygrid[248][1] c_mygrid[249][1] & - c_mygrid[250][1] c_mygrid[251][1] c_mygrid[252][1] c_mygrid[253][1] c_mygrid[254][1] c_mygrid[255][1] c_mygrid[256][1] c_mygrid[257][1] c_mygrid[258][1] c_mygrid[259][1] & - c_mygrid[260][1] c_mygrid[261][1] c_mygrid[262][1] c_mygrid[263][1] c_mygrid[264][1] c_mygrid[265][1] c_mygrid[266][1] c_mygrid[267][1] c_mygrid[268][1] c_mygrid[269][1] & - c_mygrid[270][1] c_mygrid[271][1] c_mygrid[272][1] c_mygrid[273][1] c_mygrid[274][1] c_mygrid[275][1] c_mygrid[276][1] c_mygrid[277][1] c_mygrid[278][1] c_mygrid[279][1] & - c_mygrid[280][1] c_mygrid[281][1] c_mygrid[282][1] c_mygrid[283][1] c_mygrid[284][1] c_mygrid[285][1] c_mygrid[286][1] c_mygrid[287][1] c_mygrid[288][1] c_mygrid[289][1] & - c_mygrid[290][1] c_mygrid[291][1] c_mygrid[292][1] c_mygrid[293][1] c_mygrid[294][1] c_mygrid[295][1] c_mygrid[296][1] c_mygrid[297][1] c_mygrid[298][1] c_mygrid[299][1] & - c_mygrid[301][1] c_mygrid[302][1] c_mygrid[303][1] c_mygrid[304][1] c_mygrid[305][1] c_mygrid[306][1] c_mygrid[307][1] c_mygrid[308][1] c_mygrid[309][1] & - c_mygrid[310][1] c_mygrid[311][1] c_mygrid[312][1] c_mygrid[313][1] c_mygrid[314][1] c_mygrid[315][1] c_mygrid[316][1] c_mygrid[317][1] c_mygrid[318][1] c_mygrid[319][1] & - c_mygrid[320][1] c_mygrid[321][1] c_mygrid[322][1] c_mygrid[323][1] c_mygrid[324][1] c_mygrid[325][1] c_mygrid[326][1] c_mygrid[327][1] c_mygrid[328][1] c_mygrid[329][1] & - c_mygrid[330][1] c_mygrid[331][1] c_mygrid[332][1] c_mygrid[333][1] c_mygrid[334][1] c_mygrid[335][1] c_mygrid[336][1] c_mygrid[337][1] c_mygrid[338][1] c_mygrid[339][1] & - c_mygrid[340][1] c_mygrid[341][1] c_mygrid[342][1] c_mygrid[343][1] c_mygrid[344][1] c_mygrid[345][1] c_mygrid[346][1] c_mygrid[347][1] c_mygrid[348][1] c_mygrid[349][1] & - c_mygrid[350][1] c_mygrid[351][1] c_mygrid[352][1] c_mygrid[353][1] c_mygrid[354][1] c_mygrid[355][1] c_mygrid[356][1] c_mygrid[357][1] c_mygrid[358][1] c_mygrid[359][1] & - c_mygrid[360][1] c_mygrid[361][1] c_mygrid[362][1] c_mygrid[363][1] c_mygrid[364][1] c_mygrid[365][1] c_mygrid[366][1] c_mygrid[367][1] c_mygrid[368][1] c_mygrid[369][1] & - c_mygrid[370][1] c_mygrid[371][1] c_mygrid[372][1] c_mygrid[373][1] c_mygrid[374][1] c_mygrid[375][1] c_mygrid[376][1] c_mygrid[377][1] c_mygrid[378][1] c_mygrid[379][1] & - c_mygrid[380][1] c_mygrid[381][1] c_mygrid[382][1] c_mygrid[383][1] c_mygrid[384][1] + c_mygrid[1][4] c_mygrid[2][4] c_mygrid[3][4] c_mygrid[4][4] c_mygrid[5][4] c_mygrid[6][4] c_mygrid[7][4] c_mygrid[8][4] c_mygrid[9][4] & + c_mygrid[10][4] c_mygrid[11][4] c_mygrid[12][4] c_mygrid[13][4] c_mygrid[14][4] c_mygrid[15][4] c_mygrid[16][4] c_mygrid[17][4] c_mygrid[18][4] c_mygrid[19][4] & + c_mygrid[20][4] c_mygrid[21][4] c_mygrid[22][4] c_mygrid[23][4] c_mygrid[24][4] c_mygrid[25][4] c_mygrid[26][4] c_mygrid[27][4] c_mygrid[28][4] c_mygrid[29][4] & + c_mygrid[30][4] c_mygrid[31][4] c_mygrid[32][4] c_mygrid[33][4] c_mygrid[34][4] c_mygrid[35][4] c_mygrid[36][4] c_mygrid[37][4] c_mygrid[38][4] c_mygrid[39][4] & + c_mygrid[40][4] c_mygrid[41][4] c_mygrid[42][4] c_mygrid[43][4] c_mygrid[44][4] c_mygrid[45][4] c_mygrid[46][4] c_mygrid[47][4] c_mygrid[48][4] c_mygrid[49][4] & + c_mygrid[50][4] c_mygrid[51][4] c_mygrid[52][4] c_mygrid[53][4] c_mygrid[54][4] c_mygrid[55][4] c_mygrid[56][4] c_mygrid[57][4] c_mygrid[58][4] c_mygrid[59][4] & + c_mygrid[60][4] c_mygrid[61][4] c_mygrid[62][4] c_mygrid[63][4] c_mygrid[64][4] c_mygrid[65][4] c_mygrid[66][4] c_mygrid[67][4] c_mygrid[68][4] c_mygrid[69][4] & + c_mygrid[70][4] c_mygrid[71][4] c_mygrid[72][4] c_mygrid[73][4] c_mygrid[74][4] c_mygrid[75][4] c_mygrid[76][4] c_mygrid[77][4] c_mygrid[78][4] c_mygrid[79][4] & + c_mygrid[80][4] c_mygrid[81][4] c_mygrid[82][4] c_mygrid[83][4] c_mygrid[84][4] c_mygrid[85][4] c_mygrid[86][4] c_mygrid[87][4] c_mygrid[88][4] c_mygrid[89][4] & + c_mygrid[90][4] c_mygrid[91][4] c_mygrid[92][4] c_mygrid[93][4] c_mygrid[94][4] c_mygrid[95][4] c_mygrid[96][4] c_mygrid[97][4] c_mygrid[98][4] c_mygrid[99][4] & + c_mygrid[101][4] c_mygrid[102][4] c_mygrid[103][4] c_mygrid[104][4] c_mygrid[105][4] c_mygrid[106][4] c_mygrid[107][4] c_mygrid[108][4] c_mygrid[109][4] & + c_mygrid[110][4] c_mygrid[111][4] c_mygrid[112][4] c_mygrid[113][4] c_mygrid[114][4] c_mygrid[115][4] c_mygrid[116][4] c_mygrid[117][4] c_mygrid[118][4] c_mygrid[119][4] & + c_mygrid[120][4] c_mygrid[121][4] c_mygrid[122][4] c_mygrid[123][4] c_mygrid[124][4] c_mygrid[125][4] c_mygrid[126][4] c_mygrid[127][4] c_mygrid[128][4] c_mygrid[129][4] & + c_mygrid[130][4] c_mygrid[131][4] c_mygrid[132][4] c_mygrid[133][4] c_mygrid[134][4] c_mygrid[135][4] c_mygrid[136][4] c_mygrid[137][4] c_mygrid[138][4] c_mygrid[139][4] & + c_mygrid[140][4] c_mygrid[141][4] c_mygrid[142][4] c_mygrid[143][4] c_mygrid[144][4] c_mygrid[145][4] c_mygrid[146][4] c_mygrid[147][4] c_mygrid[148][4] c_mygrid[149][4] & + c_mygrid[150][4] c_mygrid[151][4] c_mygrid[152][4] c_mygrid[153][4] c_mygrid[154][4] c_mygrid[155][4] c_mygrid[156][4] c_mygrid[157][4] c_mygrid[158][4] c_mygrid[159][4] & + c_mygrid[160][4] c_mygrid[161][4] c_mygrid[162][4] c_mygrid[163][4] c_mygrid[164][4] c_mygrid[165][4] c_mygrid[166][4] c_mygrid[167][4] c_mygrid[168][4] c_mygrid[169][4] & + c_mygrid[170][4] c_mygrid[171][4] c_mygrid[172][4] c_mygrid[173][4] c_mygrid[174][4] c_mygrid[175][4] c_mygrid[176][4] c_mygrid[177][4] c_mygrid[178][4] c_mygrid[179][4] & + c_mygrid[180][4] c_mygrid[181][4] c_mygrid[182][4] c_mygrid[183][4] c_mygrid[184][4] c_mygrid[185][4] c_mygrid[186][4] c_mygrid[187][4] c_mygrid[188][4] c_mygrid[189][4] & + c_mygrid[190][4] c_mygrid[191][4] c_mygrid[192][4] c_mygrid[193][4] c_mygrid[194][4] c_mygrid[195][4] c_mygrid[196][4] c_mygrid[197][4] c_mygrid[198][4] c_mygrid[199][4] & + c_mygrid[201][4] c_mygrid[202][4] c_mygrid[203][4] c_mygrid[204][4] c_mygrid[205][4] c_mygrid[206][4] c_mygrid[207][4] c_mygrid[208][4] c_mygrid[209][4] & + c_mygrid[210][4] c_mygrid[211][4] c_mygrid[212][4] c_mygrid[213][4] c_mygrid[214][4] c_mygrid[215][4] c_mygrid[216][4] c_mygrid[217][4] c_mygrid[218][4] c_mygrid[219][4] & + c_mygrid[220][4] c_mygrid[221][4] c_mygrid[222][4] c_mygrid[223][4] c_mygrid[224][4] c_mygrid[225][4] c_mygrid[226][4] c_mygrid[227][4] c_mygrid[228][4] c_mygrid[229][4] & + c_mygrid[230][4] c_mygrid[231][4] c_mygrid[232][4] c_mygrid[233][4] c_mygrid[234][4] c_mygrid[235][4] c_mygrid[236][4] c_mygrid[237][4] c_mygrid[238][4] c_mygrid[239][4] & + c_mygrid[240][4] c_mygrid[241][4] c_mygrid[242][4] c_mygrid[243][4] c_mygrid[244][4] c_mygrid[245][4] c_mygrid[246][4] c_mygrid[247][4] c_mygrid[248][4] c_mygrid[249][4] & + c_mygrid[250][4] c_mygrid[251][4] c_mygrid[252][4] c_mygrid[253][4] c_mygrid[254][4] c_mygrid[255][4] c_mygrid[256][4] c_mygrid[257][4] c_mygrid[258][4] c_mygrid[259][4] & + c_mygrid[260][4] c_mygrid[261][4] c_mygrid[262][4] c_mygrid[263][4] c_mygrid[264][4] c_mygrid[265][4] c_mygrid[266][4] c_mygrid[267][4] c_mygrid[268][4] c_mygrid[269][4] & + c_mygrid[270][4] c_mygrid[271][4] c_mygrid[272][4] c_mygrid[273][4] c_mygrid[274][4] c_mygrid[275][4] c_mygrid[276][4] c_mygrid[277][4] c_mygrid[278][4] c_mygrid[279][4] & + c_mygrid[280][4] c_mygrid[281][4] c_mygrid[282][4] c_mygrid[283][4] c_mygrid[284][4] c_mygrid[285][4] c_mygrid[286][4] c_mygrid[287][4] c_mygrid[288][4] c_mygrid[289][4] & + c_mygrid[290][4] c_mygrid[291][4] c_mygrid[292][4] c_mygrid[293][4] c_mygrid[294][4] c_mygrid[295][4] c_mygrid[296][4] c_mygrid[297][4] c_mygrid[298][4] c_mygrid[299][4] & + c_mygrid[301][4] c_mygrid[302][4] c_mygrid[303][4] c_mygrid[304][4] c_mygrid[305][4] c_mygrid[306][4] c_mygrid[307][4] c_mygrid[308][4] c_mygrid[309][4] & + c_mygrid[310][4] c_mygrid[311][4] c_mygrid[312][4] c_mygrid[313][4] c_mygrid[314][4] c_mygrid[315][4] c_mygrid[316][4] c_mygrid[317][4] c_mygrid[318][4] c_mygrid[319][4] & + c_mygrid[320][4] c_mygrid[321][4] c_mygrid[322][4] c_mygrid[323][4] c_mygrid[324][4] c_mygrid[325][4] c_mygrid[326][4] c_mygrid[327][4] c_mygrid[328][4] c_mygrid[329][4] & + c_mygrid[330][4] c_mygrid[331][4] c_mygrid[332][4] c_mygrid[333][4] c_mygrid[334][4] c_mygrid[335][4] c_mygrid[336][4] c_mygrid[337][4] c_mygrid[338][4] c_mygrid[339][4] & + c_mygrid[340][4] c_mygrid[341][4] c_mygrid[342][4] c_mygrid[343][4] c_mygrid[344][4] c_mygrid[345][4] c_mygrid[346][4] c_mygrid[347][4] c_mygrid[348][4] c_mygrid[349][4] & + c_mygrid[350][4] c_mygrid[351][4] c_mygrid[352][4] c_mygrid[353][4] c_mygrid[354][4] c_mygrid[355][4] c_mygrid[356][4] c_mygrid[357][4] c_mygrid[358][4] c_mygrid[359][4] & + c_mygrid[360][4] c_mygrid[361][4] c_mygrid[362][4] c_mygrid[363][4] c_mygrid[364][4] c_mygrid[365][4] c_mygrid[366][4] c_mygrid[367][4] c_mygrid[368][4] c_mygrid[369][4] & + c_mygrid[370][4] c_mygrid[371][4] c_mygrid[372][4] c_mygrid[373][4] c_mygrid[374][4] c_mygrid[375][4] c_mygrid[376][4] c_mygrid[377][4] c_mygrid[378][4] c_mygrid[379][4] & + c_mygrid[380][4] c_mygrid[381][4] c_mygrid[382][4] c_mygrid[383][4] c_mygrid[384][4] thermo_modify norm yes diff --git a/src/SNAP/compute_sna_grid.cpp b/src/SNAP/compute_sna_grid.cpp index 26bd272d9f..d9912cdc47 100644 --- a/src/SNAP/compute_sna_grid.cpp +++ b/src/SNAP/compute_sna_grid.cpp @@ -120,7 +120,7 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : rmin0,switchflag,bzeroflag); ncoeff = snaptr->ncoeff; - size_array_cols = ncoeff; + size_array_cols = size_array_cols_base + ncoeff; if (quadraticflag) size_array_cols += (ncoeff*(ncoeff+1))/2; array_flag = 1; } @@ -237,8 +237,8 @@ void ComputeSNAGrid::compute_array() snaptr->compute_zi(); snaptr->compute_bi(); for (int icoeff = 0; icoeff < ncoeff; icoeff++) - sna[igrid][icoeff] = snaptr->blist[icoeff]; - // printf("igrid = %d %g %g %g %d B0 = %g\n",igrid,xtmp,ytmp,ztmp,ninside,sna[igrid][0]); + sna[igrid][size_array_cols_base+icoeff] = snaptr->blist[icoeff]; + // printf("igrid = %d %g %g %g %d B0 = %g\n",igrid,xtmp,ytmp,ztmp,ninside,sna[igrid][size_array_cols_base+0]); if (quadraticflag) { int ncount = ncoeff; for (int icoeff = 0; icoeff < ncoeff; icoeff++) { diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp index ec969fd01f..63e4a8c885 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -51,6 +51,7 @@ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : nargbase = iarg - iarg0; size_array_rows = nx*ny*nz; + size_array_cols_base = 3; } /* ---------------------------------------------------------------------- */ @@ -131,6 +132,7 @@ void ComputeGrid::setup() z0full = boxlo[2] - mz*delz; allocate(); + assign_coords_array(); } /* ---------------------------------------------------------------------- @@ -153,6 +155,26 @@ void ComputeGrid::igridfull2x(int igrid, double *x) } +/* ---------------------------------------------------------------------- + convert array index to box coords +------------------------------------------------------------------------- */ + +void ComputeGrid::iarray2x(int iarray, double *x) +{ + int iz = iarray / (nx*ny); + iarray -= iz * (nx*ny); + int iy = iarray / nx; + iarray -= iy * nx; + int ix = iarray; + + x[0] = ix*delx; + x[1] = iy*dely; + x[2] = iz*delz; + + if (triclinic) domain->lamda2x(x, x); + +} + /* ---------------------------------------------------------------------- copy local grid to global array ------------------------------------------------------------------------- */ @@ -162,11 +184,26 @@ void ComputeGrid::copy_local_grid() int igridfull; for (int iarray = 0; iarray < size_array_rows; iarray++) { igridfull = iarray2igridfull(iarray); - for (int icol = 0; icol < size_array_cols; icol++) + for (int icol = size_array_cols_base; icol < size_array_cols; icol++) array[iarray][icol] = gridfull[igridfull][icol]; } } +/* ---------------------------------------------------------------------- + copy coords to global array +------------------------------------------------------------------------- */ + +void ComputeGrid::assign_coords_array() +{ + double x[3]; + for (int iarray = 0; iarray < size_array_rows; iarray++) { + iarray2x(iarray,x); + array[iarray][0] = x[0]; + array[iarray][1] = x[1]; + array[iarray][2] = x[2]; + } +} + /* ---------------------------------------------------------------------- gather global array from full grid ------------------------------------------------------------------------- */ @@ -181,7 +218,7 @@ void ComputeGrid::gather_global_array() // inefficient, should exploit shared ix structure iarray = igridfull2iarray(igrid); - for (int icol = 0; icol < size_array_cols; icol++) + for (int icol = size_array_cols_base; icol < size_array_cols; icol++) array[iarray][icol] += gridfull[igrid][icol]; } } diff --git a/src/compute_grid.h b/src/compute_grid.h index bc757c37bc..1468f9ad91 100644 --- a/src/compute_grid.h +++ b/src/compute_grid.h @@ -33,7 +33,7 @@ class ComputeGrid : public Compute { int nx, ny, nz; // grid dimensions int nxfull, nyfull, nzfull; // grid dimensions with ghost points int nxyfull; // nx_full*ny_full - int ngridfull; // number of full grid points + int ngridfull; // number of full grid points double **gridfull; // full grid points int mx, my, mz; // cutmax stencil dimensions int triclinic; // triclinic flag @@ -43,10 +43,13 @@ class ComputeGrid : public Compute { double x0full, y0full, z0full; // origin of full grid int nargbase; // number of base class args double cutmax; // largest cutoff distance + int size_array_cols_base; // number of columns used for coords, etc. virtual void allocate(); void igridfull2x(int, double*); // convert full grid point to coord + void iarray2x(int, double*); // convert array point to coord void gather_global_array(); // gather global array from full grid void copy_local_grid(); // copy local grid to global array + void assign_coords_array(); // assign coords to global array int igridfull2iarray(int); // convert full grid index to compute array index int iarray2igridfull(int); // convert compute array index to full grid index From 59e3b4c5ba8b3de1d7f183a1cecb793b14886b99 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Thu, 24 Oct 2019 16:12:23 -0600 Subject: [PATCH 06/19] Fixed bug in lammps.py --- python/lammps.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/python/lammps.py b/python/lammps.py index 36cf2d2fdd..f23268cd60 100644 --- a/python/lammps.py +++ b/python/lammps.py @@ -399,14 +399,9 @@ class lammps(object): ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) return ptr if type == 2: - if style == 0: - self.lib.lammps_extract_compute.restype = POINTER(c_int) - ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) - return ptr[0] - else: - self.lib.lammps_extract_compute.restype = POINTER(POINTER(c_double)) - ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) - return ptr + self.lib.lammps_extract_compute.restype = POINTER(POINTER(c_double)) + ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type) + return ptr return None # extract fix info From 2fa9e5fefb26f3e4540ed7dc6e4e1563e49432bc Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Thu, 24 Oct 2019 19:48:41 -0600 Subject: [PATCH 07/19] Completed brute force parallel implementation using MPI_Allreduce() --- src/SNAP/compute_sna_grid.cpp | 60 +++------- src/SNAP/compute_sna_grid.h | 4 - src/compute_grid.cpp | 220 ++++++++++------------------------ src/compute_grid.h | 26 ++-- 4 files changed, 93 insertions(+), 217 deletions(-) diff --git a/src/SNAP/compute_sna_grid.cpp b/src/SNAP/compute_sna_grid.cpp index d9912cdc47..a9f3b6c3fb 100644 --- a/src/SNAP/compute_sna_grid.cpp +++ b/src/SNAP/compute_sna_grid.cpp @@ -31,7 +31,7 @@ using namespace LAMMPS_NS; ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : - ComputeGrid(lmp, narg, arg), cutsq(NULL), list(NULL), sna(NULL), + ComputeGrid(lmp, narg, arg), cutsq(NULL), sna(NULL), radelem(NULL), wjelem(NULL) { double rmin0, rfac0; @@ -120,8 +120,9 @@ ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : rmin0,switchflag,bzeroflag); ncoeff = snaptr->ncoeff; - size_array_cols = size_array_cols_base + ncoeff; - if (quadraticflag) size_array_cols += (ncoeff*(ncoeff+1))/2; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff*(ncoeff+1))/2; + size_array_cols = size_array_cols_base + nvalues; array_flag = 1; } @@ -165,26 +166,10 @@ void ComputeSNAGrid::init() /* ---------------------------------------------------------------------- */ -void ComputeSNAGrid::init_list(int /*id*/, NeighList *ptr) -{ - list = ptr; -} - -/* ---------------------------------------------------------------------- */ - void ComputeSNAGrid::compute_array() { invoked_array = update->ntimestep; -// // invoke full neighbor list (will copy or build if necessary) - -// neighbor->build_one(list); - -// const int inum = list->inum; -// const int* const ilist = list->ilist; -// const int* const numneigh = list->numneigh; -// int** const firstneigh = list->firstneigh; - int * const type = atom->type; // compute sna for each gridpoint @@ -197,13 +182,12 @@ void ComputeSNAGrid::compute_array() snaptr->grow_rij(ntotal); - printf("ngridfull = %d\n",ngridfull); - for (int igrid = 0; igrid < ngridfull; igrid++) { - double rtmp[3]; - igridfull2x(igrid, rtmp); - const double xtmp = rtmp[0]; - const double ytmp = rtmp[1]; - const double ztmp = rtmp[2]; + printf("ngrid = %d\n",ngrid); + for (int igrid = 0; igrid < ngrid; igrid++) { + if (!grid_local[igrid]) continue; + const double xtmp = grid[igrid][0]; + const double ytmp = grid[igrid][1]; + const double ztmp = grid[igrid][2]; // rij[][3] = displacements between atom I and those neighbors // inside = indices of neighbors of I within cutoff @@ -237,7 +221,7 @@ void ComputeSNAGrid::compute_array() snaptr->compute_zi(); snaptr->compute_bi(); for (int icoeff = 0; icoeff < ncoeff; icoeff++) - sna[igrid][size_array_cols_base+icoeff] = snaptr->blist[icoeff]; + grid[igrid][size_array_cols_base+icoeff] = snaptr->blist[icoeff]; // printf("igrid = %d %g %g %g %d B0 = %g\n",igrid,xtmp,ytmp,ztmp,ninside,sna[igrid][size_array_cols_base+0]); if (quadraticflag) { int ncount = ncoeff; @@ -246,27 +230,16 @@ void ComputeSNAGrid::compute_array() // diagonal element of quadratic matrix - sna[igrid][ncount++] = 0.5*bi*bi; + grid[igrid][size_array_cols_base+ncount++] = 0.5*bi*bi; // upper-triangular elements of quadratic matrix for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) - sna[igrid][ncount++] = bi*snaptr->blist[jcoeff]; + grid[igrid][size_array_cols_base+ncount++] = bi*snaptr->blist[jcoeff]; } } } - // gather_global_array(); - copy_local_grid(); -} - -/* ---------------------------------------------------------------------- - allocate array in base class and then set up pointers -------------------------------------------------------------------------- */ - -void ComputeSNAGrid::allocate() -{ - ComputeGrid::allocate(); - sna = gridfull; + MPI_Allreduce(&grid[0][0],&gridall[0][0],ngrid*size_array_cols,MPI_DOUBLE,MPI_SUM,world); } /* ---------------------------------------------------------------------- @@ -275,9 +248,8 @@ void ComputeSNAGrid::allocate() double ComputeSNAGrid::memory_usage() { - double bytes = size_array_rows*size_array_cols * sizeof(double); // grid - bytes += snaptr->memory_usage(); // SNA object + double nbytes = snaptr->memory_usage(); // SNA object - return bytes; + return nbytes; } diff --git a/src/SNAP/compute_sna_grid.h b/src/SNAP/compute_sna_grid.h index 2bd18a915e..0242a2962b 100644 --- a/src/SNAP/compute_sna_grid.h +++ b/src/SNAP/compute_sna_grid.h @@ -29,16 +29,12 @@ class ComputeSNAGrid : public ComputeGrid { ComputeSNAGrid(class LAMMPS *, int, char **); ~ComputeSNAGrid(); void init(); - void init_list(int, class NeighList *); void compute_array(); double memory_usage(); - void allocate(); private: - int nmax; int ncoeff; double **cutsq; - class NeighList *list; double **sna; double rcutfac; double *radelem; diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp index 63e4a8c885..2f67d22182 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -50,7 +50,7 @@ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : nargbase = iarg - iarg0; - size_array_rows = nx*ny*nz; + size_array_rows = ngrid = nx*ny*nz; size_array_cols_base = 3; } @@ -58,6 +58,8 @@ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : ComputeGrid::~ComputeGrid() { + memory->destroy(grid); + memory->destroy(grid_local); } /* ---------------------------------------------------------------------- */ @@ -78,9 +80,13 @@ void ComputeGrid::setup() if (triclinic == 0) { prd = domain->prd; boxlo = domain->boxlo; + sublo = domain->sublo; + subhi = domain->subhi; } else { prd = domain->prd_lamda; boxlo = domain->boxlo_lamda; + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; } double xprd = prd[0]; @@ -95,78 +101,23 @@ void ComputeGrid::setup() dely = 1.0/delyinv; delz = 1.0/delzinv; - // sufficient conditions for stencil bounding rcut - - // require |delz*mz|^2 <= rcut^2 - // require |dely*my|^2 <= rcut^2 + |delyz*mz_max|^2 - // require |delx*mx|^2 <= rcut^2 + |delxz*mz_max|^2 + |delxy*my_max|^2 - - double delxy = domain->xy/ny; - double delxz = domain->xz/nz; - double delyz = domain->yz/nz; - - if (!triclinic) { - mz = cutmax*delzinv + 1; - my = sqrt(cutmax*cutmax + pow(delyz*mz,2))*delyinv + 1; - mx = sqrt(cutmax*cutmax + pow(delxz*mz,2) - + pow(delxy*my,2))*delxinv + 1; - } else { - double delxinvtmp = nx/domain->xprd; - double delyinvtmp = ny/domain->yprd; - double delzinvtmp = nz/domain->zprd; - mz = cutmax*delzinvtmp + 1; - my = sqrt(cutmax*cutmax + pow(delyz*mz,2))*delyinvtmp + 1; - mx = sqrt(cutmax*cutmax + pow(delxz*mz,2) - + pow(delxy*my,2))*delxinvtmp + 1; - } - - // size global grid to accomodate periodic interactions - - nxfull = nx + 2*mx; - nyfull = ny + 2*my; - nzfull = nz + 2*mz; - nxyfull = nxfull * nyfull; - - x0full = boxlo[0] - mx*delx; - y0full = boxlo[1] - my*dely; - z0full = boxlo[2] - mz*delz; - allocate(); - assign_coords_array(); + assign_grid_coords(); + assign_grid_local(); } /* ---------------------------------------------------------------------- - convert grid index to box coords + convert global array index to box coords ------------------------------------------------------------------------- */ -void ComputeGrid::igridfull2x(int igrid, double *x) +void ComputeGrid::grid2x(int igrid, double *x) { - int iz = igrid / nxyfull; - igrid -= iz * nxyfull; - int iy = igrid / nxfull; - igrid -= iy * nxfull; + int iz = igrid / (nx*ny); + igrid -= iz * (nx*ny); + int iy = igrid / nx; + igrid -= iy * nx; int ix = igrid; - x[0] = x0full+ix*delx; - x[1] = y0full+iy*dely; - x[2] = z0full+iz*delz; - - if (triclinic) domain->lamda2x(x, x); - -} - -/* ---------------------------------------------------------------------- - convert array index to box coords -------------------------------------------------------------------------- */ - -void ComputeGrid::iarray2x(int iarray, double *x) -{ - int iz = iarray / (nx*ny); - iarray -= iz * (nx*ny); - int iy = iarray / nx; - iarray -= iy * nx; - int ix = iarray; - x[0] = ix*delx; x[1] = iy*dely; x[2] = iz*delz; @@ -176,16 +127,43 @@ void ComputeGrid::iarray2x(int iarray, double *x) } /* ---------------------------------------------------------------------- - copy local grid to global array + check if grid point is local ------------------------------------------------------------------------- */ -void ComputeGrid::copy_local_grid() +int ComputeGrid::check_grid_local(int igrid) { - int igridfull; - for (int iarray = 0; iarray < size_array_rows; iarray++) { - igridfull = iarray2igridfull(iarray); - for (int icol = size_array_cols_base; icol < size_array_cols; icol++) - array[iarray][icol] = gridfull[igridfull][icol]; + double x[3]; + + 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; + + int islocal = + x[0] >= sublo[0] && x[0] < subhi[0] && + x[1] >= sublo[1] && x[1] < subhi[1] && + x[2] >= sublo[2] && x[2] < subhi[2]; + + return islocal; +} + +/* ---------------------------------------------------------------------- + copy coords to global array +------------------------------------------------------------------------- */ + +void ComputeGrid::assign_grid_coords() +{ + double x[3]; + for (int igrid = 0; igrid < ngrid; igrid++) { + grid2x(igrid,x); + grid[igrid][0] = x[0]; + grid[igrid][1] = x[1]; + grid[igrid][2] = x[2]; } } @@ -193,101 +171,33 @@ void ComputeGrid::copy_local_grid() copy coords to global array ------------------------------------------------------------------------- */ -void ComputeGrid::assign_coords_array() +void ComputeGrid::assign_grid_local() { double x[3]; - for (int iarray = 0; iarray < size_array_rows; iarray++) { - iarray2x(iarray,x); - array[iarray][0] = x[0]; - array[iarray][1] = x[1]; - array[iarray][2] = x[2]; + for (int igrid = 0; igrid < ngrid; igrid++) { + if (check_grid_local(igrid)) + grid_local[igrid] = 1; + else { + grid_local[igrid] = 0; + memset(grid[igrid],0,size_array_cols); + } } } -/* ---------------------------------------------------------------------- - gather global array from full grid -------------------------------------------------------------------------- */ - -void ComputeGrid::gather_global_array() -{ - int iarray; - memset(&array[0][0],0,size_array_rows*size_array_cols*sizeof(double)); - - for (int igrid = 0; igrid < ngridfull; igrid++) { - - // inefficient, should exploit shared ix structure - - iarray = igridfull2iarray(igrid); - for (int icol = size_array_cols_base; icol < size_array_cols; icol++) - array[iarray][icol] += gridfull[igrid][icol]; - } -} - -/* ---------------------------------------------------------------------- - convert full grid index to compute array index - inefficient, should exploit shared ix structure -------------------------------------------------------------------------- */ - -int ComputeGrid::igridfull2iarray(int igrid) -{ - int iz = igrid / nxyfull; - igrid -= iz*nxyfull; - int iy = igrid / nxfull; - igrid -= iy*nxfull; - int ix = igrid; - - ix -= mx; - iy -= my; - iz -= mz; - - while (ix < 0) ix += nx; - while (iy < 0) iy += ny; - while (iz < 0) iz += nz; - - while (ix >= nx) ix -= nx; - while (iy >= ny) iy -= ny; - while (iz >= nz) iz -= nz; - - int iarray = (iz * ny + iy) * nx + ix; - - return iarray; -} - -/* ---------------------------------------------------------------------- - convert compute array index to full grid index - inefficient, should exploit shared ix structure -------------------------------------------------------------------------- */ - -int ComputeGrid::iarray2igridfull(int iarray) -{ - int iz = iarray / (nx*ny); - iarray -= iz*(nx*ny); - int iy = iarray / nx; - iarray -= iy*nx; - int ix = iarray; - - ix += mx; - iy += my; - iz += mz; - - int igrid = (iz * nyfull + iy) * nxfull + ix; - - return igrid; -} - /* ---------------------------------------------------------------------- free and reallocate arrays ------------------------------------------------------------------------- */ void ComputeGrid::allocate() { - ngridfull = nxfull*nyfull*nzfull; - // grow global array if necessary - memory->destroy(array); - memory->create(array,size_array_rows,size_array_cols,"sna/grid:array"); - memory->create(gridfull,ngridfull,size_array_cols,"sna/grid:gridfull"); + memory->destroy(grid); + memory->destroy(grid_local); + memory->create(grid,size_array_rows,size_array_cols,"grid:grid"); + memory->create(gridall,size_array_rows,size_array_cols,"grid:gridall"); + memory->create(grid_local,size_array_rows,"grid:grid_local"); + array = gridall; } /* ---------------------------------------------------------------------- memory usage of local data @@ -295,6 +205,8 @@ void ComputeGrid::allocate() double ComputeGrid::memory_usage() { - int nbytes = 0; + double nbytes = size_array_rows*size_array_cols * + sizeof(double); // grid + nbytes += size_array_rows*sizeof(int); // grid_local return nbytes; } diff --git a/src/compute_grid.h b/src/compute_grid.h index 1468f9ad91..a1c938db2e 100644 --- a/src/compute_grid.h +++ b/src/compute_grid.h @@ -31,28 +31,24 @@ class ComputeGrid : public Compute { protected: int nx, ny, nz; // grid dimensions - int nxfull, nyfull, nzfull; // grid dimensions with ghost points - int nxyfull; // nx_full*ny_full - int ngridfull; // number of full grid points - double **gridfull; // full grid points - int mx, my, mz; // cutmax stencil dimensions + int ngrid; // number of grid points + int nvalues; // number of values per grid point + double **grid; // global grid + double **gridall; // global grid summed over procs 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 - double x0full, y0full, z0full; // origin of full grid int nargbase; // number of base class args double cutmax; // largest cutoff distance int size_array_cols_base; // number of columns used for coords, etc. - virtual void allocate(); - void igridfull2x(int, double*); // convert full grid point to coord - void iarray2x(int, double*); // convert array point to coord - void gather_global_array(); // gather global array from full grid - void copy_local_grid(); // copy local grid to global array - void assign_coords_array(); // assign coords to global array - int igridfull2iarray(int); // convert full grid index to compute array index - int iarray2igridfull(int); // convert compute array index to full grid index - + int *grid_local; // local flag for each grid point + void allocate(); + void grid2x(int, double*); // convert grid point to coord + void assign_grid_coords(); // assign coords for grid + void assign_grid_local(); // set local flag for each grid point + int check_grid_local(int); // check if grid point is local private: }; From 651e9c63973809d1026faac942f84998dfe7ac74 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Tue, 19 May 2020 11:01:31 -0600 Subject: [PATCH 08/19] Initialized arrays to NULL --- src/compute_grid.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp index 2f67d22182..011dafa030 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -27,7 +27,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) + Compute(lmp, narg, arg), grid(NULL), grid_local(NULL) { if (narg < 6) error->all(FLERR,"Illegal compute grid command"); From 4a261f396128a0d19beba082bd7f4ab6c54a1803 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Sat, 23 May 2020 10:18:16 -0600 Subject: [PATCH 09/19] Added python example (!) and turned on switching function --- examples/snap/grid.py | 104 +++++++++++++++++++++++++++++++++++ examples/snap/in.grid | 9 ++- examples/snap/in.grid.python | 61 ++++++++++++++++++++ examples/snap/in.grid.tri | 82 +++++++++++++-------------- 4 files changed, 214 insertions(+), 42 deletions(-) create mode 100755 examples/snap/grid.py create mode 100644 examples/snap/in.grid.python diff --git a/examples/snap/grid.py b/examples/snap/grid.py new file mode 100755 index 0000000000..9fe88fd1b4 --- /dev/null +++ b/examples/snap/grid.py @@ -0,0 +1,104 @@ +#!/Users/athomps/miniconda3/bin/python3.7 +# +# An example of SNAP grid from LAMMPS Python interface +# +# https://lammps.sandia.gov/doc/Python_library.html + + +from lammps import lammps +import lammps_utils + +# define command line input variables + +ngridx = 2 +ngridy = 3 +ngridz = 4 +twojmax = 2 + +lmp_cmdargs = ["-echo","screen"] +lmp_cmdargs = lammps_utils.set_cmdlinevars(lmp_cmdargs, + { + "ngridx":ngridx, + "ngridy":ngridy, + "ngridz":ngridz, + "twojmax":twojmax + } + ) + +# launch LAMMPS instance + +lmp = lammps(cmdargs=lmp_cmdargs) + +# run LAMMPS input script + +lmp.file("in.grid.python") + +# get quantities from LAMMPS + +num_atoms = lmp.get_natoms() + +# set things not accessible from LAMMPS + +# first 3 cols are x, y, z, coords + +ncols0 = 3 + +# analytical relation + +ncoeff = (twojmax+2)*(twojmax+3)*(twojmax+4) +ncoeff = ncoeff // 24 # integer division +ncols = ncols0+ncoeff + +# get B_0 at position (0,0,0) in 4 different ways + +# 1. from comute sna/atom + +bptr = lmp.extract_compute("b", 1, 2) # 1 = per-atom data, 2 = array +print("b = ",bptr[0][0]) + +# 2. from compute sna/grid + +bgridptr = lmp.extract_compute("bgrid", 0, 2) # 0 = style global, 2 = type array +print("bgrid = ",bgridptr[0][3]) + +# 3. from Numpy array pointing to sna/atom array + +bptr_np = lammps_utils.extract_compute_np(lmp,"b",1,2,(num_atoms,ncoeff)) +print("b_np = ",bptr_np[0][0]) + +# 4. from Numpy array pointing to sna/grid array + +bgridptr_np = lammps_utils.extract_compute_np(lmp,"bgrid",0,2,(ngridz,ngridy,ngridx,ncols)) +print("bgrid_np = ",bgridptr_np[0][0][0][ncols0+0]) + +# print out the LAMMPS array to a file + +outfile = open("bgrid.dat",'w') +igrid = 0 +for iz in range(ngridz): + for iy in range(ngridy): + for ix in range(ngridx): + outfile.write("x, y, z = %g %g %g\n" % + (bgridptr[igrid][0], + bgridptr[igrid][1], + bgridptr[igrid][2])) + for icoeff in range(ncoeff): + outfile.write("%g " % bgridptr[igrid][ncols0+icoeff]) + outfile.write("\n") + igrid += 1 +outfile.close() + +# print out the Numpy array to a file + +outfile = open("bgrid_np.dat",'w') +for iz in range(ngridz): + for iy in range(ngridy): + for ix in range(ngridx): + outfile.write("x, y, z = %g %g %g\n" % + (bgridptr_np[iz][iy][ix][0], + bgridptr_np[iz][iy][ix][1], + bgridptr_np[iz][iy][ix][2])) + for icoeff in range(ncoeff): + outfile.write("%g " % bgridptr_np[iz][iy][ix][ncols0+icoeff]) + outfile.write("\n") +outfile.close() diff --git a/examples/snap/in.grid b/examples/snap/in.grid index 5293558fab..2702fccad1 100644 --- a/examples/snap/in.grid +++ b/examples/snap/in.grid @@ -46,7 +46,7 @@ variable wj equal 1 variable radelem equal 0.5 variable bzero equal 0 variable quad equal 0 -variable switch equal 0 +variable switch equal 1 compute b all sna/atom & ${rcutfac} ${rfac0} ${twojmax} ${radelem} & @@ -60,11 +60,16 @@ compute mygrid all sna/grid grid ${ngrid} ${ngrid} ${ngrid} & # define output +# mygrid is ngrid by (3+nbis) = 8x8 thermo_style custom step temp ke pe vol & c_mygrid[1][1] c_mygrid[2][1] c_mygrid[3][1] c_mygrid[4][1] c_mygrid[5][1] c_mygrid[6][1] c_mygrid[7][1] c_mygrid[8][1] & c_mygrid[1][2] c_mygrid[2][2] c_mygrid[3][2] c_mygrid[4][2] c_mygrid[5][2] c_mygrid[6][2] c_mygrid[7][2] c_mygrid[8][2] & c_mygrid[1][3] c_mygrid[2][3] c_mygrid[3][3] c_mygrid[4][3] c_mygrid[5][3] c_mygrid[6][3] c_mygrid[7][3] c_mygrid[8][3] & - c_mygrid[1][4] c_mygrid[2][4] c_mygrid[3][4] c_mygrid[4][4] c_mygrid[5][4] c_mygrid[6][4] c_mygrid[7][4] c_mygrid[8][4] + c_mygrid[1][4] c_mygrid[2][4] c_mygrid[3][4] c_mygrid[4][4] c_mygrid[5][4] c_mygrid[6][4] c_mygrid[7][4] c_mygrid[8][4] & + c_mygrid[1][5] c_mygrid[2][5] c_mygrid[3][5] c_mygrid[4][5] c_mygrid[5][5] c_mygrid[6][5] c_mygrid[7][5] c_mygrid[8][5] & + c_mygrid[1][6] c_mygrid[2][6] c_mygrid[3][6] c_mygrid[4][6] c_mygrid[5][6] c_mygrid[6][6] c_mygrid[7][6] c_mygrid[8][6] & + c_mygrid[1][7] c_mygrid[2][7] c_mygrid[3][7] c_mygrid[4][7] c_mygrid[5][7] c_mygrid[6][7] c_mygrid[7][7] c_mygrid[8][7] & + c_mygrid[1][8] c_mygrid[2][8] c_mygrid[3][8] c_mygrid[4][8] c_mygrid[5][8] c_mygrid[6][8] c_mygrid[7][8] c_mygrid[8][8] thermo_modify norm yes dump mydump_b all custom 1000 dump_b id c_b[*] diff --git a/examples/snap/in.grid.python b/examples/snap/in.grid.python new file mode 100644 index 0000000000..13fe6d0638 --- /dev/null +++ b/examples/snap/in.grid.python @@ -0,0 +1,61 @@ +# Demonstrate bispectrum per-atom and grid computes + +# Invoked from grid.py +# pass in values for ngridx, ngridy, ngridz, twojmax + +# Initialize simulation + +variable nsteps equal 0 +variable nrep equal 1 +variable a equal 3.316 + +units metal + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrep} +variable ny equal ${nrep} +variable nz equal ${nrep} + +boundary p p p + +lattice custom ${a} a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0 0 0 basis 0.5 0.5 0.5 + +region box block 0 ${nx} 0 ${ny} 0 ${nz} +create_box 1 box +create_atoms 1 box + +mass 1 180.88 + +# define grid compute and atom compute + +group snapgroup type 1 +variable rcutfac equal 4.67637 +variable rfac0 equal 0.99363 +variable rmin0 equal 0 +variable wj equal 1 +variable radelem equal 0.5 +variable bzero equal 0 +variable quad equal 0 +variable switch equal 1 + +compute b all sna/atom ${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} bzeroflag ${bzero} quadraticflag ${quad} switchflag ${switch} + +compute bgrid all sna/grid grid ${ngridx} ${ngridy} ${ngridz} ${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} bzeroflag ${bzero} quadraticflag ${quad} switchflag ${switch} + +# create dummy potential for neighbor list + +variable rcutneigh equal 2.0*${rcutfac}*${radelem} +pair_style zero ${rcutneigh} +pair_coeff * * + +# define output + +thermo_style custom step temp ke pe vol c_bgrid[1][1] +thermo_modify norm yes + +dump mydump_b all custom 1000 dump_b id c_b[*] + +# run + +run 0 diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri index 3c471d9b66..c984799e09 100644 --- a/examples/snap/in.grid.tri +++ b/examples/snap/in.grid.tri @@ -51,7 +51,7 @@ variable wj equal 1 variable radelem equal 0.5 variable bzero equal 0 variable quad equal 0 -variable switch equal 0 +variable switch equal 1 compute b all sna/atom & ${rcutfac} ${rfac0} ${twojmax} ${radelem} & @@ -65,46 +65,48 @@ compute mygrid all sna/grid grid ${ngridx} ${ngridy} ${ngridz} & # define output +# mygrid is ngrid by (3+nbis) = 384x8 + thermo_style custom step temp ke pe vol & - c_mygrid[1][4] c_mygrid[2][4] c_mygrid[3][4] c_mygrid[4][4] c_mygrid[5][4] c_mygrid[6][4] c_mygrid[7][4] c_mygrid[8][4] c_mygrid[9][4] & - c_mygrid[10][4] c_mygrid[11][4] c_mygrid[12][4] c_mygrid[13][4] c_mygrid[14][4] c_mygrid[15][4] c_mygrid[16][4] c_mygrid[17][4] c_mygrid[18][4] c_mygrid[19][4] & - c_mygrid[20][4] c_mygrid[21][4] c_mygrid[22][4] c_mygrid[23][4] c_mygrid[24][4] c_mygrid[25][4] c_mygrid[26][4] c_mygrid[27][4] c_mygrid[28][4] c_mygrid[29][4] & - c_mygrid[30][4] c_mygrid[31][4] c_mygrid[32][4] c_mygrid[33][4] c_mygrid[34][4] c_mygrid[35][4] c_mygrid[36][4] c_mygrid[37][4] c_mygrid[38][4] c_mygrid[39][4] & - c_mygrid[40][4] c_mygrid[41][4] c_mygrid[42][4] c_mygrid[43][4] c_mygrid[44][4] c_mygrid[45][4] c_mygrid[46][4] c_mygrid[47][4] c_mygrid[48][4] c_mygrid[49][4] & - c_mygrid[50][4] c_mygrid[51][4] c_mygrid[52][4] c_mygrid[53][4] c_mygrid[54][4] c_mygrid[55][4] c_mygrid[56][4] c_mygrid[57][4] c_mygrid[58][4] c_mygrid[59][4] & - c_mygrid[60][4] c_mygrid[61][4] c_mygrid[62][4] c_mygrid[63][4] c_mygrid[64][4] c_mygrid[65][4] c_mygrid[66][4] c_mygrid[67][4] c_mygrid[68][4] c_mygrid[69][4] & - c_mygrid[70][4] c_mygrid[71][4] c_mygrid[72][4] c_mygrid[73][4] c_mygrid[74][4] c_mygrid[75][4] c_mygrid[76][4] c_mygrid[77][4] c_mygrid[78][4] c_mygrid[79][4] & - c_mygrid[80][4] c_mygrid[81][4] c_mygrid[82][4] c_mygrid[83][4] c_mygrid[84][4] c_mygrid[85][4] c_mygrid[86][4] c_mygrid[87][4] c_mygrid[88][4] c_mygrid[89][4] & - c_mygrid[90][4] c_mygrid[91][4] c_mygrid[92][4] c_mygrid[93][4] c_mygrid[94][4] c_mygrid[95][4] c_mygrid[96][4] c_mygrid[97][4] c_mygrid[98][4] c_mygrid[99][4] & - c_mygrid[101][4] c_mygrid[102][4] c_mygrid[103][4] c_mygrid[104][4] c_mygrid[105][4] c_mygrid[106][4] c_mygrid[107][4] c_mygrid[108][4] c_mygrid[109][4] & - c_mygrid[110][4] c_mygrid[111][4] c_mygrid[112][4] c_mygrid[113][4] c_mygrid[114][4] c_mygrid[115][4] c_mygrid[116][4] c_mygrid[117][4] c_mygrid[118][4] c_mygrid[119][4] & - c_mygrid[120][4] c_mygrid[121][4] c_mygrid[122][4] c_mygrid[123][4] c_mygrid[124][4] c_mygrid[125][4] c_mygrid[126][4] c_mygrid[127][4] c_mygrid[128][4] c_mygrid[129][4] & - c_mygrid[130][4] c_mygrid[131][4] c_mygrid[132][4] c_mygrid[133][4] c_mygrid[134][4] c_mygrid[135][4] c_mygrid[136][4] c_mygrid[137][4] c_mygrid[138][4] c_mygrid[139][4] & - c_mygrid[140][4] c_mygrid[141][4] c_mygrid[142][4] c_mygrid[143][4] c_mygrid[144][4] c_mygrid[145][4] c_mygrid[146][4] c_mygrid[147][4] c_mygrid[148][4] c_mygrid[149][4] & - c_mygrid[150][4] c_mygrid[151][4] c_mygrid[152][4] c_mygrid[153][4] c_mygrid[154][4] c_mygrid[155][4] c_mygrid[156][4] c_mygrid[157][4] c_mygrid[158][4] c_mygrid[159][4] & - c_mygrid[160][4] c_mygrid[161][4] c_mygrid[162][4] c_mygrid[163][4] c_mygrid[164][4] c_mygrid[165][4] c_mygrid[166][4] c_mygrid[167][4] c_mygrid[168][4] c_mygrid[169][4] & - c_mygrid[170][4] c_mygrid[171][4] c_mygrid[172][4] c_mygrid[173][4] c_mygrid[174][4] c_mygrid[175][4] c_mygrid[176][4] c_mygrid[177][4] c_mygrid[178][4] c_mygrid[179][4] & - c_mygrid[180][4] c_mygrid[181][4] c_mygrid[182][4] c_mygrid[183][4] c_mygrid[184][4] c_mygrid[185][4] c_mygrid[186][4] c_mygrid[187][4] c_mygrid[188][4] c_mygrid[189][4] & - c_mygrid[190][4] c_mygrid[191][4] c_mygrid[192][4] c_mygrid[193][4] c_mygrid[194][4] c_mygrid[195][4] c_mygrid[196][4] c_mygrid[197][4] c_mygrid[198][4] c_mygrid[199][4] & - c_mygrid[201][4] c_mygrid[202][4] c_mygrid[203][4] c_mygrid[204][4] c_mygrid[205][4] c_mygrid[206][4] c_mygrid[207][4] c_mygrid[208][4] c_mygrid[209][4] & - c_mygrid[210][4] c_mygrid[211][4] c_mygrid[212][4] c_mygrid[213][4] c_mygrid[214][4] c_mygrid[215][4] c_mygrid[216][4] c_mygrid[217][4] c_mygrid[218][4] c_mygrid[219][4] & - c_mygrid[220][4] c_mygrid[221][4] c_mygrid[222][4] c_mygrid[223][4] c_mygrid[224][4] c_mygrid[225][4] c_mygrid[226][4] c_mygrid[227][4] c_mygrid[228][4] c_mygrid[229][4] & - c_mygrid[230][4] c_mygrid[231][4] c_mygrid[232][4] c_mygrid[233][4] c_mygrid[234][4] c_mygrid[235][4] c_mygrid[236][4] c_mygrid[237][4] c_mygrid[238][4] c_mygrid[239][4] & - c_mygrid[240][4] c_mygrid[241][4] c_mygrid[242][4] c_mygrid[243][4] c_mygrid[244][4] c_mygrid[245][4] c_mygrid[246][4] c_mygrid[247][4] c_mygrid[248][4] c_mygrid[249][4] & - c_mygrid[250][4] c_mygrid[251][4] c_mygrid[252][4] c_mygrid[253][4] c_mygrid[254][4] c_mygrid[255][4] c_mygrid[256][4] c_mygrid[257][4] c_mygrid[258][4] c_mygrid[259][4] & - c_mygrid[260][4] c_mygrid[261][4] c_mygrid[262][4] c_mygrid[263][4] c_mygrid[264][4] c_mygrid[265][4] c_mygrid[266][4] c_mygrid[267][4] c_mygrid[268][4] c_mygrid[269][4] & - c_mygrid[270][4] c_mygrid[271][4] c_mygrid[272][4] c_mygrid[273][4] c_mygrid[274][4] c_mygrid[275][4] c_mygrid[276][4] c_mygrid[277][4] c_mygrid[278][4] c_mygrid[279][4] & - c_mygrid[280][4] c_mygrid[281][4] c_mygrid[282][4] c_mygrid[283][4] c_mygrid[284][4] c_mygrid[285][4] c_mygrid[286][4] c_mygrid[287][4] c_mygrid[288][4] c_mygrid[289][4] & - c_mygrid[290][4] c_mygrid[291][4] c_mygrid[292][4] c_mygrid[293][4] c_mygrid[294][4] c_mygrid[295][4] c_mygrid[296][4] c_mygrid[297][4] c_mygrid[298][4] c_mygrid[299][4] & - c_mygrid[301][4] c_mygrid[302][4] c_mygrid[303][4] c_mygrid[304][4] c_mygrid[305][4] c_mygrid[306][4] c_mygrid[307][4] c_mygrid[308][4] c_mygrid[309][4] & - c_mygrid[310][4] c_mygrid[311][4] c_mygrid[312][4] c_mygrid[313][4] c_mygrid[314][4] c_mygrid[315][4] c_mygrid[316][4] c_mygrid[317][4] c_mygrid[318][4] c_mygrid[319][4] & - c_mygrid[320][4] c_mygrid[321][4] c_mygrid[322][4] c_mygrid[323][4] c_mygrid[324][4] c_mygrid[325][4] c_mygrid[326][4] c_mygrid[327][4] c_mygrid[328][4] c_mygrid[329][4] & - c_mygrid[330][4] c_mygrid[331][4] c_mygrid[332][4] c_mygrid[333][4] c_mygrid[334][4] c_mygrid[335][4] c_mygrid[336][4] c_mygrid[337][4] c_mygrid[338][4] c_mygrid[339][4] & - c_mygrid[340][4] c_mygrid[341][4] c_mygrid[342][4] c_mygrid[343][4] c_mygrid[344][4] c_mygrid[345][4] c_mygrid[346][4] c_mygrid[347][4] c_mygrid[348][4] c_mygrid[349][4] & - c_mygrid[350][4] c_mygrid[351][4] c_mygrid[352][4] c_mygrid[353][4] c_mygrid[354][4] c_mygrid[355][4] c_mygrid[356][4] c_mygrid[357][4] c_mygrid[358][4] c_mygrid[359][4] & - c_mygrid[360][4] c_mygrid[361][4] c_mygrid[362][4] c_mygrid[363][4] c_mygrid[364][4] c_mygrid[365][4] c_mygrid[366][4] c_mygrid[367][4] c_mygrid[368][4] c_mygrid[369][4] & - c_mygrid[370][4] c_mygrid[371][4] c_mygrid[372][4] c_mygrid[373][4] c_mygrid[374][4] c_mygrid[375][4] c_mygrid[376][4] c_mygrid[377][4] c_mygrid[378][4] c_mygrid[379][4] & - c_mygrid[380][4] c_mygrid[381][4] c_mygrid[382][4] c_mygrid[383][4] c_mygrid[384][4] + c_mygrid[1][7] c_mygrid[2][7] c_mygrid[3][7] c_mygrid[4][7] c_mygrid[5][7] c_mygrid[6][7] c_mygrid[7][7] c_mygrid[8][7] c_mygrid[9][7] & + c_mygrid[10][7] c_mygrid[11][7] c_mygrid[12][7] c_mygrid[13][7] c_mygrid[14][7] c_mygrid[15][7] c_mygrid[16][7] c_mygrid[17][7] c_mygrid[18][7] c_mygrid[19][7] & + c_mygrid[20][7] c_mygrid[21][7] c_mygrid[22][7] c_mygrid[23][7] c_mygrid[24][7] c_mygrid[25][7] c_mygrid[26][7] c_mygrid[27][7] c_mygrid[28][7] c_mygrid[29][7] & + c_mygrid[30][7] c_mygrid[31][7] c_mygrid[32][7] c_mygrid[33][7] c_mygrid[34][7] c_mygrid[35][7] c_mygrid[36][7] c_mygrid[37][7] c_mygrid[38][7] c_mygrid[39][7] & + c_mygrid[40][7] c_mygrid[41][7] c_mygrid[42][7] c_mygrid[43][7] c_mygrid[44][7] c_mygrid[45][7] c_mygrid[46][7] c_mygrid[47][7] c_mygrid[48][7] c_mygrid[49][7] & + c_mygrid[50][7] c_mygrid[51][7] c_mygrid[52][7] c_mygrid[53][7] c_mygrid[54][7] c_mygrid[55][7] c_mygrid[56][7] c_mygrid[57][7] c_mygrid[58][7] c_mygrid[59][7] & + c_mygrid[60][7] c_mygrid[61][7] c_mygrid[62][7] c_mygrid[63][7] c_mygrid[64][7] c_mygrid[65][7] c_mygrid[66][7] c_mygrid[67][7] c_mygrid[68][7] c_mygrid[69][7] & + c_mygrid[70][7] c_mygrid[71][7] c_mygrid[72][7] c_mygrid[73][7] c_mygrid[74][7] c_mygrid[75][7] c_mygrid[76][7] c_mygrid[77][7] c_mygrid[78][7] c_mygrid[79][7] & + c_mygrid[80][7] c_mygrid[81][7] c_mygrid[82][7] c_mygrid[83][7] c_mygrid[84][7] c_mygrid[85][7] c_mygrid[86][7] c_mygrid[87][7] c_mygrid[88][7] c_mygrid[89][7] & + c_mygrid[90][7] c_mygrid[91][7] c_mygrid[92][7] c_mygrid[93][7] c_mygrid[94][7] c_mygrid[95][7] c_mygrid[96][7] c_mygrid[97][7] c_mygrid[98][7] c_mygrid[99][7] & + c_mygrid[101][7] c_mygrid[102][7] c_mygrid[103][7] c_mygrid[104][7] c_mygrid[105][7] c_mygrid[106][7] c_mygrid[107][7] c_mygrid[108][7] c_mygrid[109][7] & + c_mygrid[110][7] c_mygrid[111][7] c_mygrid[112][7] c_mygrid[113][7] c_mygrid[114][7] c_mygrid[115][7] c_mygrid[116][7] c_mygrid[117][7] c_mygrid[118][7] c_mygrid[119][7] & + c_mygrid[120][7] c_mygrid[121][7] c_mygrid[122][7] c_mygrid[123][7] c_mygrid[124][7] c_mygrid[125][7] c_mygrid[126][7] c_mygrid[127][7] c_mygrid[128][7] c_mygrid[129][7] & + c_mygrid[130][7] c_mygrid[131][7] c_mygrid[132][7] c_mygrid[133][7] c_mygrid[134][7] c_mygrid[135][7] c_mygrid[136][7] c_mygrid[137][7] c_mygrid[138][7] c_mygrid[139][7] & + c_mygrid[140][7] c_mygrid[141][7] c_mygrid[142][7] c_mygrid[143][7] c_mygrid[144][7] c_mygrid[145][7] c_mygrid[146][7] c_mygrid[147][7] c_mygrid[148][7] c_mygrid[149][7] & + c_mygrid[150][7] c_mygrid[151][7] c_mygrid[152][7] c_mygrid[153][7] c_mygrid[154][7] c_mygrid[155][7] c_mygrid[156][7] c_mygrid[157][7] c_mygrid[158][7] c_mygrid[159][7] & + c_mygrid[160][7] c_mygrid[161][7] c_mygrid[162][7] c_mygrid[163][7] c_mygrid[164][7] c_mygrid[165][7] c_mygrid[166][7] c_mygrid[167][7] c_mygrid[168][7] c_mygrid[169][7] & + c_mygrid[170][7] c_mygrid[171][7] c_mygrid[172][7] c_mygrid[173][7] c_mygrid[174][7] c_mygrid[175][7] c_mygrid[176][7] c_mygrid[177][7] c_mygrid[178][7] c_mygrid[179][7] & + c_mygrid[180][7] c_mygrid[181][7] c_mygrid[182][7] c_mygrid[183][7] c_mygrid[184][7] c_mygrid[185][7] c_mygrid[186][7] c_mygrid[187][7] c_mygrid[188][7] c_mygrid[189][7] & + c_mygrid[190][7] c_mygrid[191][7] c_mygrid[192][7] c_mygrid[193][7] c_mygrid[194][7] c_mygrid[195][7] c_mygrid[196][7] c_mygrid[197][7] c_mygrid[198][7] c_mygrid[199][7] & + c_mygrid[201][7] c_mygrid[202][7] c_mygrid[203][7] c_mygrid[204][7] c_mygrid[205][7] c_mygrid[206][7] c_mygrid[207][7] c_mygrid[208][7] c_mygrid[209][7] & + c_mygrid[210][7] c_mygrid[211][7] c_mygrid[212][7] c_mygrid[213][7] c_mygrid[214][7] c_mygrid[215][7] c_mygrid[216][7] c_mygrid[217][7] c_mygrid[218][7] c_mygrid[219][7] & + c_mygrid[220][7] c_mygrid[221][7] c_mygrid[222][7] c_mygrid[223][7] c_mygrid[224][7] c_mygrid[225][7] c_mygrid[226][7] c_mygrid[227][7] c_mygrid[228][7] c_mygrid[229][7] & + c_mygrid[230][7] c_mygrid[231][7] c_mygrid[232][7] c_mygrid[233][7] c_mygrid[234][7] c_mygrid[235][7] c_mygrid[236][7] c_mygrid[237][7] c_mygrid[238][7] c_mygrid[239][7] & + c_mygrid[240][7] c_mygrid[241][7] c_mygrid[242][7] c_mygrid[243][7] c_mygrid[244][7] c_mygrid[245][7] c_mygrid[246][7] c_mygrid[247][7] c_mygrid[248][7] c_mygrid[249][7] & + c_mygrid[250][7] c_mygrid[251][7] c_mygrid[252][7] c_mygrid[253][7] c_mygrid[254][7] c_mygrid[255][7] c_mygrid[256][7] c_mygrid[257][7] c_mygrid[258][7] c_mygrid[259][7] & + c_mygrid[260][7] c_mygrid[261][7] c_mygrid[262][7] c_mygrid[263][7] c_mygrid[264][7] c_mygrid[265][7] c_mygrid[266][7] c_mygrid[267][7] c_mygrid[268][7] c_mygrid[269][7] & + c_mygrid[270][7] c_mygrid[271][7] c_mygrid[272][7] c_mygrid[273][7] c_mygrid[274][7] c_mygrid[275][7] c_mygrid[276][7] c_mygrid[277][7] c_mygrid[278][7] c_mygrid[279][7] & + c_mygrid[280][7] c_mygrid[281][7] c_mygrid[282][7] c_mygrid[283][7] c_mygrid[284][7] c_mygrid[285][7] c_mygrid[286][7] c_mygrid[287][7] c_mygrid[288][7] c_mygrid[289][7] & + c_mygrid[290][7] c_mygrid[291][7] c_mygrid[292][7] c_mygrid[293][7] c_mygrid[294][7] c_mygrid[295][7] c_mygrid[296][7] c_mygrid[297][7] c_mygrid[298][7] c_mygrid[299][7] & + c_mygrid[301][7] c_mygrid[302][7] c_mygrid[303][7] c_mygrid[304][7] c_mygrid[305][7] c_mygrid[306][7] c_mygrid[307][7] c_mygrid[308][7] c_mygrid[309][7] & + c_mygrid[310][7] c_mygrid[311][7] c_mygrid[312][7] c_mygrid[313][7] c_mygrid[314][7] c_mygrid[315][7] c_mygrid[316][7] c_mygrid[317][7] c_mygrid[318][7] c_mygrid[319][7] & + c_mygrid[320][7] c_mygrid[321][7] c_mygrid[322][7] c_mygrid[323][7] c_mygrid[324][7] c_mygrid[325][7] c_mygrid[326][7] c_mygrid[327][7] c_mygrid[328][7] c_mygrid[329][7] & + c_mygrid[330][7] c_mygrid[331][7] c_mygrid[332][7] c_mygrid[333][7] c_mygrid[334][7] c_mygrid[335][7] c_mygrid[336][7] c_mygrid[337][7] c_mygrid[338][7] c_mygrid[339][7] & + c_mygrid[340][7] c_mygrid[341][7] c_mygrid[342][7] c_mygrid[343][7] c_mygrid[344][7] c_mygrid[345][7] c_mygrid[346][7] c_mygrid[347][7] c_mygrid[348][7] c_mygrid[349][7] & + c_mygrid[350][7] c_mygrid[351][7] c_mygrid[352][7] c_mygrid[353][7] c_mygrid[354][7] c_mygrid[355][7] c_mygrid[356][7] c_mygrid[357][7] c_mygrid[358][7] c_mygrid[359][7] & + c_mygrid[360][7] c_mygrid[361][7] c_mygrid[362][7] c_mygrid[363][7] c_mygrid[364][7] c_mygrid[365][7] c_mygrid[366][7] c_mygrid[367][7] c_mygrid[368][7] c_mygrid[369][7] & + c_mygrid[370][7] c_mygrid[371][7] c_mygrid[372][7] c_mygrid[373][7] c_mygrid[374][7] c_mygrid[375][7] c_mygrid[376][7] c_mygrid[377][7] c_mygrid[378][7] c_mygrid[379][7] & + c_mygrid[380][7] c_mygrid[381][7] c_mygrid[382][7] c_mygrid[383][7] c_mygrid[384][7] thermo_modify norm yes From 4bfb5051237421cdef4501f9a3edb23f52589629 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Mon, 25 May 2020 16:35:50 -0600 Subject: [PATCH 10/19] Added reviewer response doc --- src/ReviewerResponse.tex | 59 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100644 src/ReviewerResponse.tex diff --git a/src/ReviewerResponse.tex b/src/ReviewerResponse.tex new file mode 100644 index 0000000000..e29d756d78 --- /dev/null +++ b/src/ReviewerResponse.tex @@ -0,0 +1,59 @@ +\pdfoutput=1 +\documentclass[onecolumn,prb,showpacs,aps,superscriptaddress]{revtex4-1} +%\documentclass[prb,showpacs,aps,superscriptaddress]{revtex4-1} + + +\usepackage[english]{babel} +\usepackage[ansinew]{inputenc} +\usepackage{times} +\usepackage{graphicx} +\usepackage{graphics} +\usepackage{amsmath} +\usepackage{amsfonts} +\usepackage{amssymb} +\usepackage{amsbsy} +\usepackage{dsfont} +\usepackage{epstopdf} +\usepackage{makeidx} +%\usepackage{subfigure} +\usepackage{color} +\usepackage{pgf} +\usepackage{bm} +\usepackage{tikz} +%\usepackage[normalem]{ulem} +\usepackage{graphicx} % needed for figures +\usepackage{dcolumn} % needed for some tables +\usepackage{bm} % for math +\usepackage{amssymb} % for math + + +\begin{document} + +{\bf Title Here} + + +\vspace{1cm} + +We thank the reviewers for their careful reading of the manuscript and their suggestions for improvement. We have updated the manuscript accordingly. Please see below for our detailed responses to the review comments. + +\vspace{2cm} + + +\textbf{\underline{Reviewer 1:}} +\begin{itemize} + \item \textbf{Review Comment 1:} Blah blah + \item \textbf{Author Reply 1:} Thank you for your review. We agree making the blha blah is a good idea. We added the following sentence to the Conclusions to address this point \\ + ``Blah blah'' + +\end{itemize} + +\textbf{\underline{Reviewer 2:}} +\begin{itemize} + \item \textbf{Review Comment 1:} Blah blah + \item \textbf{Author Reply 1:} Blah blah \\ + +\end{itemize} + +\end{document} + + From 4295dd2dbcc9aae860607e823c440169c1cc23ae Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Mon, 25 May 2020 17:09:32 -0600 Subject: [PATCH 11/19] Wrong number --- src/ReviewerResponse.tex | 59 ---------------------------------------- 1 file changed, 59 deletions(-) delete mode 100644 src/ReviewerResponse.tex diff --git a/src/ReviewerResponse.tex b/src/ReviewerResponse.tex deleted file mode 100644 index e29d756d78..0000000000 --- a/src/ReviewerResponse.tex +++ /dev/null @@ -1,59 +0,0 @@ -\pdfoutput=1 -\documentclass[onecolumn,prb,showpacs,aps,superscriptaddress]{revtex4-1} -%\documentclass[prb,showpacs,aps,superscriptaddress]{revtex4-1} - - -\usepackage[english]{babel} -\usepackage[ansinew]{inputenc} -\usepackage{times} -\usepackage{graphicx} -\usepackage{graphics} -\usepackage{amsmath} -\usepackage{amsfonts} -\usepackage{amssymb} -\usepackage{amsbsy} -\usepackage{dsfont} -\usepackage{epstopdf} -\usepackage{makeidx} -%\usepackage{subfigure} -\usepackage{color} -\usepackage{pgf} -\usepackage{bm} -\usepackage{tikz} -%\usepackage[normalem]{ulem} -\usepackage{graphicx} % needed for figures -\usepackage{dcolumn} % needed for some tables -\usepackage{bm} % for math -\usepackage{amssymb} % for math - - -\begin{document} - -{\bf Title Here} - - -\vspace{1cm} - -We thank the reviewers for their careful reading of the manuscript and their suggestions for improvement. We have updated the manuscript accordingly. Please see below for our detailed responses to the review comments. - -\vspace{2cm} - - -\textbf{\underline{Reviewer 1:}} -\begin{itemize} - \item \textbf{Review Comment 1:} Blah blah - \item \textbf{Author Reply 1:} Thank you for your review. We agree making the blha blah is a good idea. We added the following sentence to the Conclusions to address this point \\ - ``Blah blah'' - -\end{itemize} - -\textbf{\underline{Reviewer 2:}} -\begin{itemize} - \item \textbf{Review Comment 1:} Blah blah - \item \textbf{Author Reply 1:} Blah blah \\ - -\end{itemize} - -\end{document} - - From cf6570d81283a356915bd4379ca807837988ccc9 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 2 Jul 2021 15:06:19 -0600 Subject: [PATCH 12/19] Tweaked comments in grid examples --- examples/snap/in.grid | 3 +++ examples/snap/in.grid.tri | 26 ++++++++++++++++++++------ 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/examples/snap/in.grid b/examples/snap/in.grid index 2702fccad1..44c42673e0 100644 --- a/examples/snap/in.grid +++ b/examples/snap/in.grid @@ -1,5 +1,8 @@ # Demonstrate bispectrum computes +# CORRECTNESS: thermo output for c_mygrid[*][1] and c_mygrid[*][8] should +# match the values in dump_b: 108.173 3.21778 0.712238 7.06634 1.04273 + # Initialize simulation variable nsteps index 0 diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri index c984799e09..fc71c3a43f 100644 --- a/examples/snap/in.grid.tri +++ b/examples/snap/in.grid.tri @@ -1,5 +1,15 @@ # Demonstrate bispectrum computes +# This triclinic cell has 6 times the volume of a single bcc cube +# and contains 12 atoms. It is a 3x2x1 supercell +# with each unit cell containing 2 atoms and the +# reduced lattice vectors are [1 0 0], [1 1 0], and [1 1 1]. +# The grid is listed in x-fastest order +# CORRECTNESS: thermo output for c_mygrid[*][7] should contain +# the same 2 values that appear in c_mygrid[*][7] for in.grid, +# 7.0663376 and 13.808803, corresponding to atom and insterstitial sites, +# respectively, and with occurrences 12 and 36, respectively. + # Initialize simulation variable nsteps index 0 @@ -11,13 +21,17 @@ variable ngridx equal 3*${ngrid} variable ngridy equal 2*${ngrid} variable ngridz equal 1*${ngrid} +variable nrepx equal 1*${nrep} +variable nrepy equal 1*${nrep} +variable nrepz equal 1*${nrep} + units metal # generate the box and atom positions using a BCC lattice -variable nx equal ${nrep} -variable ny equal ${nrep} -variable nz equal ${nrep} +variable nx equal ${nrepx} +variable ny equal ${nrepy} +variable nz equal ${nrepz} boundary p p p @@ -78,7 +92,7 @@ thermo_style custom step temp ke pe vol & c_mygrid[70][7] c_mygrid[71][7] c_mygrid[72][7] c_mygrid[73][7] c_mygrid[74][7] c_mygrid[75][7] c_mygrid[76][7] c_mygrid[77][7] c_mygrid[78][7] c_mygrid[79][7] & c_mygrid[80][7] c_mygrid[81][7] c_mygrid[82][7] c_mygrid[83][7] c_mygrid[84][7] c_mygrid[85][7] c_mygrid[86][7] c_mygrid[87][7] c_mygrid[88][7] c_mygrid[89][7] & c_mygrid[90][7] c_mygrid[91][7] c_mygrid[92][7] c_mygrid[93][7] c_mygrid[94][7] c_mygrid[95][7] c_mygrid[96][7] c_mygrid[97][7] c_mygrid[98][7] c_mygrid[99][7] & - c_mygrid[101][7] c_mygrid[102][7] c_mygrid[103][7] c_mygrid[104][7] c_mygrid[105][7] c_mygrid[106][7] c_mygrid[107][7] c_mygrid[108][7] c_mygrid[109][7] & + c_mygrid[100][7] c_mygrid[101][7] c_mygrid[102][7] c_mygrid[103][7] c_mygrid[104][7] c_mygrid[105][7] c_mygrid[106][7] c_mygrid[107][7] c_mygrid[108][7] c_mygrid[109][7] & c_mygrid[110][7] c_mygrid[111][7] c_mygrid[112][7] c_mygrid[113][7] c_mygrid[114][7] c_mygrid[115][7] c_mygrid[116][7] c_mygrid[117][7] c_mygrid[118][7] c_mygrid[119][7] & c_mygrid[120][7] c_mygrid[121][7] c_mygrid[122][7] c_mygrid[123][7] c_mygrid[124][7] c_mygrid[125][7] c_mygrid[126][7] c_mygrid[127][7] c_mygrid[128][7] c_mygrid[129][7] & c_mygrid[130][7] c_mygrid[131][7] c_mygrid[132][7] c_mygrid[133][7] c_mygrid[134][7] c_mygrid[135][7] c_mygrid[136][7] c_mygrid[137][7] c_mygrid[138][7] c_mygrid[139][7] & @@ -88,7 +102,7 @@ thermo_style custom step temp ke pe vol & c_mygrid[170][7] c_mygrid[171][7] c_mygrid[172][7] c_mygrid[173][7] c_mygrid[174][7] c_mygrid[175][7] c_mygrid[176][7] c_mygrid[177][7] c_mygrid[178][7] c_mygrid[179][7] & c_mygrid[180][7] c_mygrid[181][7] c_mygrid[182][7] c_mygrid[183][7] c_mygrid[184][7] c_mygrid[185][7] c_mygrid[186][7] c_mygrid[187][7] c_mygrid[188][7] c_mygrid[189][7] & c_mygrid[190][7] c_mygrid[191][7] c_mygrid[192][7] c_mygrid[193][7] c_mygrid[194][7] c_mygrid[195][7] c_mygrid[196][7] c_mygrid[197][7] c_mygrid[198][7] c_mygrid[199][7] & - c_mygrid[201][7] c_mygrid[202][7] c_mygrid[203][7] c_mygrid[204][7] c_mygrid[205][7] c_mygrid[206][7] c_mygrid[207][7] c_mygrid[208][7] c_mygrid[209][7] & + c_mygrid[200][7] c_mygrid[201][7] c_mygrid[202][7] c_mygrid[203][7] c_mygrid[204][7] c_mygrid[205][7] c_mygrid[206][7] c_mygrid[207][7] c_mygrid[208][7] c_mygrid[209][7] & c_mygrid[210][7] c_mygrid[211][7] c_mygrid[212][7] c_mygrid[213][7] c_mygrid[214][7] c_mygrid[215][7] c_mygrid[216][7] c_mygrid[217][7] c_mygrid[218][7] c_mygrid[219][7] & c_mygrid[220][7] c_mygrid[221][7] c_mygrid[222][7] c_mygrid[223][7] c_mygrid[224][7] c_mygrid[225][7] c_mygrid[226][7] c_mygrid[227][7] c_mygrid[228][7] c_mygrid[229][7] & c_mygrid[230][7] c_mygrid[231][7] c_mygrid[232][7] c_mygrid[233][7] c_mygrid[234][7] c_mygrid[235][7] c_mygrid[236][7] c_mygrid[237][7] c_mygrid[238][7] c_mygrid[239][7] & @@ -98,7 +112,7 @@ thermo_style custom step temp ke pe vol & c_mygrid[270][7] c_mygrid[271][7] c_mygrid[272][7] c_mygrid[273][7] c_mygrid[274][7] c_mygrid[275][7] c_mygrid[276][7] c_mygrid[277][7] c_mygrid[278][7] c_mygrid[279][7] & c_mygrid[280][7] c_mygrid[281][7] c_mygrid[282][7] c_mygrid[283][7] c_mygrid[284][7] c_mygrid[285][7] c_mygrid[286][7] c_mygrid[287][7] c_mygrid[288][7] c_mygrid[289][7] & c_mygrid[290][7] c_mygrid[291][7] c_mygrid[292][7] c_mygrid[293][7] c_mygrid[294][7] c_mygrid[295][7] c_mygrid[296][7] c_mygrid[297][7] c_mygrid[298][7] c_mygrid[299][7] & - c_mygrid[301][7] c_mygrid[302][7] c_mygrid[303][7] c_mygrid[304][7] c_mygrid[305][7] c_mygrid[306][7] c_mygrid[307][7] c_mygrid[308][7] c_mygrid[309][7] & + c_mygrid[300][7] c_mygrid[301][7] c_mygrid[302][7] c_mygrid[303][7] c_mygrid[304][7] c_mygrid[305][7] c_mygrid[306][7] c_mygrid[307][7] c_mygrid[308][7] c_mygrid[309][7] & c_mygrid[310][7] c_mygrid[311][7] c_mygrid[312][7] c_mygrid[313][7] c_mygrid[314][7] c_mygrid[315][7] c_mygrid[316][7] c_mygrid[317][7] c_mygrid[318][7] c_mygrid[319][7] & c_mygrid[320][7] c_mygrid[321][7] c_mygrid[322][7] c_mygrid[323][7] c_mygrid[324][7] c_mygrid[325][7] c_mygrid[326][7] c_mygrid[327][7] c_mygrid[328][7] c_mygrid[329][7] & c_mygrid[330][7] c_mygrid[331][7] c_mygrid[332][7] c_mygrid[333][7] c_mygrid[334][7] c_mygrid[335][7] c_mygrid[336][7] c_mygrid[337][7] c_mygrid[338][7] c_mygrid[339][7] & From 442585313c4e8c9ba03419987a44266412acaa33 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 2 Jul 2021 15:11:14 -0600 Subject: [PATCH 13/19] Test --- examples/snap/in.grid.tri | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri index fc71c3a43f..da87774107 100644 --- a/examples/snap/in.grid.tri +++ b/examples/snap/in.grid.tri @@ -1,4 +1,4 @@ -# Demonstrate bispectrum computes +# Demonstrate bispectrum computes for triclinic cell # This triclinic cell has 6 times the volume of a single bcc cube # and contains 12 atoms. It is a 3x2x1 supercell From e102864c2dd9bf56b9a87808e98ac9b7042cff41 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 2 Jul 2021 15:13:40 -0600 Subject: [PATCH 14/19] Test2 --- examples/snap/in.grid.tri | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri index da87774107..a5b2453c3b 100644 --- a/examples/snap/in.grid.tri +++ b/examples/snap/in.grid.tri @@ -1,6 +1,7 @@ # Demonstrate bispectrum computes for triclinic cell -# This triclinic cell has 6 times the volume of a single bcc cube +# This triclinic cell has 6 times the volume of the single +# unit cell used by in.grid # and contains 12 atoms. It is a 3x2x1 supercell # with each unit cell containing 2 atoms and the # reduced lattice vectors are [1 0 0], [1 1 0], and [1 1 1]. From 01475cb3a84005ddde51411c97c4f5f1eb353d71 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 2 Jul 2021 15:35:59 -0600 Subject: [PATCH 15/19] Test3 --- examples/snap/in.grid.tri | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri index a5b2453c3b..201f763eeb 100644 --- a/examples/snap/in.grid.tri +++ b/examples/snap/in.grid.tri @@ -131,3 +131,4 @@ dump mydump_b all custom 1000 dump_b id c_b[*] run 0 + From 39039d261f39966bad202ded3c5244065f323aa3 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 2 Jul 2021 15:37:49 -0600 Subject: [PATCH 16/19] Test6 --- examples/snap/in.grid.tri | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri index 201f763eeb..f375c18c42 100644 --- a/examples/snap/in.grid.tri +++ b/examples/snap/in.grid.tri @@ -132,3 +132,4 @@ dump mydump_b all custom 1000 dump_b id c_b[*] run 0 + From e17ace385ddc6d00a485c1c8fd04989526462d6f Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 2 Jul 2021 17:47:45 -0600 Subject: [PATCH 17/19] First pass at distributed memory for grid --- src/SNAP/compute_sna_grid.cpp | 93 +++++++------- src/compute_grid.cpp | 224 ++++++++++++++++++++-------------- src/compute_grid.h | 17 ++- 3 files changed, 194 insertions(+), 140 deletions(-) diff --git a/src/SNAP/compute_sna_grid.cpp b/src/SNAP/compute_sna_grid.cpp index a9f3b6c3fb..196819008c 100644 --- a/src/SNAP/compute_sna_grid.cpp +++ b/src/SNAP/compute_sna_grid.cpp @@ -182,66 +182,75 @@ void ComputeSNAGrid::compute_array() snaptr->grow_rij(ntotal); - printf("ngrid = %d\n",ngrid); - for (int igrid = 0; igrid < ngrid; igrid++) { - if (!grid_local[igrid]) continue; - const double xtmp = grid[igrid][0]; - const double ytmp = grid[igrid][1]; - const double ztmp = grid[igrid][2]; + 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 double xtmp = grid[igrid][0]; + const double ytmp = grid[igrid][1]; + const double ztmp = grid[igrid][2]; - // rij[][3] = displacements between atom I and those neighbors - // inside = indices of neighbors of I within cutoff - // typej = types of neighbors of I within cutoff + // rij[][3] = displacements between atom I and those neighbors + // inside = indices of neighbors of I within cutoff + // typej = types of neighbors of I within cutoff - int ninside = 0; - for (int j = 0; j < ntotal; j++) { + int ninside = 0; + for (int j = 0; j < ntotal; j++) { - // check that j is in compute group + // check that j is in compute group - if (!(mask[j] & groupbit)) continue; + if (!(mask[j] & groupbit)) continue; - 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; - int jtype = type[j]; - if (rsq < cutsq[jtype][jtype] && rsq>1e-20) { - // printf("ninside = %d\n",ninside); - 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; - ninside++; - } - } + 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; + int jtype = type[j]; + 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; + ninside++; + } + } - snaptr->compute_ui(ninside); - snaptr->compute_zi(); - snaptr->compute_bi(); - for (int icoeff = 0; icoeff < ncoeff; icoeff++) - grid[igrid][size_array_cols_base+icoeff] = snaptr->blist[icoeff]; - // printf("igrid = %d %g %g %g %d B0 = %g\n",igrid,xtmp,ytmp,ztmp,ninside,sna[igrid][size_array_cols_base+0]); - if (quadraticflag) { - int ncount = ncoeff; - for (int icoeff = 0; icoeff < ncoeff; icoeff++) { - double bi = snaptr->blist[icoeff]; + snaptr->compute_ui(ninside); + snaptr->compute_zi(); + snaptr->compute_bi(); + for (int icoeff = 0; icoeff < ncoeff; icoeff++) + gridlocal[size_array_cols_base+icoeff][iz][iy][ix] = snaptr->blist[icoeff]; + if (quadraticflag) { + int ncount = ncoeff; + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + double bi = snaptr->blist[icoeff]; - // diagonal element of quadratic matrix + // diagonal element of quadratic matrix - grid[igrid][size_array_cols_base+ncount++] = 0.5*bi*bi; + gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = 0.5*bi*bi; // upper-triangular elements of quadratic matrix for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) - grid[igrid][size_array_cols_base+ncount++] = bi*snaptr->blist[jcoeff]; + gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = bi*snaptr->blist[jcoeff]; } } } + + 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; + 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],ngrid*size_array_cols,MPI_DOUBLE,MPI_SUM,world); + } + /* ---------------------------------------------------------------------- memory usage ------------------------------------------------------------------------- */ diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp index 011dafa030..b16a51b739 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -21,13 +21,14 @@ #include "force.h" #include "memory.h" #include "error.h" +#include "comm.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), grid(NULL), grid_local(NULL) + Compute(lmp, narg, arg), grid(NULL), local_flags(NULL), gridlocal(NULL) { if (narg < 6) error->all(FLERR,"Illegal compute grid command"); @@ -59,7 +60,8 @@ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : ComputeGrid::~ComputeGrid() { memory->destroy(grid); - memory->destroy(grid_local); + memory->destroy(local_flags); + memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo); } /* ---------------------------------------------------------------------- */ @@ -72,7 +74,118 @@ void ComputeGrid::init() void ComputeGrid::setup() { - + + set_grid_global(); + set_grid_local(); + allocate(); + assign_coords(); + assign_local_flags(); +} + +/* ---------------------------------------------------------------------- + convert global array index to box coords +------------------------------------------------------------------------- */ + +void ComputeGrid::grid2x(int igrid, double *x) +{ + 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; + + if (triclinic) domain->lamda2x(x, x); + +} + +/* ---------------------------------------------------------------------- + check if grid point is local +------------------------------------------------------------------------- */ + +int ComputeGrid::check_local(int igrid) +{ + double x[3]; + + 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; + + int islocal = + x[0] >= sublo[0] && x[0] < subhi[0] && + x[1] >= sublo[1] && x[1] < subhi[1] && + x[2] >= sublo[2] && x[2] < subhi[2]; + + return islocal; +} + +/* ---------------------------------------------------------------------- + copy coords to global array +------------------------------------------------------------------------- */ + +void ComputeGrid::assign_coords() +{ + double x[3]; + for (int igrid = 0; igrid < ngrid; igrid++) { + grid2x(igrid,x); + grid[igrid][0] = x[0]; + grid[igrid][1] = x[1]; + grid[igrid][2] = x[2]; + } +} + +/* ---------------------------------------------------------------------- + copy coords to global array +------------------------------------------------------------------------- */ + +void ComputeGrid::assign_local_flags() +{ + double x[3]; + for (int igrid = 0; igrid < ngrid; igrid++) { + if (check_local(igrid)) + local_flags[igrid] = 1; + else { + local_flags[igrid] = 0; + memset(grid[igrid],0,size_array_cols); + } + } +} + +/* ---------------------------------------------------------------------- + free and reallocate arrays +------------------------------------------------------------------------- */ + +void ComputeGrid::allocate() +{ + // allocate arrays + + memory->destroy(grid); + memory->destroy(local_flags); + memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo); + memory->create(grid,size_array_rows,size_array_cols,"grid:grid"); + memory->create(gridall,size_array_rows,size_array_cols,"grid:gridall"); + memory->create(local_flags,size_array_rows,"grid:local_flags"); + memory->create4d_offset(gridlocal,size_array_cols,nzlo,nzhi,nylo,nyhi, + nxlo,nxhi,"grid:gridlocal"); + array = gridall; +} + + +/* ---------------------------------------------------------------------- + set global grid +------------------------------------------------------------------------- */ + +void ComputeGrid::set_grid_global() +{ // calculate grid layout triclinic = domain->triclinic; @@ -100,105 +213,31 @@ void ComputeGrid::setup() delx = 1.0/delxinv; dely = 1.0/delyinv; delz = 1.0/delzinv; - - allocate(); - assign_grid_coords(); - assign_grid_local(); } /* ---------------------------------------------------------------------- - convert global array index to box coords + set local subset of grid that I own + n xyz lo/hi = 3d brick that I own (inclusive) ------------------------------------------------------------------------- */ -void ComputeGrid::grid2x(int igrid, double *x) +void ComputeGrid::set_grid_local() { - int iz = igrid / (nx*ny); - igrid -= iz * (nx*ny); - int iy = igrid / nx; - igrid -= iy * nx; - int ix = igrid; + // global indices of grid range from 0 to N-1 + // nlo,nhi = lower/upper limits of the 3d sub-brick of + // global grid that I own without ghost cells - x[0] = ix*delx; - x[1] = iy*dely; - x[2] = iz*delz; + nxlo = static_cast (comm->xsplit[comm->myloc[0]] * nx); + nxhi = static_cast (comm->xsplit[comm->myloc[0]+1] * nx) - 1; - if (triclinic) domain->lamda2x(x, x); + nylo = static_cast (comm->ysplit[comm->myloc[1]] * ny); + nyhi = static_cast (comm->ysplit[comm->myloc[1]+1] * ny) - 1; + nzlo = static_cast (comm->zsplit[comm->myloc[2]] * nz); + nzhi = static_cast (comm->zsplit[comm->myloc[2]+1] * nz) - 1; + + ngridlocal = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1); } -/* ---------------------------------------------------------------------- - check if grid point is local -------------------------------------------------------------------------- */ - -int ComputeGrid::check_grid_local(int igrid) -{ - double x[3]; - - 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; - - int islocal = - x[0] >= sublo[0] && x[0] < subhi[0] && - x[1] >= sublo[1] && x[1] < subhi[1] && - x[2] >= sublo[2] && x[2] < subhi[2]; - - return islocal; -} - -/* ---------------------------------------------------------------------- - copy coords to global array -------------------------------------------------------------------------- */ - -void ComputeGrid::assign_grid_coords() -{ - double x[3]; - for (int igrid = 0; igrid < ngrid; igrid++) { - grid2x(igrid,x); - grid[igrid][0] = x[0]; - grid[igrid][1] = x[1]; - grid[igrid][2] = x[2]; - } -} - -/* ---------------------------------------------------------------------- - copy coords to global array -------------------------------------------------------------------------- */ - -void ComputeGrid::assign_grid_local() -{ - double x[3]; - for (int igrid = 0; igrid < ngrid; igrid++) { - if (check_grid_local(igrid)) - grid_local[igrid] = 1; - else { - grid_local[igrid] = 0; - memset(grid[igrid],0,size_array_cols); - } - } -} - -/* ---------------------------------------------------------------------- - free and reallocate arrays -------------------------------------------------------------------------- */ - -void ComputeGrid::allocate() -{ - // grow global array if necessary - - memory->destroy(grid); - memory->destroy(grid_local); - memory->create(grid,size_array_rows,size_array_cols,"grid:grid"); - memory->create(gridall,size_array_rows,size_array_cols,"grid:gridall"); - memory->create(grid_local,size_array_rows,"grid:grid_local"); - array = gridall; -} /* ---------------------------------------------------------------------- memory usage of local data ------------------------------------------------------------------------- */ @@ -207,6 +246,7 @@ double ComputeGrid::memory_usage() { double nbytes = size_array_rows*size_array_cols * sizeof(double); // grid - nbytes += size_array_rows*sizeof(int); // grid_local + nbytes += size_array_rows*sizeof(int); // local_flags + nbytes += size_array_cols*ngridlocal*sizeof(double); // gridlocal return nbytes; } diff --git a/src/compute_grid.h b/src/compute_grid.h index a1c938db2e..7e1a4d706b 100644 --- a/src/compute_grid.h +++ b/src/compute_grid.h @@ -30,11 +30,14 @@ class ComputeGrid : public Compute { double memory_usage(); protected: - int nx, ny, nz; // grid dimensions - int ngrid; // number of grid points + int nx, ny, nz; // global grid dimensions + int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive + int ngrid; // number of global grid points + 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) @@ -43,12 +46,14 @@ class ComputeGrid : public Compute { 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 *grid_local; // local flag for each grid point + int *local_flags; // local flag for each grid point void allocate(); void grid2x(int, double*); // convert grid point to coord - void assign_grid_coords(); // assign coords for grid - void assign_grid_local(); // set local flag for each grid point - int check_grid_local(int); // check if grid point is local + void assign_coords(); // assign coords for grid + void assign_local_flags(); // set local flag for each grid point + int check_local(int); // check if grid point is local + void set_grid_global(); // set global grid + void set_grid_local(); // set bounds for local grid private: }; From 4c22f094de1c0a5632c3514974f75bc6e2c85c78 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 2 Jul 2021 18:15:55 -0600 Subject: [PATCH 18/19] Minor tweak --- src/compute_grid.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp index b16a51b739..1939dd544e 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -74,7 +74,6 @@ void ComputeGrid::init() void ComputeGrid::setup() { - set_grid_global(); set_grid_local(); allocate(); @@ -149,7 +148,6 @@ void ComputeGrid::assign_coords() void ComputeGrid::assign_local_flags() { - double x[3]; for (int igrid = 0; igrid < ngrid; igrid++) { if (check_local(igrid)) local_flags[igrid] = 1; From 07db7a40955ebb356ea16166985b1b296d1d8bb3 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Wed, 14 Jul 2021 13:35:05 -0600 Subject: [PATCH 19/19] Changed to different check_local() --- src/SNAP/compute_sna_grid.cpp | 6 +-- src/compute_grid.cpp | 75 ++++++++++++++++++++++++++--------- src/compute_grid.h | 3 ++ 3 files changed, 63 insertions(+), 21 deletions(-) diff --git a/src/SNAP/compute_sna_grid.cpp b/src/SNAP/compute_sna_grid.cpp index 196819008c..d2c2ae74ca 100644 --- a/src/SNAP/compute_sna_grid.cpp +++ b/src/SNAP/compute_sna_grid.cpp @@ -231,10 +231,10 @@ void ComputeSNAGrid::compute_array() gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = 0.5*bi*bi; - // upper-triangular elements of quadratic matrix + // upper-triangular elements of quadratic matrix - for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) - gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = bi*snaptr->blist[jcoeff]; + for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) + gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = bi*snaptr->blist[jcoeff]; } } } diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp index 1939dd544e..3cbcb2c2db 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -53,6 +53,7 @@ ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : size_array_rows = ngrid = nx*ny*nz; size_array_cols_base = 3; + gridlocal_allocated = 0; } /* ---------------------------------------------------------------------- */ @@ -61,7 +62,8 @@ ComputeGrid::~ComputeGrid() { memory->destroy(grid); memory->destroy(local_flags); - memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo); + if (gridlocal_allocated) + memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo); } /* ---------------------------------------------------------------------- */ @@ -101,28 +103,59 @@ void ComputeGrid::grid2x(int igrid, double *x) } +/* ---------------------------------------------------------------------- + convert global array index to box coords +------------------------------------------------------------------------- */ + +void ComputeGrid::grid2ix(int igrid, int& ix, int& iy, int& iz) +{ + iz = igrid / (nx*ny); + igrid -= iz * (nx*ny); + iy = igrid / nx; + igrid -= iy * nx; + ix = igrid; +} + +// /* ---------------------------------------------------------------------- +// check if grid point is local +// ------------------------------------------------------------------------- */ + +// int ComputeGrid::check_local(int igrid) +// { +// double x[3]; + +// 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; + +// int islocal = +// x[0] >= sublo[0] && x[0] < subhi[0] && +// x[1] >= sublo[1] && x[1] < subhi[1] && +// x[2] >= sublo[2] && x[2] < subhi[2]; + +// return islocal; +// } + /* ---------------------------------------------------------------------- check if grid point is local ------------------------------------------------------------------------- */ int ComputeGrid::check_local(int igrid) { - double x[3]; + int ix, iy, iz; - 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; + grid2ix(igrid, ix, iy, iz); int islocal = - x[0] >= sublo[0] && x[0] < subhi[0] && - x[1] >= sublo[1] && x[1] < subhi[1] && - x[2] >= sublo[2] && x[2] < subhi[2]; + ix >= nxlo && ix <= nxhi && + iy >= nylo && iy <= nyhi && + iz >= nzlo && iz <= nzhi; return islocal; } @@ -153,7 +186,7 @@ void ComputeGrid::assign_local_flags() local_flags[igrid] = 1; else { local_flags[igrid] = 0; - memset(grid[igrid],0,size_array_cols); + memset(grid[igrid],0,size_array_cols * sizeof(double)); } } } @@ -168,12 +201,16 @@ void ComputeGrid::allocate() memory->destroy(grid); memory->destroy(local_flags); - memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo); + if (gridlocal_allocated) + memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo); memory->create(grid,size_array_rows,size_array_cols,"grid:grid"); memory->create(gridall,size_array_rows,size_array_cols,"grid:gridall"); memory->create(local_flags,size_array_rows,"grid:local_flags"); - memory->create4d_offset(gridlocal,size_array_cols,nzlo,nzhi,nylo,nyhi, - nxlo,nxhi,"grid:gridlocal"); + 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"); + } array = gridall; } @@ -234,6 +271,8 @@ void ComputeGrid::set_grid_local() nzhi = static_cast (comm->zsplit[comm->myloc[2]+1] * nz) - 1; ngridlocal = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1); + + printf("me = %d n = %d x %d %d y %d %d z %d %d \n", comm->me, ngridlocal, nxlo, nxhi, nylo, nyhi, nzlo, nzhi); } /* ---------------------------------------------------------------------- diff --git a/src/compute_grid.h b/src/compute_grid.h index 7e1a4d706b..1b2797732d 100644 --- a/src/compute_grid.h +++ b/src/compute_grid.h @@ -47,8 +47,11 @@ class ComputeGrid : public Compute { double cutmax; // largest cutoff distance int size_array_cols_base; // number of columns used for coords, etc. int *local_flags; // local flag for each grid point + int gridlocal_allocated; // shows if gridlocal allocated + void allocate(); void grid2x(int, double*); // convert grid point to coord + void grid2ix(int, int&, int&, int&); // convert grid point to ix, iy, iz void assign_coords(); // assign coords for grid void assign_local_flags(); // set local flag for each grid point int check_local(int); // check if grid point is local