First successful compile of pair style

This commit is contained in:
Aidan Thompson
2021-10-03 17:55:55 -06:00
parent 0a4e85a1f3
commit 8bddc801df
4 changed files with 690 additions and 0 deletions

View File

@ -0,0 +1,182 @@
/* ----------------------------------------------------------------------
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 "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 "domain.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "tokenizer.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
PairSNAGrid::PairSNAGrid(LAMMPS *lmp) :
PairGrid(lmp), cutsq(nullptr),
radelem(nullptr), wjelem(nullptr)
{
}
/* ---------------------------------------------------------------------- */
PairSNAGrid::~PairSNAGrid()
{
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(cutsq);
delete snaptr;
if (chemflag) memory->destroy(map);
}
/* ---------------------------------------------------------------------- */
void PairSNAGrid::init()
{
if (force->pair == nullptr)
error->all(FLERR,"Pair sna/grid requires a pair style be defined");
if (cutmax > force->pair->cutforce)
error->all(FLERR,"Pair 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;
snaptr->init();
}
/* ---------------------------------------------------------------------- */
void PairSNAGrid::init_list(int /*id*/, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void PairSNAGrid::compute(int eflag, int vflag)
{
// compute sna for each gridpoint
double** const x = atom->x;
const int* const mask = atom->mask;
int * const type = atom->type;
const int ntotal = atom->nlocal + atom->nghost;
// insure rij, inside, and typej are of size jnum
snaptr->grow_rij(ntotal);
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;
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
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];
}
}
}
}
/* ----------------------------------------------------------------------
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;
}

View File

@ -0,0 +1,76 @@
/* -*- 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();
void init_list(int, class NeighList *);
void compute(int, int);
double memory_usage();
private:
int ncoeff;
double **cutsq;
class NeighList *list;
double rcutfac;
double *radelem;
double *wjelem;
int *map; // map types to [0,nelements)
int nelements, chemflag;
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/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.
*/

344
src/pair_grid.cpp Normal file
View File

@ -0,0 +1,344 @@
/* ----------------------------------------------------------------------
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 <mpi.h>
#include <cstring>
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "comm.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairGrid::PairGrid(LAMMPS *lmp) :
Pair(lmp), gridlocal(nullptr), alocal(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;
memory->destroy(beta);
deallocate_grid();
}
/* ---------------------------------------------------------------------- */
void PairGrid::init()
{
}
/* ---------------------------------------------------------------------- */
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");
}
}
/* ----------------------------------------------------------------------
free arrays
------------------------------------------------------------------------- */
void PairGrid::deallocate_grid()
{
if (gridlocal_allocated) {
gridlocal_allocated = 0;
memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo);
memory->destroy(alocal);
}
}
/* ----------------------------------------------------------------------
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<int> (xfraclo * nx);
if (1.0*nxlo != xfraclo*nx) nxlo++;
nxhi = static_cast<int> (xfrachi * nx);
if (1.0*nxhi == xfrachi*nx) nxhi--;
nylo = static_cast<int> (yfraclo * ny);
if (1.0*nylo != yfraclo*ny) nylo++;
nyhi = static_cast<int> (yfrachi * ny);
if (1.0*nyhi == yfrachi*ny) nyhi--;
nzlo = static_cast<int> (zfraclo * nz);
if (1.0*nzlo != zfraclo*nz) nzlo++;
nzhi = static_cast<int> (zfrachi * nz);
if (1.0*nzhi == zfrachi*nz) nzhi--;
ngridlocal = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1);
}
/* ----------------------------------------------------------------------
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++;
}
}
/* ----------------------------------------------------------------------
get beta from someplace
------------------------------------------------------------------------- */
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 = ndesc_base; icol < ndesc; icol++)
beta[igrid][icol] = 1.0;
igrid++;
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairGrid::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
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 < 6) error->all(FLERR,"Illegal pair style command");
int iarg0 = 3;
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 + atom->ntypes) error->all(FLERR,"Incorrect args for pair coefficients");
map_element2type(narg-4,arg+4);
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairGrid::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style grid requires newton pair on");
// no neighbor list
// int irequest = neighbor->request(this,instance_me);
// neighbor->requests[irequest]->half = 0;
// neighbor->requests[irequest]->full = 1;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairGrid::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
return 0.0;
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double PairGrid::memory_usage()
{
int nbytes = ndesc*ngridlocal*sizeof(double); // gridlocal
return nbytes;
}

88
src/pair_grid.h Normal file
View File

@ -0,0 +1,88 @@
/* -*- 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();
void init();
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 void init_style();
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
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.
*/