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

Merge pull request #3305 from athomps/compute-grid-new

Browse Source 
Compute grid for ML-SNAP
This commit is contained in:
Axel Kohlmeyer
2022-06-28 09:00:38 -04:00
committed by GitHub
parent ae18e1e01c a9f3108f29
commit f65bc76326
29 changed files with 2732 additions and 358 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
src/.gitignore vendored
View File

@ -173,12 +173,20 @@
/pair_tdpd.cpp
/pair_tdpd.h
/compute_grid.cpp
/compute_grid.h
/compute_grid_local.cpp
/compute_grid_local.h
/compute_sna_atom.cpp
/compute_sna_atom.h
/compute_snad_atom.cpp
/compute_snad_atom.h
/compute_snav_atom.cpp
/compute_snav_atom.h
/compute_sna_grid.cpp
/compute_sna_grid.h
/compute_sna_grid_local.cpp
/compute_sna_grid_local.h
/compute_snap.cpp
/compute_snap.h
/openmp_snap.h

242
src/ML-SNAP/compute_grid.cpp Normal file
View File

@ -0,0 +1,242 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "compute_grid.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), grid(nullptr), gridall(nullptr), gridlocal(nullptr)
{
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 = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
ny = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
nz = utils::inumeric(FLERR, arg[iarg + 3], false, lmp);
if (nx <= 0 || ny <= 0 || nz <= 0) error->all(FLERR, "All grid dimensions must be positive");
iarg += 4;
} else
error->all(FLERR, "Illegal compute grid command");
nargbase = iarg - iarg0;
size_array_rows = nx * ny * nz;
size_array_cols_base = 3;
gridlocal_allocated = 0;
}
/* ---------------------------------------------------------------------- */
ComputeGrid::~ComputeGrid()
{
deallocate();
}
/* ---------------------------------------------------------------------- */
void ComputeGrid::setup()
{
deallocate();
set_grid_global();
set_grid_local();
allocate();
}
/* ----------------------------------------------------------------------
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);
}
/* ----------------------------------------------------------------------
copy coords to global array
------------------------------------------------------------------------- */
void ComputeGrid::assign_coords_all()
{
double x[3];
for (int igrid = 0; igrid < size_array_rows; igrid++) {
grid2x(igrid, x);
gridall[igrid][0] = x[0];
gridall[igrid][1] = x[1];
gridall[igrid][2] = x[2];
}
}
/* ----------------------------------------------------------------------
create arrays
------------------------------------------------------------------------- */
void ComputeGrid::allocate()
{
// allocate arrays
memory->create(grid, size_array_rows, size_array_cols, "grid:grid");
memory->create(gridall, size_array_rows, size_array_cols, "grid:gridall");
if (nxlo <= nxhi && nylo <= nyhi && nzlo <= nzhi) {
gridlocal_allocated = 1;
memory->create4d_offset(gridlocal, size_array_cols, nzlo, nzhi, nylo, nyhi, nxlo, nxhi,
"grid:gridlocal");
}
array = gridall;
}
/* ----------------------------------------------------------------------
free arrays
------------------------------------------------------------------------- */
void ComputeGrid::deallocate()
{
memory->destroy(grid);
memory->destroy(gridall);
if (gridlocal_allocated) {
gridlocal_allocated = 0;
memory->destroy4d_offset(gridlocal, nzlo, nylo, nxlo);
}
array = nullptr;
}
/* ----------------------------------------------------------------------
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()
{
// 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);
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeGrid::memory_usage()
{
double nbytes = size_array_rows * size_array_cols * sizeof(double); // grid
nbytes += size_array_rows * size_array_cols * sizeof(double); // gridall
nbytes += size_array_cols * ngridlocal * sizeof(double); // gridlocal
return nbytes;
}

58
src/ML-SNAP/compute_grid.h Normal file
View File

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

270
src/ML-SNAP/compute_grid_local.cpp Normal file
View File

@ -0,0 +1,270 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "compute_grid_local.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include <cstring>
// For the subdomain test below; grid-points and subdomain boundaries
// sometimes differ by minimal amounts (in the order of 2e-17).
static constexpr double EPSILON = 1.0e-10;
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeGridLocal::ComputeGridLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), alocal(nullptr)
{
if (narg < 6) error->all(FLERR, "Illegal compute grid/local command");
local_flag = 1;
size_local_cols = 0;
size_local_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/local command");
nx = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
ny = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
nz = utils::inumeric(FLERR, arg[iarg + 3], false, lmp);
if (nx <= 0 || ny <= 0 || nz <= 0)
error->all(FLERR, "All grid/local dimensions must be positive");
iarg += 4;
} else
error->all(FLERR, "Illegal compute grid/local command");
nargbase = iarg - iarg0;
size_local_cols_base = 6;
gridlocal_allocated = 0;
}
/* ---------------------------------------------------------------------- */
ComputeGridLocal::~ComputeGridLocal()
{
deallocate();
}
/* ---------------------------------------------------------------------- */
void ComputeGridLocal::setup()
{
deallocate();
set_grid_global();
set_grid_local();
allocate();
assign_coords();
}
/* ----------------------------------------------------------------------
convert global array indexes to box coords
------------------------------------------------------------------------- */
void ComputeGridLocal::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);
}
/* ----------------------------------------------------------------------
convert global array indexes to lamda coords; for orthorombic
cells defaults to grid2x.
------------------------------------------------------------------------- */
void ComputeGridLocal::grid2lamda(int ix, int iy, int iz, double *x)
{
x[0] = ix * delx;
x[1] = iy * dely;
x[2] = iz * delz;
}
/* ----------------------------------------------------------------------
create arrays
------------------------------------------------------------------------- */
void ComputeGridLocal::allocate()
{
if (nxlo <= nxhi && nylo <= nyhi && nzlo <= nzhi) {
gridlocal_allocated = 1;
memory->create(alocal, size_local_rows, size_local_cols, "compute/grid/local:alocal");
array_local = alocal;
}
}
/* ----------------------------------------------------------------------
free arrays
------------------------------------------------------------------------- */
void ComputeGridLocal::deallocate()
{
if (gridlocal_allocated) {
gridlocal_allocated = 0;
memory->destroy(alocal);
}
array_local = nullptr;
}
/* ----------------------------------------------------------------------
set global grid
------------------------------------------------------------------------- */
void ComputeGridLocal::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 ComputeGridLocal::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--;
size_local_rows = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1);
}
/* ----------------------------------------------------------------------
copy coords to local array
------------------------------------------------------------------------- */
void ComputeGridLocal::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];
// for triclinic: create gridpoint in lamda coordinates and transform after check.
// for orthorombic: create gridpoint in box coordinates.
if (triclinic)
grid2lamda(ix, iy, iz, xgrid);
else
grid2x(ix, iy, iz, xgrid);
// ensure gridpoint is not strictly outside subdomain
if ((sublo[0] - xgrid[0]) > EPSILON || (xgrid[0] - subhi[0]) > EPSILON ||
(sublo[1] - xgrid[1]) > EPSILON || (xgrid[1] - subhi[1]) > EPSILON ||
(sublo[2] - xgrid[2]) > EPSILON || (xgrid[2] - subhi[2]) > EPSILON)
error->one(FLERR, "Invalid gridpoint position in compute grid/local");
// convert lamda to x, y, z, after sudomain check
if (triclinic) domain->lamda2x(xgrid, xgrid);
alocal[igrid][3] = xgrid[0];
alocal[igrid][4] = xgrid[1];
alocal[igrid][5] = xgrid[2];
igrid++;
}
}
/* ----------------------------------------------------------------------
memory usage of local data
------------------------------------------------------------------------- */
double ComputeGridLocal::memory_usage()
{
int nbytes = size_local_rows * size_local_cols * sizeof(double); // gridlocal
return nbytes;
}

