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

Completed merge to current master

Browse Source
This commit is contained in:
Aidan Thompson
2021-07-14 13:47:55 -06:00
parent 38b156a08a 07db7a4095
commit 2cf00a382e
8 changed files with 1082 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

264
src/SNAP/compute_sna_grid.cpp Normal file
View File

@ -0,0 +1,264 @@
/* ----------------------------------------------------------------------
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 "compute_sna_grid.h"
#include <cstring>
#include <cstdlib>
#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) :
ComputeGrid(lmp, narg, arg), cutsq(NULL), sna(NULL),
radelem(NULL), wjelem(NULL)
{
double rmin0, rfac0;
int twojmax, switchflag, bzeroflag;
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;
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;
nvalues = ncoeff;
if (quadraticflag) nvalues += (ncoeff*(ncoeff+1))/2;
size_array_cols = size_array_cols_base + nvalues;
array_flag = 1;
}
/* ---------------------------------------------------------------------- */
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::compute_array()
{
invoked_array = update->ntimestep;
int * const type = atom->type;
// compute sna for each gridpoint
double** const x = atom->x;
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);
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
int ninside = 0;
for (int j = 0; j < ntotal; j++) {
// check that j is in compute group
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) {
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++)
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
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++)
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
------------------------------------------------------------------------- */
double ComputeSNAGrid::memory_usage()
{
double nbytes = snaptr->memory_usage(); // SNA object
return nbytes;
}

71
src/SNAP/compute_sna_grid.h Normal file
View File

@ -0,0 +1,71 @@
/* -*- 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_grid.h"
namespace LAMMPS_NS {
class ComputeSNAGrid : public ComputeGrid {
public:
ComputeSNAGrid(class LAMMPS *, int, char **);
~ComputeSNAGrid();
void init();
void compute_array();
double memory_usage();
private:
int ncoeff;
double **cutsq;
double **sna;
double rcutfac;
double *radelem;
double *wjelem;
class SNA* snaptr;
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.
*/

289
src/compute_grid.cpp Normal file
View File

@ -0,0 +1,289 @@
/* ----------------------------------------------------------------------
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 <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;
/* ---------------------------------------------------------------------- */
ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), grid(NULL), local_flags(NULL), gridlocal(NULL)
{
if (narg < 6) error->all(FLERR,"Illegal compute grid command");
array_flag = 1;
size_array_cols = 0;
size_array_rows = 0;
extarray = 0;
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");
nargbase = iarg - iarg0;
size_array_rows = ngrid = nx*ny*nz;
size_array_cols_base = 3;
gridlocal_allocated = 0;
}
/* ---------------------------------------------------------------------- */
ComputeGrid::~ComputeGrid()
{
memory->destroy(grid);
memory->destroy(local_flags);
if (gridlocal_allocated)
memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo);
}
/* ---------------------------------------------------------------------- */
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);
}
/* ----------------------------------------------------------------------
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)
{
int ix, iy, iz;
grid2ix(igrid, ix, iy, iz);
int islocal =
ix >= nxlo && ix <= nxhi &&
iy >= nylo && iy <= nyhi &&
iz >= nzlo && iz <= nzhi;
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()
{
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 * sizeof(double));
}
}
}
/* ----------------------------------------------------------------------
free and reallocate arrays
------------------------------------------------------------------------- */
void ComputeGrid::allocate()
{
// allocate arrays
memory->destroy(grid);
memory->destroy(local_flags);
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");
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;
}
/* ----------------------------------------------------------------------
set global grid
------------------------------------------------------------------------- */
void ComputeGrid::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 ComputeGrid::set_grid_local()
{
// 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
nxlo = static_cast<int> (comm->xsplit[comm->myloc[0]] * nx);
nxhi = static_cast<int> (comm->xsplit[comm->myloc[0]+1] * nx) - 1;
nylo = static_cast<int> (comm->ysplit[comm->myloc[1]] * ny);
nyhi = static_cast<int> (comm->ysplit[comm->myloc[1]+1] * ny) - 1;
nzlo = static_cast<int> (comm->zsplit[comm->myloc[2]] * nz);
nzhi = static_cast<int> (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);
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeGrid::memory_usage()
{
double nbytes = size_array_rows*size_array_cols *
sizeof(double); // grid
nbytes += size_array_rows*sizeof(int); // local_flags
nbytes += size_array_cols*ngridlocal*sizeof(double); // gridlocal
return nbytes;
}

75
src/compute_grid.h Normal file
View File

@ -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.
------------------------------------------------------------------------- */
#ifndef LMP_COMPUTE_GRID_H
#define LMP_COMPUTE_GRID_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeGrid : public Compute {
public:
ComputeGrid(class LAMMPS *, int, char **);
virtual ~ComputeGrid();
void init();
void setup();
virtual void compute_array() = 0;
double memory_usage();
protected:
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)
double delxinv,delyinv,delzinv; // inverse grid spacing
double delx,dely,delz; // grid spacing
int nargbase; // number of base class args
double cutmax; // largest cutoff distance
int size_array_cols_base; // number of columns used for coords, etc.
int *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
void set_grid_global(); // set global grid
void set_grid_local(); // set bounds for local grid
private:
};
}
#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.
*/