From d972bf9010f0217fb10ead0c49f7a76607aa1cc0 Mon Sep 17 00:00:00 2001 From: Lenz Fiedler Date: Thu, 18 Apr 2024 21:07:21 +0200 Subject: [PATCH] Added/moved files originally authored by Aidan Thompson; small adjustment made to get it running Co-authored-by: Aidan Thompson --- src/ML-SNAP/pair_grid.cpp | 330 +++++++++++++++++++++++++ src/ML-SNAP/pair_grid.h | 87 +++++++ src/ML-SNAP/pair_sna_grid.cpp | 449 ++++++++++++++++++++++++++++++++++ src/ML-SNAP/pair_sna_grid.h | 79 ++++++ 4 files changed, 945 insertions(+) create mode 100644 src/ML-SNAP/pair_grid.cpp create mode 100644 src/ML-SNAP/pair_grid.h create mode 100644 src/ML-SNAP/pair_sna_grid.cpp create mode 100644 src/ML-SNAP/pair_sna_grid.h diff --git a/src/ML-SNAP/pair_grid.cpp b/src/ML-SNAP/pair_grid.cpp new file mode 100644 index 00000000000..0391461cefd --- /dev/null +++ b/src/ML-SNAP/pair_grid.cpp @@ -0,0 +1,330 @@ +/* ---------------------------------------------------------------------- + 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 "pair_grid.h" +#include +#include +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "domain.h" +#include "force.h" +#include "memory.h" +#include "error.h" +#include "comm.h" + +#define BETA_CONST 1.0e-2 +#define BETA_IGRID 1.0e-2 +#define BETA_ICOL 1.0e-2 + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairGrid::PairGrid(LAMMPS *lmp) : + Pair(lmp), gridlocal(nullptr), alocal(nullptr), beta(nullptr) +{ + single_enable = 0; + restartinfo = 0; + one_coeff = 1; + manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; + + ndesc = 0; + ngridlocal = 0; + + ndesc_base = 6; + gridlocal_allocated = 0; + beta_max = 0; + beta = nullptr; +} + +/* ---------------------------------------------------------------------- */ + +PairGrid::~PairGrid() +{ + if (copymode) return; + deallocate_grid(); +} + +/* ---------------------------------------------------------------------- */ + +void PairGrid::setup() +{ + deallocate_grid(); + set_grid_global(); + set_grid_local(); + allocate_grid(); + assign_coords(); +} + +/* ---------------------------------------------------------------------- + convert global array indexes to box coords +------------------------------------------------------------------------- */ + +void PairGrid::grid2x(int ix, int iy, int iz, double *x) +{ + x[0] = ix*delx; + x[1] = iy*dely; + x[2] = iz*delz; + + if (triclinic) domain->lamda2x(x, x); +} + +/* ---------------------------------------------------------------------- + create arrays +------------------------------------------------------------------------- */ + +void PairGrid::allocate_grid() +{ + if (nxlo <= nxhi && nylo <= nyhi && nzlo <= nzhi) { + gridlocal_allocated = 1; + memory->create4d_offset(gridlocal,ndesc,nzlo,nzhi,nylo,nyhi, + nxlo,nxhi,"pair/grid:gridlocal"); + memory->create(alocal, ngridlocal, ndesc, "pair/grid:alocal"); + memory->create(beta, ngridlocal, ndesc-ndesc_base, "pair/grid:beta"); + } +} + +/* ---------------------------------------------------------------------- + free arrays +------------------------------------------------------------------------- */ + +void PairGrid::deallocate_grid() +{ + if (gridlocal_allocated) { + gridlocal_allocated = 0; + memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo); + memory->destroy(alocal); + memory->destroy(beta); + } +} + +/* ---------------------------------------------------------------------- + set global grid +------------------------------------------------------------------------- */ + +void PairGrid::set_grid_global() +{ + // calculate grid layout + + triclinic = domain->triclinic; + + 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]; + 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; +} + +/* ---------------------------------------------------------------------- + set local subset of grid that I own + n xyz lo/hi = 3d brick that I own (inclusive) +------------------------------------------------------------------------- */ + +void PairGrid::set_grid_local() +{ + // nx,ny,nz = extent of global grid + // indices into the global grid range from 0 to N-1 in each dim + // if grid point is inside my sub-domain I own it, + // this includes sub-domain lo boundary but excludes hi boundary + // ixyz lo/hi = inclusive lo/hi bounds of global grid sub-brick I own + // if proc owns no grid cells in a dim, then ilo > ihi + // if 2 procs share a boundary a grid point is exactly on, + // the 2 equality if tests insure a consistent decision + // as to which proc owns it + + double xfraclo,xfrachi,yfraclo,yfrachi,zfraclo,zfrachi; + + if (comm->layout != Comm::LAYOUT_TILED) { + xfraclo = comm->xsplit[comm->myloc[0]]; + xfrachi = comm->xsplit[comm->myloc[0]+1]; + yfraclo = comm->ysplit[comm->myloc[1]]; + yfrachi = comm->ysplit[comm->myloc[1]+1]; + zfraclo = comm->zsplit[comm->myloc[2]]; + zfrachi = comm->zsplit[comm->myloc[2]+1]; + } else { + xfraclo = comm->mysplit[0][0]; + xfrachi = comm->mysplit[0][1]; + yfraclo = comm->mysplit[1][0]; + yfrachi = comm->mysplit[1][1]; + zfraclo = comm->mysplit[2][0]; + zfrachi = comm->mysplit[2][1]; + } + + nxlo = static_cast (xfraclo * nx); + if (1.0*nxlo != xfraclo*nx) nxlo++; + nxhi = static_cast (xfrachi * nx); + if (1.0*nxhi == xfrachi*nx) nxhi--; + + nylo = static_cast (yfraclo * ny); + if (1.0*nylo != yfraclo*ny) nylo++; + nyhi = static_cast (yfrachi * ny); + if (1.0*nyhi == yfrachi*ny) nyhi--; + + nzlo = static_cast (zfraclo * nz); + if (1.0*nzlo != zfraclo*nz) nzlo++; + nzhi = static_cast (zfrachi * nz); + if (1.0*nzhi == zfrachi*nz) nzhi--; + + ngridlocal = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1); +} + +/* ---------------------------------------------------------------------- + copy coords to local array +------------------------------------------------------------------------- */ + +void PairGrid::assign_coords() +{ + int igrid = 0; + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + alocal[igrid][0] = ix; + alocal[igrid][1] = iy; + alocal[igrid][2] = iz; + double xgrid[3]; + grid2x(ix, iy, iz, xgrid); + alocal[igrid][3] = xgrid[0]; + alocal[igrid][4] = xgrid[1]; + alocal[igrid][5] = xgrid[2]; + igrid++; + } +} + +/* ---------------------------------------------------------------------- + copy the 4d gridlocal array values to the 2d local array +------------------------------------------------------------------------- */ + +void PairGrid::copy_gridlocal_to_local_array() +{ + int igrid = 0; + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + for (int icol = ndesc_base; icol < ndesc; icol++) + alocal[igrid][icol] = gridlocal[icol][iz][iy][ix]; + igrid++; + } +} + +/* ---------------------------------------------------------------------- + calculate beta +------------------------------------------------------------------------- */ + +// this is a proxy for a call to the energy model +// beta is dE/dB^i, the derivative of the total +// energy w.r.t. to descriptors of grid point i + +void PairGrid::compute_beta() +{ + int igrid = 0; + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + for (int icol = 0; icol < ndesc-ndesc_base; icol++) + beta[igrid][icol] = BETA_CONST + BETA_IGRID*igrid + BETA_ICOL*icol; + igrid++; + } +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairGrid::allocate() +{ + allocated = 1; + int n = atom->ntypes; + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + map = new int[n+1]; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairGrid::settings(int narg, char ** arg) +{ + if (narg < 4) error->all(FLERR,"Illegal pair style command"); + int iarg0 = 0; + int iarg = iarg0; + if (strcmp(arg[iarg],"grid") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal pair grid command"); + nx = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + ny = utils::inumeric(FLERR,arg[iarg+2],false,lmp); + nz = utils::inumeric(FLERR,arg[iarg+3],false,lmp); + if (nx <= 0 || ny <= 0 || nz <= 0) + error->all(FLERR,"All grid/local dimensions must be positive"); + iarg += 4; + } else error->all(FLERR,"Illegal pair grid command"); + nargbase = iarg - iarg0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairGrid::coeff(int narg, char **arg) +{ + if (!allocated) allocate(); + if (narg < 2) error->all(FLERR,"Incorrect args for pair coefficients"); + // if (narg != 2 + atom->ntypes) error->all(FLERR,"Incorrect args for pair coefficients"); + + map_element2type(narg-2,arg+2); + // map_element2type(0,nullptr); + +} + +/* ---------------------------------------------------------------------- + return maximum force cut off distance +------------------------------------------------------------------------- */ + +double PairGrid::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + return cutmax; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double PairGrid::memory_usage() +{ + int nbytes = ndesc*ngridlocal*sizeof(double); // gridlocal + return nbytes; +} \ No newline at end of file diff --git a/src/ML-SNAP/pair_grid.h b/src/ML-SNAP/pair_grid.h new file mode 100644 index 00000000000..ccf0a7b4b83 --- /dev/null +++ b/src/ML-SNAP/pair_grid.h @@ -0,0 +1,87 @@ +/* -*- 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 PAIR_CLASS +// clang-format off +PairStyle(grid,PairGrid); +// clang-format on +#else + +#ifndef LMP_PAIR_GRID_H +#define LMP_PAIR_GRID_H + +#include "pair.h" + +namespace LAMMPS_NS { + + class PairGrid : public Pair { + public: + PairGrid(class LAMMPS *); + virtual ~PairGrid(); + virtual void init_style(){}; + void setup(); + virtual void compute(int, int) { + printf("DANGER! This function should always be overridden by child\n"); + }; + + void settings(int, char **); + virtual void coeff(int, char **); + virtual double init_one(int, int); + double memory_usage(); + + protected: + int nx, ny, nz; // global grid dimensions + int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive + int ngridlocal; // number of local grid points + int nvalues; // number of values per grid point + double ****gridlocal; // local grid, redundant w.r.t. alocal + double **alocal; // pointer to Compute::array_local + int triclinic; // triclinic flag + double *boxlo, *prd; // box info (units real/ortho or reduced/tri) + double *sublo, *subhi; // subdomain info (units real/ortho or reduced/tri) + double delxinv,delyinv,delzinv; // inverse grid spacing + double delx,dely,delz; // grid spacing + int nargbase; // number of base class args + double cutmax; // largest cutoff distance + int ndesc; // number of descriptors + int ndesc_base; // number of columns used for coords, etc. + int gridlocal_allocated; // shows if gridlocal allocated + double **beta; // betas for all local grid points in list + int beta_max; // length of beta + + void allocate(); // allocate pairstyle arrays + void allocate_grid(); // create grid arrays + void deallocate_grid(); // free grid arrays + void grid2x(int, int, int, double*); // convert global indices to coordinates + void set_grid_global(); // set global grid + void set_grid_local(); // set bounds for local grid + void assign_coords(); // assign coords for grid + void copy_gridlocal_to_local_array();// copy 4d gridlocal array to 2d local array + void compute_beta(); // get betas from someplace + private: + }; + +} + +#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. + +*/ \ No newline at end of file diff --git a/src/ML-SNAP/pair_sna_grid.cpp b/src/ML-SNAP/pair_sna_grid.cpp new file mode 100644 index 00000000000..75fba8bc39e --- /dev/null +++ b/src/ML-SNAP/pair_sna_grid.cpp @@ -0,0 +1,449 @@ +/* ---------------------------------------------------------------------- + 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 "pair_grid.h" +#include "pair_sna_grid.h" +#include "sna.h" +#include "../../build2/includes/lammps/atom.h" +#include "../../build/includes/lammps/update.h" +#include "../modify.h" +#include "../neighbor.h" +#include "../../build/includes/lammps/neigh_list.h" +#include "../neigh_request.h" +#include "../../build2/includes/lammps/force.h" +#include "../../build/includes/lammps/pair.h" +#include "../../build/includes/lammps/domain.h" +#include "../../build2/includes/lammps/comm.h" +#include "../memory.h" +#include "error.h" +#include "../tokenizer.h" + +#include +#include + +using namespace LAMMPS_NS; + +PairSNAGrid::PairSNAGrid(LAMMPS *lmp) : + PairGrid(lmp), + radelem(nullptr), wjelem(nullptr) +{ + snaptr = nullptr; +} + +/* ---------------------------------------------------------------------- */ + +PairSNAGrid::~PairSNAGrid() +{ + memory->destroy(radelem); + memory->destroy(wjelem); + memory->destroy(cutsq); + delete snaptr; + + if (chemflag) memory->destroy(map); +} + +/* ---------------------------------------------------------------------- */ + +void PairSNAGrid::init_style() +{ + if (force->newton_pair == 0) + error->all(FLERR,"Pair style sna/grid requires newton pair on"); + + snaptr = new SNA(lmp, rfac0, twojmax, + rmin0, switchflag, bzeroflag, + chemflag, bnormflag, wselfallflag, + nelements, switchinnerflag); + ncoeff = snaptr->ncoeff; + ndesc = ndesc_base + ncoeff; + snaptr->init(); + +} + +/* ---------------------------------------------------------------------- */ + +void PairSNAGrid::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void PairSNAGrid::compute(int eflag, int vflag) +{ + double fij[3]; + + ev_init(eflag,vflag); + + // compute sna for each gridpoint + + double** const x = atom->x; + double **f = atom->f; + const int* const mask = atom->mask; + int * const type = atom->type; + const int ntotal = atom->nlocal + atom->nghost; + + // insure rij, inside, and typej are of size ntotal + + snaptr->grow_rij(ntotal); + + // first generate fingerprint, + // which allows calculation of beta + + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + double xgrid[3]; + grid2x(ix, iy, iz, xgrid); + const double xtmp = xgrid[0]; + const double ytmp = xgrid[1]; + const double ztmp = xgrid[2]; + + // currently, all grid points are type 1 + + const int itype = 1; + int ielem = 0; + if (chemflag) + ielem = map[itype]; + const double radi = radelem[itype]; + + // 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++) { + + 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]; + int jelem = 0; + if (chemflag) + jelem = map[jtype]; + if (rsq < cutsq[jtype][jtype] && rsq > 1e-20) { + snaptr->rij[ninside][0] = delx; + snaptr->rij[ninside][1] = dely; + snaptr->rij[ninside][2] = delz; + snaptr->inside[ninside] = j; + snaptr->wj[ninside] = wjelem[jtype]; + snaptr->rcutij[ninside] = 2.0*radelem[jtype]*rcutfac; + if (chemflag) snaptr->element[ninside] = jelem; // element index for multi-element snap + ninside++; + } + } + + snaptr->compute_ui(ninside, ielem); + snaptr->compute_zi(); + snaptr->compute_bi(ielem); + + // linear contributions + + for (int icoeff = 0; icoeff < ncoeff; icoeff++) + gridlocal[ndesc_base+icoeff][iz][iy][ix] = + snaptr->blist[icoeff]; + + // quadratic contributions + // untested + + if (quadraticflag) { + int ncount = ncoeff; + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + double bveci = snaptr->blist[icoeff]; + gridlocal[ndesc_base+ncount++][iz][iy][ix] = + 0.5*bveci*bveci; + for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) + gridlocal[ndesc_base+ncount++][iz][iy][ix] = + bveci*snaptr->blist[jcoeff]; + } + } + } + + // this is a proxy for a call to the energy model + // beta is dE/dB^i, the derivative of the total + // energy w.r.t. to descriptors of grid point i + + compute_beta(); + + // second compute forces using beta + + int igrid = 0; + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + double xgrid[3]; + grid2x(ix, iy, iz, xgrid); + const double xtmp = xgrid[0]; + const double ytmp = xgrid[1]; + const double ztmp = xgrid[2]; + + // currently, all grid points are type 1 + + const int itype = 1; + int ielem = 0; + if (chemflag) + ielem = map[itype]; + const double radi = radelem[itype]; + + // 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++) { + + 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]; + int jelem = 0; + jelem = map[jtype]; + + if (rsq < cutsq[jtype][jtype] && rsq > 1e-20) { + snaptr->rij[ninside][0] = delx; + snaptr->rij[ninside][1] = dely; + snaptr->rij[ninside][2] = delz; + snaptr->inside[ninside] = j; + snaptr->wj[ninside] = wjelem[jtype]; + snaptr->rcutij[ninside] = 2.0*radelem[jtype]*rcutfac; + if (switchinnerflag) { + snaptr->sinnerij[ninside] = 0.5*(sinnerelem[ielem]+sinnerelem[jelem]); + snaptr->dinnerij[ninside] = 0.5*(dinnerelem[ielem]+dinnerelem[jelem]); + } + if (chemflag) snaptr->element[ninside] = jelem; + ninside++; + } + } + + // compute Ui, Yi for atom I + + if (chemflag) + snaptr->compute_ui(ninside, ielem); + else + snaptr->compute_ui(ninside, 0); + + // for neighbors of I within cutoff: + // compute Fij = dEi/dRj = -dEi/dRi + // add to Fi, subtract from Fj + // scaling is that for type I + + snaptr->compute_yi(beta[igrid]); + + for (int jj = 0; jj < ninside; jj++) { + int j = snaptr->inside[jj]; + snaptr->compute_duidrj(jj); + + snaptr->compute_deidrj(fij); + + f[j][0] += fij[0]; + f[j][1] += fij[1]; + f[j][2] += fij[2]; + + // tally per-atom virial contribution + + if (vflag) + ev_tally_xyz(-1,j,atom->nlocal,force->newton_pair,0.0,0.0, + fij[0],fij[1],fij[2], + -snaptr->rij[jj][0],-snaptr->rij[jj][1], + -snaptr->rij[jj][2]); + } + + // tally energy contribution + + if (eflag) { + + // get descriptors again + + snaptr->compute_zi(); + snaptr->compute_bi(ielem); + + // evdwl = energy of atom I, sum over coeffs_k * Bi_k + + double evdwl = 0.0; + + // E = beta.B + + for (int icoeff = 0; icoeff < ncoeff; icoeff++) + evdwl += beta[igrid][icoeff]*snaptr->blist[icoeff]; + + ev_tally_full(-1,2.0*evdwl,0.0,0.0,0.0,0.0,0.0); + + } + igrid++; + } + + if (vflag_fdotr) virial_fdotr_compute(); + +} + + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairSNAGrid::settings(int narg, char ** arg) +{ + + // call base class first + + PairGrid::settings(narg, arg); + + // 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 = 3+2*ntypes; + + if (narg < nargmin) error->all(FLERR,"Illegal pair sna/grid command"); + + // default values + + rmin0 = 0.0; + switchflag = 1; + bzeroflag = 1; + quadraticflag = 0; + chemflag = 0; + bnormflag = 0; + wselfallflag = 0; + switchinnerflag = 0; + nelements = 1; + + // process required arguments + + memory->create(radelem,ntypes+1,"pair:sna/grid:radelem"); // offset by 1 to match up with types + memory->create(wjelem,ntypes+1,"pair:sna/grid:wjelem"); + + rcutfac = atof(arg[0]); + rfac0 = atof(arg[1]); + twojmax = atoi(arg[2]); + + for(int i = 0; i < ntypes; i++) + radelem[i+1] = atof(arg[3+i]); + for(int i = 0; i < ntypes; i++) + wjelem[i+1] = atof(arg[3+ntypes+i]); + + // construct cutsq + + double cut; + cutmax = 0.0; + memory->create(cutsq,ntypes+1,ntypes+1,"pair: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; + } + } + + // set local input checks + + int sinnerflag = 0; + int dinnerflag = 0; + + // process optional args + + int iarg = nargmin; + + while (iarg < narg) { + if (strcmp(arg[iarg],"rmin0") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal pair 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 pair 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 pair 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 pair sna/grid command"); + quadraticflag = atoi(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"chem") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal pair sna/grid command"); + chemflag = 1; + memory->create(map,ntypes+1,"pair:sna/grid:map"); + nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + for (int i = 0; i < ntypes; i++) { + int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp); + if (jelem < 0 || jelem >= nelements) + error->all(FLERR,"Illegal pair sna/grid command"); + map[i+1] = jelem; + } + iarg += 2+ntypes; + } else if (strcmp(arg[iarg],"bnormflag") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal pair sna/grid command"); + bnormflag = atoi(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"wselfallflag") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal pair sna/grid command"); + wselfallflag = atoi(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal pair sna/grid command"); + switchinnerflag = atoi(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"sinner") == 0) { + iarg++; + if (iarg+ntypes > narg) + error->all(FLERR,"Illegal pair sna/grid command"); + memory->create(sinnerelem,ntypes+1,"snap:sinnerelem"); + for (int i = 0; i < ntypes; i++) + sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerflag = 1; + iarg += ntypes; + } else if (strcmp(arg[iarg],"dinner") == 0) { + iarg++; + if (iarg+ntypes > narg) + error->all(FLERR,"Illegal pair sna/grid command"); + memory->create(dinnerelem,ntypes+1,"snap:dinnerelem"); + for (int i = 0; i < ntypes; i++) + dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerflag = 1; + iarg += ntypes; + } else error->all(FLERR,"Illegal pair sna/grid command"); + + } + +} + +/* ---------------------------------------------------------------------- + memory usage +------------------------------------------------------------------------- */ + +double PairSNAGrid::memory_usage() +{ + double nbytes = snaptr->memory_usage(); // SNA object + int n = atom->ntypes+1; + nbytes += (double)n*sizeof(int); // map + + return nbytes; +} + diff --git a/src/ML-SNAP/pair_sna_grid.h b/src/ML-SNAP/pair_sna_grid.h new file mode 100644 index 00000000000..653bb9d1c13 --- /dev/null +++ b/src/ML-SNAP/pair_sna_grid.h @@ -0,0 +1,79 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 PAIR_CLASS +// clang-format off +PairStyle(sna/grid, PairSNAGrid); +// clang-format on +#else + +#ifndef LMP_PAIR_SNA_GRID_H +#define LMP_PAIR_SNA_GRID_H + +#include "pair_grid.h" + +namespace LAMMPS_NS { + +class PairSNAGrid : public PairGrid { + public: + PairSNAGrid(class LAMMPS *); + ~PairSNAGrid(); + + void init_style(); + void init_list(int, class NeighList *); + void settings(int, char **); + void compute(int, int); + double memory_usage(); + + private: + int ncoeff; + class NeighList *list; + double rcutfac; + double *radelem; + double *wjelem; + class SNA *snaptr; + int quadraticflag; + int twojmax, switchflag, bzeroflag, bnormflag; + int chemflag, wselfallflag; + int switchinnerflag; + double *sinnerelem; + double *dinnerelem; + double rfac0, rmin0; +}; + +} + +#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/local requires a pair style be defined + +Self-explanatory. + +E: Compute sna/grid/local cutoff is longer than pairwise cutoff + +Self-explanatory. + +W: More than one compute sna/grid/local + +Self-explanatory. + +*/