56
src/ML-SNAP/compute_grid_local.h Normal file
View File

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

151
src/ML-SNAP/compute_sna_atom.cpp
View File

@ -35,20 +35,21 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr)
{
double rmin0, rfac0;
// begin code common to all SNAP computes
double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
int ntypes = atom->ntypes;
int nargmin = 6+2*ntypes;
int nargmin = 6 + 2 * ntypes;
if (narg < nargmin) error->all(FLERR,"Illegal compute sna/atom command");
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
// default values
rmin0 = 0.0;
switchflag = 1;
bzeroflag = 1;
bnormflag = 0;
quadraticflag = 0;
chemflag = 0;
bnormflag = 0;
@ -56,32 +57,34 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
switchinnerflag = 0;
nelements = 1;
// offset by 1 to match up with types
// process required arguments
memory->create(radelem,ntypes+1,"sna/atom:radelem");
memory->create(wjelem,ntypes+1,"sna/atom:wjelem");
memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types
memory->create(wjelem, ntypes + 1, "sna/atom:wjelem");
rcutfac = atof(arg[3]);
rfac0 = atof(arg[4]);
twojmax = atoi(arg[5]);
rcutfac = utils::numeric(FLERR, arg[3], false, lmp);
rfac0 = utils::numeric(FLERR, arg[4], false, lmp);
twojmax = utils::inumeric(FLERR, arg[5], false, lmp);
for (int i = 0; i < ntypes; i++)
radelem[i+1] = atof(arg[6+i]);
radelem[i + 1] =
utils::numeric(FLERR, arg[6 + i], false, lmp);
for (int i = 0; i < ntypes; i++)
wjelem[i+1] = atof(arg[6+ntypes+i]);
wjelem[i + 1] =
utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp);
// construct cutsq
double cut;
cutmax = 0.0;
memory->create(cutsq,ntypes+1,ntypes+1,"sna/atom:cutsq");
memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq");
for (int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
cut = 2.0 * radelem[i] * rcutfac;
if (cut > cutmax) cutmax = cut;
cutsq[i][i] = cut*cut;
for (int j = i+1; j <= ntypes; j++) {
cut = (radelem[i]+radelem[j])*rcutfac;
cutsq[i][j] = cutsq[j][i] = cut*cut;
cutsq[i][i] = cut * cut;
for (int j = i + 1; j <= ntypes; j++) {
cut = (radelem[i] + radelem[j]) * rcutfac;
cutsq[i][j] = cutsq[j][i] = cut * cut;
}
}
@ -95,89 +98,87 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
rmin0 = atof(arg[iarg+1]);
if (strcmp(arg[iarg], "rmin0") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"switchflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
switchflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"bzeroflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
bzeroflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "bzeroflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"quadraticflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
quadraticflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "quadraticflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2+ntypes > narg)
error->all(FLERR,"Illegal compute sna/atom command");
} else if (strcmp(arg[iarg], "chem") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
chemflag = 1;
memory->create(map,ntypes+1,"compute_sna_atom:map");
nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
memory->create(map, ntypes + 1, "compute_sna_grid:map");
nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
for (int i = 0; i < ntypes; i++) {
int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp);
if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute sna/atom command");
map[i+1] = jelem;
int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp);
if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style);
map[i + 1] = jelem;
}
iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2 + ntypes;
} else if (strcmp(arg[iarg], "bnormflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
wselfallflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "wselfallflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
switchinnerflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchinnerflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"sinner") == 0) {
} else if (strcmp(arg[iarg], "sinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute sna/atom command");
memory->create(sinnerelem,ntypes+1,"sna/atom:sinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style);
memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem");
for (int i = 0; i < ntypes; i++)
sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
sinnerflag = 1;
iarg += ntypes;
} else if (strcmp(arg[iarg],"dinner") == 0) {
} else if (strcmp(arg[iarg], "dinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute sna/atom command");
memory->create(dinnerelem,ntypes+1,"sna/atom:dinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style);
memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem");
for (int i = 0; i < ntypes; i++)
dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
dinnerflag = 1;
iarg += ntypes;
} else error->all(FLERR,"Illegal compute sna/atom command");
} else
error->all(FLERR, "Illegal compute {} command", style);
}
if (switchinnerflag && !(sinnerflag && dinnerflag))
error->all(FLERR,"Illegal compute sna/atom command: switchinnerflag = 1, missing sinner/dinner keyword");
error->all(
FLERR,
"Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword",
style);
if (!switchinnerflag && (sinnerflag || dinnerflag))
error->all(FLERR,"Illegal compute sna/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword");
error->all(
FLERR,
"Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword",
style);
snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag,
chemflag, bnormflag, wselfallflag,
nelements, switchinnerflag);
snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag,
wselfallflag, nelements, switchinnerflag);
ncoeff = snaptr->ncoeff;
size_peratom_cols = ncoeff;
if (quadraticflag) size_peratom_cols += (ncoeff*(ncoeff+1))/2;
nvalues = ncoeff;
if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2;
// end code common to all SNAP computes
size_peratom_cols = nvalues;
peratom_flag = 1;
nmax = 0;

1
src/ML-SNAP/compute_sna_atom.h
View File

@ -50,6 +50,7 @@ class ComputeSNAAtom : public Compute {
class SNA *snaptr;
double cutmax;
int quadraticflag;
int nvalues;
};
} // namespace LAMMPS_NS

320
src/ML-SNAP/compute_sna_grid.cpp Normal file
View File

@ -0,0 +1,320 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "compute_sna_grid.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "sna.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) :
ComputeGrid(lmp, narg, arg), cutsq(nullptr), radelem(nullptr), wjelem(nullptr)
{
// skip over arguments used by base class
// so that argument positions are identical to
// regular per-atom compute
arg += nargbase;
narg -= nargbase;
// begin code common to all SNAP computes
double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
int ntypes = atom->ntypes;
int nargmin = 6 + 2 * ntypes;
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
// 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, "sna/atom:radelem"); // offset by 1 to match up with types
memory->create(wjelem, ntypes + 1, "sna/atom:wjelem");
rcutfac = utils::numeric(FLERR, arg[3], false, lmp);
rfac0 = utils::numeric(FLERR, arg[4], false, lmp);
twojmax = utils::inumeric(FLERR, arg[5], false, lmp);
for (int i = 0; i < ntypes; i++) radelem[i + 1] = utils::numeric(FLERR, arg[6 + i], false, lmp);
for (int i = 0; i < ntypes; i++)
wjelem[i + 1] = utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp);
// construct cutsq
double cut;
cutmax = 0.0;
memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom: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 compute {} command", style);
rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "switchflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "bzeroflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "quadraticflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "chem") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
chemflag = 1;
memory->create(map, ntypes + 1, "compute_sna_grid:map");
nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
for (int i = 0; i < ntypes; i++) {
int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp);
if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style);
map[i + 1] = jelem;
}
iarg += 2 + ntypes;
} else if (strcmp(arg[iarg], "bnormflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "wselfallflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "switchinnerflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "sinner") == 0) {
iarg++;
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style);
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 compute {} command", style);
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 compute {} command", style);
}
if (switchinnerflag && !(sinnerflag && dinnerflag))
error->all(FLERR,
"Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword",
style);
if (!switchinnerflag && (sinnerflag || dinnerflag))
error->all(FLERR,
"Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword",
style);
snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag,
wselfallflag, nelements, switchinnerflag);
ncoeff = snaptr->ncoeff;
nvalues = ncoeff;
if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2;
// end code common to all SNAP computes
size_array_cols = size_array_cols_base + nvalues;
array_flag = 1;
}
/* ---------------------------------------------------------------------- */
ComputeSNAGrid::~ComputeSNAGrid()
{
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(cutsq);
delete snaptr;
if (chemflag) memory->destroy(map);
}
/* ---------------------------------------------------------------------- */
void ComputeSNAGrid::init()
{
if ((modify->get_compute_by_style("^sna/grid$").size() > 1) && (comm->me == 0))
error->warning(FLERR, "More than one instance of compute sna/grid");
snaptr->init();
}
/* ---------------------------------------------------------------------- */
void ComputeSNAGrid::compute_array()
{
invoked_array = update->ntimestep;
// 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];
const int igrid = iz * (nx * ny) + iy * nx + ix;
grid2x(igrid, xgrid);
const double xtmp = xgrid[0];
const double ytmp = xgrid[1];
const double ztmp = xgrid[2];
// currently, all grid points are type 1
// not clear what a better choice would be
const int itype = 1;
int ielem = 0;
if (chemflag) ielem = map[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++) {
// 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];
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 (switchinnerflag) {
snaptr->sinnerij[ninside] = sinnerelem[jelem];
snaptr->dinnerij[ninside] = dinnerelem[jelem];
}
if (chemflag) snaptr->element[ninside] = jelem;
ninside++;
}
}
snaptr->compute_ui(ninside, ielem);
snaptr->compute_zi();
snaptr->compute_bi(ielem);
// linear contributions
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
gridlocal[size_array_cols_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[size_array_cols_base + ncount++][iz][iy][ix] = 0.5 * bveci * bveci;
for (int jcoeff = icoeff + 1; jcoeff < ncoeff; jcoeff++)
gridlocal[size_array_cols_base + ncount++][iz][iy][ix] =
bveci * snaptr->blist[jcoeff];
}
}
}
memset(&grid[0][0], 0, size_array_rows * size_array_cols * sizeof(double));
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], size_array_rows * size_array_cols, MPI_DOUBLE, MPI_SUM,
world);
assign_coords_all();
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double ComputeSNAGrid::memory_usage()
{
double nbytes = snaptr->memory_usage(); // SNA object
int n = atom->ntypes + 1;
nbytes += (double) n * sizeof(int); // map
return nbytes;
}

54
src/ML-SNAP/compute_sna_grid.h Normal file
View File

@ -0,0 +1,54 @@
/* -*- 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 COMPUTE_CLASS
// clang-format off
ComputeStyle(sna/grid,ComputeSNAGrid);
// clang-format on
#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() override;
void init() override;
void compute_array() override;
double memory_usage() override;
private:
int ncoeff;
double **cutsq;
double rcutfac;
double *radelem;
double *wjelem;
int *map; // map types to [0,nelements)
int nelements, chemflag;
int switchinnerflag;
double *sinnerelem;
double *dinnerelem;
class SNA *snaptr;
double cutmax;
int quadraticflag;
};
} // namespace LAMMPS_NS
#endif
#endif

306
src/ML-SNAP/compute_sna_grid_local.cpp Normal file
View File

@ -0,0 +1,306 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "compute_sna_grid_local.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "sna.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
ComputeSNAGridLocal::ComputeSNAGridLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeGridLocal(lmp, narg, arg), cutsq(nullptr), radelem(nullptr), wjelem(nullptr)
{
// skip over arguments used by base class
// so that argument positions are identical to
// regular per-atom compute
arg += nargbase;
narg -= nargbase;
// begin code common to all SNAP computes
double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
int ntypes = atom->ntypes;
int nargmin = 6 + 2 * ntypes;
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
// 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, "sna/atom:radelem"); // offset by 1 to match up with types
memory->create(wjelem, ntypes + 1, "sna/atom:wjelem");
rcutfac = utils::numeric(FLERR, arg[3], false, lmp);
rfac0 = utils::numeric(FLERR, arg[4], false, lmp);
twojmax = utils::inumeric(FLERR, arg[5], false, lmp);
for (int i = 0; i < ntypes; i++) radelem[i + 1] = utils::numeric(FLERR, arg[6 + i], false, lmp);
for (int i = 0; i < ntypes; i++)
wjelem[i + 1] = utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp);
// construct cutsq
double cut;
cutmax = 0.0;
memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom: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 compute {} command", style);
rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "switchflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "bzeroflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "quadraticflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "chem") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
chemflag = 1;
memory->create(map, ntypes + 1, "compute_sna_grid:map");
nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
for (int i = 0; i < ntypes; i++) {
int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp);
if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style);
map[i + 1] = jelem;
}
iarg += 2 + ntypes;
} else if (strcmp(arg[iarg], "bnormflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "wselfallflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "switchinnerflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "sinner") == 0) {
iarg++;
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style);
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 compute {} command", style);
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 compute {} command", style);
}
if (switchinnerflag && !(sinnerflag && dinnerflag))
error->all(FLERR,
"Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword",
style);
if (!switchinnerflag && (sinnerflag || dinnerflag))
error->all(FLERR,
"Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword",
style);
snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag,
wselfallflag, nelements, switchinnerflag);
ncoeff = snaptr->ncoeff;
nvalues = ncoeff;
if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2;
// end code common to all SNAP computes
size_local_cols = size_local_cols_base + nvalues;
}
/* ---------------------------------------------------------------------- */
ComputeSNAGridLocal::~ComputeSNAGridLocal()
{
memory->destroy(radelem);
memory->destroy(wjelem);
memory->destroy(cutsq);
delete snaptr;
if (chemflag) memory->destroy(map);
}
/* ---------------------------------------------------------------------- */
void ComputeSNAGridLocal::init()
{
if ((modify->get_compute_by_style("^sna/grid/local$").size() > 1) && (comm->me == 0))
error->warning(FLERR, "More than one instance of compute sna/grid/local");
snaptr->init();
}
/* ---------------------------------------------------------------------- */
void ComputeSNAGridLocal::compute_local()
{
invoked_array = update->ntimestep;
// 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);
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
// not clear what a better choice would be
const int itype = 1;
int ielem = 0;
if (chemflag) ielem = map[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++) {
// 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];
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 (switchinnerflag) {
snaptr->sinnerij[ninside] = sinnerelem[jelem];
snaptr->dinnerij[ninside] = dinnerelem[jelem];
}
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++)
alocal[igrid][size_local_cols_base + icoeff] = snaptr->blist[icoeff];
// quadratic contributions
if (quadraticflag) {
int ncount = ncoeff;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bveci = snaptr->blist[icoeff];
alocal[igrid][size_local_cols_base + ncount++] = 0.5 * bveci * bveci;
for (int jcoeff = icoeff + 1; jcoeff < ncoeff; jcoeff++)
alocal[igrid][size_local_cols_base + ncount++] = bveci * snaptr->blist[jcoeff];
}
}
igrid++;
}
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double ComputeSNAGridLocal::memory_usage()
{
double nbytes = snaptr->memory_usage(); // SNA object
int n = atom->ntypes + 1;
nbytes += (double) n * sizeof(int); // map
return nbytes;
}

54
src/ML-SNAP/compute_sna_grid_local.h Normal file
View File

@ -0,0 +1,54 @@
/* -*- 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 COMPUTE_CLASS
// clang-format off
ComputeStyle(sna/grid/local,ComputeSNAGridLocal);
// clang-format on
#else
#ifndef LMP_COMPUTE_SNA_GRID_LOCAL_H
#define LMP_COMPUTE_SNA_GRID_LOCAL_H
#include "compute_grid_local.h"
namespace LAMMPS_NS {
class ComputeSNAGridLocal : public ComputeGridLocal {
public:
ComputeSNAGridLocal(class LAMMPS *, int, char **);
~ComputeSNAGridLocal() override;
void init() override;
void compute_local() override;
double memory_usage() override;
private:
int ncoeff;
double **cutsq;
double rcutfac;
double *radelem;
double *wjelem;
int *map; // map types to [0,nelements)
int nelements, chemflag;
int switchinnerflag;
double *sinnerelem;
double *dinnerelem;
class SNA *snaptr;
double cutmax;
int quadraticflag;
};
} // namespace LAMMPS_NS
#endif
#endif

158
src/ML-SNAP/compute_snad_atom.cpp
View File

@ -34,20 +34,22 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snad(nullptr),
radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr)
{
// begin code common to all SNAP computes
double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
int ntypes = atom->ntypes;
int nargmin = 6+2*ntypes;
int nargmin = 6 + 2 * ntypes;
if (narg < nargmin) error->all(FLERR,"Illegal compute snad/atom command");
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
// default values
rmin0 = 0.0;
switchflag = 1;
bzeroflag = 1;
bnormflag = 0;
quadraticflag = 0;
chemflag = 0;
bnormflag = 0;
@ -57,28 +59,32 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
// process required arguments
memory->create(radelem,ntypes+1,"snad/atom:radelem"); // offset by 1 to match up with types
memory->create(wjelem,ntypes+1,"snad/atom:wjelem");
rcutfac = atof(arg[3]);
rfac0 = atof(arg[4]);
twojmax = atoi(arg[5]);
memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types
memory->create(wjelem, ntypes + 1, "sna/atom:wjelem");
rcutfac = utils::numeric(FLERR, arg[3], false, lmp);
rfac0 = utils::numeric(FLERR, arg[4], false, lmp);
twojmax = utils::inumeric(FLERR, arg[5], false, lmp);
for (int i = 0; i < ntypes; i++)
radelem[i+1] = atof(arg[6+i]);
radelem[i + 1] =
utils::numeric(FLERR, arg[6 + i], false, lmp);
for (int i = 0; i < ntypes; i++)
wjelem[i+1] = atof(arg[6+ntypes+i]);
wjelem[i + 1] =
utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp);
// construct cutsq
double cut;
cutmax = 0.0;
memory->create(cutsq,ntypes+1,ntypes+1,"snad/atom:cutsq");
memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq");
for (int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
cut = 2.0 * radelem[i] * rcutfac;
if (cut > cutmax) cutmax = cut;
cutsq[i][i] = cut*cut;
for (int j = i+1; j <= ntypes; j++) {
cut = (radelem[i]+radelem[j])*rcutfac;
cutsq[i][j] = cutsq[j][i] = cut*cut;
cutsq[i][i] = cut * cut;
for (int j = i + 1; j <= ntypes; j++) {
cut = (radelem[i] + radelem[j]) * rcutfac;
cutsq[i][j] = cutsq[j][i] = cut * cut;
}
}
@ -92,93 +98,89 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
rmin0 = atof(arg[iarg+1]);
if (strcmp(arg[iarg], "rmin0") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"bzeroflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
bzeroflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"switchflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
switchflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "bzeroflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"quadraticflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
quadraticflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "quadraticflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2+ntypes > narg)
error->all(FLERR,"Illegal compute snad/atom command");
} else if (strcmp(arg[iarg], "chem") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
chemflag = 1;
memory->create(map,ntypes+1,"compute_snad_atom:map");
nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
memory->create(map, ntypes + 1, "compute_sna_grid:map");
nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
for (int i = 0; i < ntypes; i++) {
int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp);
if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute snad/atom command");
map[i+1] = jelem;
int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp);
if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style);
map[i + 1] = jelem;
}
iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2 + ntypes;
} else if (strcmp(arg[iarg], "bnormflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
wselfallflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "wselfallflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
switchinnerflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchinnerflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"sinner") == 0) {
} else if (strcmp(arg[iarg], "sinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute snad/atom command");
memory->create(sinnerelem,ntypes+1,"snad/atom:sinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style);
memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem");
for (int i = 0; i < ntypes; i++)
sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
sinnerflag = 1;
iarg += ntypes;
} else if (strcmp(arg[iarg],"dinner") == 0) {
} else if (strcmp(arg[iarg], "dinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute snad/atom command");
memory->create(dinnerelem,ntypes+1,"snad/atom:dinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style);
memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem");
for (int i = 0; i < ntypes; i++)
dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
dinnerflag = 1;
iarg += ntypes;
} else error->all(FLERR,"Illegal compute snad/atom command");
} else
error->all(FLERR, "Illegal compute {} command", style);
}
if (switchinnerflag && !(sinnerflag && dinnerflag))
error->all(FLERR,"Illegal compute snad/atom command: switchinnerflag = 1, missing sinner/dinner keyword");
error->all(
FLERR,
"Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword",
style);
if (!switchinnerflag && (sinnerflag || dinnerflag))
error->all(FLERR,"Illegal compute snad/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword");
error->all(
FLERR,
"Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword",
style);
snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag,
chemflag, bnormflag, wselfallflag,
nelements, switchinnerflag);
snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag,
wselfallflag, nelements, switchinnerflag);
ncoeff = snaptr->ncoeff;
nperdim = ncoeff;
if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2;
yoffset = nperdim;
zoffset = 2*nperdim;
size_peratom_cols = 3*nperdim*atom->ntypes;
nvalues = ncoeff;
if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2;
// end code common to all SNAP computes
yoffset = nvalues;
zoffset = 2*nvalues;
size_peratom_cols = 3*nvalues*atom->ntypes;
comm_reverse = size_peratom_cols;
peratom_flag = 1;
@ -289,7 +291,7 @@ void ComputeSNADAtom::compute_peratom()
// const int typeoffset = threencoeff*(atom->type[i]-1);
// const int quadraticoffset = threencoeff*atom->ntypes +
// threencoeffq*(atom->type[i]-1);
const int typeoffset = 3*nperdim*(atom->type[i]-1);
const int typeoffset = 3*nvalues*(atom->type[i]-1);
// insure rij, inside, and typej are of size jnum

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

@ -37,7 +37,7 @@ class ComputeSNADAtom : public Compute {
private:
int nmax;
int ncoeff, nperdim, yoffset, zoffset;
int ncoeff, nvalues, yoffset, zoffset;
double **cutsq;
class NeighList *list;
double **snad;

177
src/ML-SNAP/compute_snap.cpp
View File

@ -41,13 +41,15 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
array_flag = 1;
extarray = 0;
// begin code common to all SNAP computes
double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
int ntypes = atom->ntypes;
int nargmin = 6+2*ntypes;
int nargmin = 6 + 2 * ntypes;
if (narg < nargmin) error->all(FLERR,"Illegal compute snap command");
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
// default values
@ -55,7 +57,6 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
switchflag = 1;
bzeroflag = 1;
quadraticflag = 0;
bikflag = 0;
chemflag = 0;
bnormflag = 0;
wselfallflag = 0;
@ -64,28 +65,32 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
// process required arguments
memory->create(radelem,ntypes+1,"snap:radelem"); // offset by 1 to match up with types
memory->create(wjelem,ntypes+1,"snap:wjelem");
rcutfac = atof(arg[3]);
rfac0 = atof(arg[4]);
twojmax = atoi(arg[5]);
memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types
memory->create(wjelem, ntypes + 1, "sna/atom:wjelem");
rcutfac = utils::numeric(FLERR, arg[3], false, lmp);
rfac0 = utils::numeric(FLERR, arg[4], false, lmp);
twojmax = utils::inumeric(FLERR, arg[5], false, lmp);
for (int i = 0; i < ntypes; i++)
radelem[i+1] = atof(arg[6+i]);
radelem[i + 1] =
utils::numeric(FLERR, arg[6 + i], false, lmp);
for (int i = 0; i < ntypes; i++)
wjelem[i+1] = atof(arg[6+ntypes+i]);
wjelem[i + 1] =
utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp);
// construct cutsq
double cut;
cutmax = 0.0;
memory->create(cutsq,ntypes+1,ntypes+1,"snap:cutsq");
memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq");
for (int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
cut = 2.0 * radelem[i] * rcutfac;
if (cut > cutmax) cutmax = cut;
cutsq[i][i] = cut*cut;
for (int j = i+1; j <= ntypes; j++) {
cut = (radelem[i]+radelem[j])*rcutfac;
cutsq[i][j] = cutsq[j][i] = cut*cut;
cutsq[i][i] = cut * cut;
for (int j = i + 1; j <= ntypes; j++) {
cut = (radelem[i] + radelem[j]) * rcutfac;
cutsq[i][j] = cutsq[j][i] = cut * cut;
}
}
@ -99,107 +104,99 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
rmin0 = atof(arg[iarg+1]);
if (strcmp(arg[iarg], "rmin0") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"bzeroflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
bzeroflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"switchflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
switchflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "bzeroflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"quadraticflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
quadraticflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "quadraticflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2+ntypes > narg)
error->all(FLERR,"Illegal compute snap command");
} else if (strcmp(arg[iarg], "chem") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
chemflag = 1;
memory->create(map,ntypes+1,"compute_snap:map");
nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
memory->create(map, ntypes + 1, "compute_sna_grid:map");
nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
for (int i = 0; i < ntypes; i++) {
int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp);
if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute snap command");
map[i+1] = jelem;
int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp);
if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style);
map[i + 1] = jelem;
}
iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2 + ntypes;
} else if (strcmp(arg[iarg], "bnormflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
wselfallflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "wselfallflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"bikflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
bikflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchinnerflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
switchinnerflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"sinner") == 0) {
} else if (strcmp(arg[iarg], "sinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute snap command");
memory->create(sinnerelem,ntypes+1,"snap:sinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style);
memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem");
for (int i = 0; i < ntypes; i++)
sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
sinnerflag = 1;
iarg += ntypes;
} else if (strcmp(arg[iarg],"dinner") == 0) {
} else if (strcmp(arg[iarg], "dinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute snap command");
memory->create(dinnerelem,ntypes+1,"snap:dinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style);
memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem");
for (int i = 0; i < ntypes; i++)
dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
dinnerflag = 1;
iarg += ntypes;
} else error->all(FLERR,"Illegal compute snap command");
} else
error->all(FLERR, "Illegal compute {} command", style);
}
if (switchinnerflag && !(sinnerflag && dinnerflag))
error->all(FLERR,"Illegal compute snap command: switchinnerflag = 1, missing sinner/dinner keyword");
error->all(
FLERR,
"Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword",
style);
if (!switchinnerflag && (sinnerflag || dinnerflag))
error->all(FLERR,"Illegal compute snap command: switchinnerflag = 0, unexpected sinner/dinner keyword");
error->all(
FLERR,
"Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword",
style);
snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag,
chemflag, bnormflag, wselfallflag,
nelements, switchinnerflag);
snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag,
wselfallflag, nelements, switchinnerflag);
ncoeff = snaptr->ncoeff;
nperdim = ncoeff;
if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2;
nvalues = ncoeff;
if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2;
// end code common to all SNAP computes
ndims_force = 3;
ndims_virial = 6;
yoffset = nperdim;
zoffset = 2*nperdim;
yoffset = nvalues;
zoffset = 2*nvalues;
natoms = atom->natoms;
bik_rows = 1;
if (bikflag) bik_rows = natoms;
size_array_rows = bik_rows+ndims_force*natoms+ndims_virial;
size_array_cols = nperdim*atom->ntypes+1;
size_array_cols = nvalues*atom->ntypes+1;
lastcol = size_array_cols-1;
ndims_peratom = ndims_force;
size_peratom = ndims_peratom*nperdim*atom->ntypes;
size_peratom = ndims_peratom*nvalues*atom->ntypes;
nmax = 0;
}
@ -341,8 +338,8 @@ void ComputeSnap::compute_array()
const double radi = radelem[itype];
const int* const jlist = firstneigh[i];
const int jnum = numneigh[i];
const int typeoffset_local = ndims_peratom*nperdim*(itype-1);
const int typeoffset_global = nperdim*(itype-1);
const int typeoffset_local = ndims_peratom*nvalues*(itype-1);
const int typeoffset_global = nvalues*(itype-1);
// insure rij, inside, and typej are of size jnum
@ -481,9 +478,9 @@ void ComputeSnap::compute_array()
// accumulate bispectrum force contributions to global array
for (int itype = 0; itype < atom->ntypes; itype++) {
const int typeoffset_local = ndims_peratom*nperdim*itype;
const int typeoffset_global = nperdim*itype;
for (int icoeff = 0; icoeff < nperdim; icoeff++) {
const int typeoffset_local = ndims_peratom*nvalues*itype;
const int typeoffset_global = nvalues*itype;
for (int icoeff = 0; icoeff < nvalues; icoeff++) {
for (int i = 0; i < ntotal; i++) {
double *snadi = snap_peratom[i]+typeoffset_local;
int iglobal = atom->tag[i];
@ -549,10 +546,10 @@ void ComputeSnap::dbdotr_compute()
int nall = atom->nlocal + atom->nghost;
for (int i = 0; i < nall; i++)
for (int itype = 0; itype < atom->ntypes; itype++) {
const int typeoffset_local = ndims_peratom*nperdim*itype;
const int typeoffset_global = nperdim*itype;
const int typeoffset_local = ndims_peratom*nvalues*itype;
const int typeoffset_global = nvalues*itype;
double *snadi = snap_peratom[i]+typeoffset_local;
for (int icoeff = 0; icoeff < nperdim; icoeff++) {
for (int icoeff = 0; icoeff < nvalues; icoeff++) {
double dbdx = snadi[icoeff];
double dbdy = snadi[icoeff+yoffset];
double dbdz = snadi[icoeff+zoffset];

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

@ -35,7 +35,7 @@ class ComputeSnap : public Compute {
private:
int natoms, nmax, size_peratom, lastcol;
int ncoeff, nperdim, yoffset, zoffset;
int ncoeff, nvalues, yoffset, zoffset;
int ndims_peratom, ndims_force, ndims_virial;
double **cutsq;
class NeighList *list;

225
src/ML-SNAP/compute_snav_atom.cpp
View File

@ -21,6 +21,7 @@
#include "neighbor.h"
#include "neigh_list.h"
#include "force.h"
#include "pair.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
@ -33,20 +34,22 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snav(nullptr),
radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr)
{
// begin code common to all SNAP computes
double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
int ntypes = atom->ntypes;
int nargmin = 6+2*ntypes;
int nargmin = 6 + 2 * ntypes;
if (narg < nargmin) error->all(FLERR,"Illegal compute snav/atom command");
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
// default values
rmin0 = 0.0;
switchflag = 1;
bzeroflag = 1;
bnormflag = 0;
quadraticflag = 0;
chemflag = 0;
bnormflag = 0;
@ -56,24 +59,32 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
// process required arguments
memory->create(radelem,ntypes+1,"snav/atom:radelem"); // offset by 1 to match up with types
memory->create(wjelem,ntypes+1,"snav/atom:wjelem");
rcutfac = atof(arg[3]);
rfac0 = atof(arg[4]);
twojmax = atoi(arg[5]);
memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types
memory->create(wjelem, ntypes + 1, "sna/atom:wjelem");
rcutfac = utils::numeric(FLERR, arg[3], false, lmp);
rfac0 = utils::numeric(FLERR, arg[4], false, lmp);
twojmax = utils::inumeric(FLERR, arg[5], false, lmp);
for (int i = 0; i < ntypes; i++)
radelem[i+1] = atof(arg[6+i]);
radelem[i + 1] =
utils::numeric(FLERR, arg[6 + i], false, lmp);
for (int i = 0; i < ntypes; i++)
wjelem[i+1] = atof(arg[6+ntypes+i]);
wjelem[i + 1] =
utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp);
// construct cutsq
double cut;
memory->create(cutsq,ntypes+1,ntypes+1,"snav/atom:cutsq");
cutmax = 0.0;
memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq");
for (int i = 1; i <= ntypes; i++) {
cut = 2.0*radelem[i]*rcutfac;
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;
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;
}
}
@ -87,90 +98,87 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
rmin0 = atof(arg[iarg+1]);
if (strcmp(arg[iarg], "rmin0") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"switchflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
switchflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"bzeroflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
bzeroflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "bzeroflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"quadraticflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
quadraticflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "quadraticflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2+ntypes > narg)
error->all(FLERR,"Illegal compute snav/atom command");
} else if (strcmp(arg[iarg], "chem") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
chemflag = 1;
memory->create(map,ntypes+1,"compute_sna_atom:map");
nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
memory->create(map, ntypes + 1, "compute_sna_grid:map");
nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
for (int i = 0; i < ntypes; i++) {
int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp);
if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute snav/atom command");
map[i+1] = jelem;
int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp);
if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style);
map[i + 1] = jelem;
}
iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2 + ntypes;
} else if (strcmp(arg[iarg], "bnormflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
wselfallflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "wselfallflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
switchinnerflag = atoi(arg[iarg+1]);
} else if (strcmp(arg[iarg], "switchinnerflag") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"sinner") == 0) {
} else if (strcmp(arg[iarg], "sinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute snav/atom command");
memory->create(sinnerelem,ntypes+1,"snav/atom:sinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style);
memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem");
for (int i = 0; i < ntypes; i++)
sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
sinnerflag = 1;
iarg += ntypes;
} else if (strcmp(arg[iarg],"dinner") == 0) {
} else if (strcmp(arg[iarg], "dinner") == 0) {
iarg++;
if (iarg+ntypes > narg)
error->all(FLERR,"Illegal compute snav/atom command");
memory->create(dinnerelem,ntypes+1,"snav/atom:dinnerelem");
if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style);
memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem");
for (int i = 0; i < ntypes; i++)
dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp);
dinnerflag = 1;
iarg += ntypes;
} else error->all(FLERR,"Illegal compute snav/atom command");
} else
error->all(FLERR, "Illegal compute {} command", style);
}
if (switchinnerflag && !(sinnerflag && dinnerflag))
error->all(FLERR,"Illegal compute snav/atom command: switchinnerflag = 1, missing sinner/dinner keyword");
error->all(
FLERR,
"Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword",
style);
if (!switchinnerflag && (sinnerflag || dinnerflag))
error->all(FLERR,"Illegal compute snav/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword");
error->all(
FLERR,
"Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword",
style);
snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag,
chemflag, bnormflag, wselfallflag,
nelements, switchinnerflag);
snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag,
wselfallflag, nelements, switchinnerflag);
ncoeff = snaptr->ncoeff;
nperdim = ncoeff;
if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2;
size_peratom_cols = 6*nperdim*atom->ntypes;
nvalues = ncoeff;
if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2;
// end code common to all SNAP computes
size_peratom_cols = 6*nvalues*atom->ntypes;
comm_reverse = size_peratom_cols;
peratom_flag = 1;
@ -203,10 +211,9 @@ void ComputeSNAVAtom::init()
{
if (force->pair == nullptr)
error->all(FLERR,"Compute snav/atom requires a pair style be defined");
// TODO: Not sure what to do with this error check since cutoff radius is not
// a single number
//if (sqrt(cutsq) > force->pair->cutforce)
// error->all(FLERR,"Compute snav/atom cutoff is longer than pairwise cutoff");
if (cutmax > force->pair->cutforce)
error->all(FLERR,"Compute snav/atom cutoff is longer than pairwise cutoff");
// need an occasional full neighbor list
@ -280,7 +287,7 @@ void ComputeSNAVAtom::compute_peratom()
const int* const jlist = firstneigh[i];
const int jnum = numneigh[i];
const int typeoffset = 6*nperdim*(atom->type[i]-1);
const int typeoffset = 6*nvalues*(atom->type[i]-1);
// insure rij, inside, and typej are of size jnum
@ -339,17 +346,17 @@ void ComputeSNAVAtom::compute_peratom()
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
snavi[icoeff] += snaptr->dblist[icoeff][0]*xtmp;
snavi[icoeff+nperdim] += snaptr->dblist[icoeff][1]*ytmp;
snavi[icoeff+2*nperdim] += snaptr->dblist[icoeff][2]*ztmp;
snavi[icoeff+3*nperdim] += snaptr->dblist[icoeff][1]*ztmp;
snavi[icoeff+4*nperdim] += snaptr->dblist[icoeff][0]*ztmp;
snavi[icoeff+5*nperdim] += snaptr->dblist[icoeff][0]*ytmp;
snavi[icoeff+nvalues] += snaptr->dblist[icoeff][1]*ytmp;
snavi[icoeff+2*nvalues] += snaptr->dblist[icoeff][2]*ztmp;
snavi[icoeff+3*nvalues] += snaptr->dblist[icoeff][1]*ztmp;
snavi[icoeff+4*nvalues] += snaptr->dblist[icoeff][0]*ztmp;
snavi[icoeff+5*nvalues] += snaptr->dblist[icoeff][0]*ytmp;
snavj[icoeff] -= snaptr->dblist[icoeff][0]*x[j][0];
snavj[icoeff+nperdim] -= snaptr->dblist[icoeff][1]*x[j][1];
snavj[icoeff+2*nperdim] -= snaptr->dblist[icoeff][2]*x[j][2];
snavj[icoeff+3*nperdim] -= snaptr->dblist[icoeff][1]*x[j][2];
snavj[icoeff+4*nperdim] -= snaptr->dblist[icoeff][0]*x[j][2];
snavj[icoeff+5*nperdim] -= snaptr->dblist[icoeff][0]*x[j][1];
snavj[icoeff+nvalues] -= snaptr->dblist[icoeff][1]*x[j][1];
snavj[icoeff+2*nvalues] -= snaptr->dblist[icoeff][2]*x[j][2];
snavj[icoeff+3*nvalues] -= snaptr->dblist[icoeff][1]*x[j][2];
snavj[icoeff+4*nvalues] -= snaptr->dblist[icoeff][0]*x[j][2];
snavj[icoeff+5*nvalues] -= snaptr->dblist[icoeff][0]*x[j][1];
}
if (quadraticflag) {
@ -369,17 +376,17 @@ void ComputeSNAVAtom::compute_peratom()
double dbytmp = bi*biy;
double dbztmp = bi*biz;
snavi[ncount] += dbxtmp*xtmp;
snavi[ncount+nperdim] += dbytmp*ytmp;
snavi[ncount+2*nperdim] += dbztmp*ztmp;
snavi[ncount+3*nperdim] += dbytmp*ztmp;
snavi[ncount+4*nperdim] += dbxtmp*ztmp;
snavi[ncount+5*nperdim] += dbxtmp*ytmp;
snavi[ncount+nvalues] += dbytmp*ytmp;
snavi[ncount+2*nvalues] += dbztmp*ztmp;
snavi[ncount+3*nvalues] += dbytmp*ztmp;
snavi[ncount+4*nvalues] += dbxtmp*ztmp;
snavi[ncount+5*nvalues] += dbxtmp*ytmp;
snavj[ncount] -= dbxtmp*x[j][0];
snavj[ncount+nperdim] -= dbytmp*x[j][1];
snavj[ncount+2*nperdim] -= dbztmp*x[j][2];
snavj[ncount+3*nperdim] -= dbytmp*x[j][2];
snavj[ncount+4*nperdim] -= dbxtmp*x[j][2];
snavj[ncount+5*nperdim] -= dbxtmp*x[j][1];
snavj[ncount+nvalues] -= dbytmp*x[j][1];
snavj[ncount+2*nvalues] -= dbztmp*x[j][2];
snavj[ncount+3*nvalues] -= dbytmp*x[j][2];
snavj[ncount+4*nvalues] -= dbxtmp*x[j][2];
snavj[ncount+5*nvalues] -= dbxtmp*x[j][1];
ncount++;
// upper-triangular elements of quadratic matrix
@ -392,17 +399,17 @@ void ComputeSNAVAtom::compute_peratom()
double dbztmp = bi*snaptr->dblist[jcoeff][2]
+ biz*snaptr->blist[jcoeff];
snavi[ncount] += dbxtmp*xtmp;
snavi[ncount+nperdim] += dbytmp*ytmp;
snavi[ncount+2*nperdim] += dbztmp*ztmp;
snavi[ncount+3*nperdim] += dbytmp*ztmp;
snavi[ncount+4*nperdim] += dbxtmp*ztmp;
snavi[ncount+5*nperdim] += dbxtmp*ytmp;
snavi[ncount+nvalues] += dbytmp*ytmp;
snavi[ncount+2*nvalues] += dbztmp*ztmp;
snavi[ncount+3*nvalues] += dbytmp*ztmp;
snavi[ncount+4*nvalues] += dbxtmp*ztmp;
snavi[ncount+5*nvalues] += dbxtmp*ytmp;
snavj[ncount] -= dbxtmp*x[j][0];
snavj[ncount+nperdim] -= dbytmp*x[j][1];
snavj[ncount+2*nperdim] -= dbztmp*x[j][2];
snavj[ncount+3*nperdim] -= dbytmp*x[j][2];
snavj[ncount+4*nperdim] -= dbxtmp*x[j][2];
snavj[ncount+5*nperdim] -= dbxtmp*x[j][1];
snavj[ncount+nvalues] -= dbytmp*x[j][1];
snavj[ncount+2*nvalues] -= dbztmp*x[j][2];
snavj[ncount+3*nvalues] -= dbytmp*x[j][2];
snavj[ncount+4*nvalues] -= dbxtmp*x[j][2];
snavj[ncount+5*nvalues] -= dbxtmp*x[j][1];
ncount++;
}
}

3
src/ML-SNAP/compute_snav_atom.h
View File

@ -37,7 +37,7 @@ class ComputeSNAVAtom : public Compute {
private:
int nmax;
int ncoeff, nperdim;
int ncoeff, nvalues;
double **cutsq;
class NeighList *list;
double **snav;
@ -50,6 +50,7 @@ class ComputeSNAVAtom : public Compute {
double *sinnerelem;
double *dinnerelem;
class SNA *snaptr;
double cutmax;
int quadraticflag;
};