jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Completed back-porting of Zombie SNAP improvements, particularly noteworthy is reduction in memory footprint, elimination of most multidimensional arrays, elimination of diagonal_style, elimination of Z array in force calculation.

This commit is contained in:
Aidan Thompson 2019-06-12 16:42:28 -06:00
parent 0559e155f2
commit a973700295
14 changed files with 130 additions and 154 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
doc/src/compute_sna_atom.txt
View File

@ -24,12 +24,8 @@ twojmax = band limit for bispectrum components (non-negative integer) :l
R_1, R_2,... = list of cutoff radii, one for each type (distance units) :l
w_1, w_2,... = list of neighbor weights, one for each type :l
zero or more keyword/value pairs may be appended :l
keyword = {diagonal} or {rmin0} or {switchflag} or {bzeroflag} or {quadraticflag} :l
{diagonal} value = {0} or {1} or {2} or {3}
{0} = all j1, j2, j <= twojmax, j2 <= j1
{1} = subset satisfying j1 == j2
{2} = subset satisfying j1 == j2 == j3
{3} = subset satisfying j2 <= j1 <= j
keyword = {rmin0} or {switchflag} or {bzeroflag} or {quadraticflag} :l
// {3} = subset satisfying j2 <= j1 <= j
{rmin0} value = parameter in distance to angle conversion (distance units)
{switchflag} value = {0} or {1}
{0} = do not use switching function
@ -44,7 +40,7 @@ keyword = {diagonal} or {rmin0} or {switchflag} or {bzeroflag} or {quadraticflag
[Examples:]
compute b all sna/atom 1.4 0.99363 6 2.0 2.4 0.75 1.0 diagonal 3 rmin0 0.0
compute b all sna/atom 1.4 0.99363 6 2.0 2.4 0.75 1.0 rmin0 0.0
compute db all sna/atom 1.4 0.95 6 2.0 1.0
compute vb all sna/atom 1.4 0.95 6 2.0 1.0 :pre
@ -151,7 +147,7 @@ The argument {rfac0} and the optional keyword {rmin0} define the
linear mapping from radial distance to polar angle {theta0} on the
3-sphere.
The argument {twojmax} and the keyword {diagonal} define which
The argument {twojmax} defines which
bispectrum components are generated. See section below on output for a
detailed explanation of the number of bispectrum components and the
ordered in which they are listed.
@ -192,25 +188,20 @@ command that includes all pairs in the neighbor list.
Compute {sna/atom} calculates a per-atom array, each column
corresponding to a particular bispectrum component. The total number
of columns and the identity of the bispectrum component contained in
each column depend on the values of {twojmax} and {diagonal}, as
each column depend of the value of {twojmax}, as
described by the following piece of python code:
for j1 in range(0,twojmax+1):
if(diagonal==2):
print j1/2.,j1/2.,j1/2.
elif(diagonal==1):
for j in range(0,min(twojmax,2*j1)+1,2):
print j1/2.,j1/2.,j/2.
elif(diagonal==0):
for j2 in range(0,j1+1):
for j in range(j1-j2,min(twojmax,j1+j2)+1,2):
print j1/2.,j2/2.,j/2.
elif(diagonal==3):
for j2 in range(0,j1+1):
for j in range(j1-j2,min(twojmax,j1+j2)+1,2):
if (j>=j1): print j1/2.,j2/2.,j/2. :pre
for j2 in range(0,j1+1):
for j in range(j1-j2,min(twojmax,j1+j2)+1,2):
if (j>=j1): print j1/2.,j2/2.,j/2. :pre
Compute {snad/atom} evaluates a per-atom array. The columns are
NOTE: the {diagonal} keyword allowing other possible choices
for the number of bispectrum components was removed in 2019,
since all potentials use the value of 3, corresponding to the
above set of bispectrum components.
ompute {snad/atom} evaluates a per-atom array. The columns are
arranged into {ntypes} blocks, listed in order of atom type {I}. Each
block contains three sub-blocks corresponding to the {x}, {y}, and {z}
components of the atom position. Each of these sub-blocks contains
@ -259,7 +250,7 @@ package"_Build_package.html doc page for more info.
[Default:]
The optional keyword defaults are {diagonal} = 0, {rmin0} = 0,
The optional keyword defaults are {rmin0} = 0,
{switchflag} = 1, {bzeroflag} = 1, {quadraticflag} = 0,
:line

8
doc/src/pair_snap.txt
View File

@ -38,7 +38,7 @@ where {B_k^i} is the {k}-th bispectrum component of atom {i},
and {beta_k^alpha_i} is the corresponding linear coefficient
that depends on {alpha_i}, the SNAP element of atom {i}. The
number of bispectrum components used and their definitions
depend on the values of {twojmax} and {diagonalstyle}
depend on the value of {twojmax}
defined in the SNAP parameter file described below.
The bispectrum calculation is described in more detail
in "compute sna/atom"_compute_sna_atom.html.
@ -125,14 +125,13 @@ This line is followed by {ncoeff} coefficients, one per line.
The SNAP parameter file can contain blank and comment lines (start
with #) anywhere. Each non-blank non-comment line must contain one
keyword/value pair. The required keywords are {rcutfac} and
{twojmax}. Optional keywords are {rfac0}, {rmin0}, {diagonalstyle},
{twojmax}. Optional keywords are {rfac0}, {rmin0},
{switchflag}, and {bzeroflag}.
The default values for these keywords are
{rfac0} = 0.99363
{rmin0} = 0.0
{diagonalstyle} = 3
{switchflag} = 0
{bzeroflag} = 1
{quadraticflag} = 1 :ul
@ -144,6 +143,9 @@ If {quadraticflag} is set to 1, then the SNAP energy expression includes the qua
The SNAP element file should contain {K}({K}+1)/2 additional coefficients
for each element, the upper-triangular elements of alpha.
NOTE: The previously used {diagonalstyle} keyword was removed in 2019,
since all known SNAP potentials use the default value of 3.
:line
[Mixing, shift, table, tail correction, restart, rRESPA info]:

1
potentials/Ta06A.snapparam
View File

@ -10,6 +10,5 @@ twojmax 6
rfac0 0.99363
rmin0 0
diagonalstyle 3
bzeroflag 0
quadraticflag 0

1
potentials/W_2940_2017_2.snapparam
View File

@ -8,6 +8,5 @@ twojmax 8
rfac0 0.99363
rmin0 0
diagonalstyle 3
bzeroflag 0
quadraticflag 0

19
src/SNAP/compute_sna_atom.cpp
View File

@ -44,7 +44,6 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
// default values
diagonalstyle = 0;
rmin0 = 0.0;
switchflag = 1;
bzeroflag = 1;
@ -84,14 +83,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"diagonal") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
diagonalstyle = atoi(arg[iarg+1]);
if (diagonalstyle < 0 || diagonalstyle > 3)
error->all(FLERR,"Illegal compute sna/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"rmin0") == 0) {
if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
rmin0 = atof(arg[iarg+1]);
@ -114,7 +106,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
} else error->all(FLERR,"Illegal compute sna/atom command");
}
snaptr = new SNA(lmp,rfac0,twojmax,diagonalstyle,
snaptr = new SNA(lmp,rfac0,twojmax,
rmin0,switchflag,bzeroflag);
ncoeff = snaptr->ncoeff;
@ -123,7 +115,6 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
peratom_flag = 1;
nmax = 0;
njmax = 0;
sna = NULL;
}
@ -277,9 +268,9 @@ void ComputeSNAAtom::compute_peratom()
double ComputeSNAAtom::memory_usage()
{
double bytes = nmax*size_peratom_cols * sizeof(double);
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
double bytes = nmax*size_peratom_cols * sizeof(double); // sna
bytes += snaptr->memory_usage(); // SNA object
return bytes;
}

2
src/SNAP/compute_sna_atom.h
View File

@ -34,7 +34,7 @@ class ComputeSNAAtom : public Compute {
double memory_usage();
private:
int nmax, njmax, diagonalstyle;
int nmax;
int ncoeff;
double **cutsq;
class NeighList *list;

21
src/SNAP/compute_snad_atom.cpp
View File

@ -44,7 +44,6 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
// default values
diagonalstyle = 0;
rmin0 = 0.0;
switchflag = 1;
bzeroflag = 1;
@ -82,14 +81,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"diagonal") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
diagonalstyle = atof(arg[iarg+1]);
if (diagonalstyle < 0 || diagonalstyle > 3)
error->all(FLERR,"Illegal compute snad/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"rmin0") == 0) {
if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
rmin0 = atof(arg[iarg+1]);
@ -112,7 +104,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
} else error->all(FLERR,"Illegal compute snad/atom command");
}
snaptr = new SNA(lmp,rfac0,twojmax,diagonalstyle,
snaptr = new SNA(lmp,rfac0,twojmax,
rmin0,switchflag,bzeroflag);
ncoeff = snaptr->ncoeff;
@ -125,7 +117,6 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
peratom_flag = 1;
nmax = 0;
njmax = 0;
snad = NULL;
}
@ -378,9 +369,9 @@ void ComputeSNADAtom::unpack_reverse_comm(int n, int *list, double *buf)
double ComputeSNADAtom::memory_usage()
{
double bytes = nmax*size_peratom_cols * sizeof(double);
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
bytes += 3*nperdim*atom->ntypes;
double bytes = nmax*size_peratom_cols * sizeof(double); // snad
bytes += snaptr->memory_usage(); // SNA object
return bytes;
}

2
src/SNAP/compute_snad_atom.h
View File

@ -36,7 +36,7 @@ class ComputeSNADAtom : public Compute {
double memory_usage();
private:
int nmax, njmax, diagonalstyle;
int nmax;
int ncoeff, nperdim, yoffset, zoffset;
double **cutsq;
class NeighList *list;

21
src/SNAP/compute_snav_atom.cpp
View File

@ -44,7 +44,6 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
// default values
diagonalstyle = 0;
rmin0 = 0.0;
switchflag = 1;
bzeroflag = 1;
@ -78,14 +77,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = nargmin;
while (iarg < narg) {
if (strcmp(arg[iarg],"diagonal") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
diagonalstyle = atof(arg[iarg+1]);
if (diagonalstyle < 0 || diagonalstyle > 3)
error->all(FLERR,"Illegal compute snav/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"rmin0") == 0) {
if (strcmp(arg[iarg],"rmin0") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
rmin0 = atof(arg[iarg+1]);
@ -108,7 +100,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
} else error->all(FLERR,"Illegal compute snav/atom command");
}
snaptr = new SNA(lmp,rfac0,twojmax,diagonalstyle,
snaptr = new SNA(lmp,rfac0,twojmax,
rmin0,switchflag,bzeroflag);
ncoeff = snaptr->ncoeff;
@ -119,7 +111,6 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
peratom_flag = 1;
nmax = 0;
njmax = 0;
snav = NULL;
}
@ -389,10 +380,8 @@ void ComputeSNAVAtom::unpack_reverse_comm(int n, int *list, double *buf)
double ComputeSNAVAtom::memory_usage()
{
double bytes = nmax*size_peratom_cols * sizeof(double);
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
bytes += 6*nperdim*atom->ntypes;
if (quadraticflag) bytes += 6*nperdim*atom->ntypes;
double bytes = nmax*size_peratom_cols * sizeof(double); // snav
bytes += snaptr->memory_usage(); // SNA object
return bytes;
}

2
src/SNAP/compute_snav_atom.h
View File

@ -36,7 +36,7 @@ class ComputeSNAVAtom : public Compute {
double memory_usage();
private:
int nmax, njmax, diagonalstyle;
int nmax;
int ncoeff, nperdim;
double **cutsq;
class NeighList *list;

75
src/SNAP/pair_snap.cpp
View File

@ -34,10 +34,6 @@ using namespace LAMMPS_NS;
#define MAXLINE 1024
#define MAXWORD 3
// Outstanding issues with quadratic term
// 1. there seems to a problem with compute_optimized energy calc
// it does not match compute_regular, even when quadratic coeffs = 0
/* ---------------------------------------------------------------------- */
PairSNAP::PairSNAP(LAMMPS *lmp) : Pair(lmp)
@ -53,8 +49,6 @@ PairSNAP::PairSNAP(LAMMPS *lmp) : Pair(lmp)
wjelem = NULL;
coeffelem = NULL;
nmax = 0;
beta_max = 0;
beta = NULL;
bispectrum = NULL;
@ -74,6 +68,7 @@ PairSNAP::~PairSNAP()
memory->destroy(wjelem);
memory->destroy(coeffelem);
}
memory->destroy(beta);
memory->destroy(bispectrum);
@ -115,7 +110,8 @@ void PairSNAP::compute(int eflag, int vflag)
// compute dE_i/dB_i = beta_i for all i in list
compute_bispectrum();
if (quadraticflag || eflag)
compute_bispectrum();
compute_beta();
numneigh = list->numneigh;
@ -209,15 +205,25 @@ void PairSNAP::compute(int eflag, int vflag)
evdwl = coeffi[0];
// E = beta.B + 0.5*B^t.alpha.B
// coeff[k] = beta[k-1] or
// coeff[k] = alpha_ii or
// coeff[k] = alpha_ij = alpha_ji, j != i
// linear contributions
for (int k = 0; k < ncoeff; k++)
evdwl += beta[ii][k]*bispectrum[ii][k];
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
evdwl += coeffi[icoeff+1]*bispectrum[ii][icoeff];
// quadratic contributions
if (quadraticflag) {
int k = ncoeff+1;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bveci = bispectrum[ii][icoeff];
evdwl += 0.5*coeffi[k++]*bveci*bveci;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
double bvecj = bispectrum[ii][jcoeff];
evdwl += coeffi[k++]*bveci*bvecj;
}
}
}
ev_tally_full(i,2.0*evdwl,0.0,0.0,0.0,0.0,0.0);
}
@ -241,8 +247,23 @@ void PairSNAP::compute_beta()
const int ielem = map[itype];
double* coeffi = coeffelem[ielem];
for (int k = 1; k <= ncoeff; k++)
beta[ii][k-1] = coeffi[k];
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
beta[ii][icoeff] = coeffi[icoeff+1];
if (quadraticflag) {
int k = ncoeff+1;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bveci = bispectrum[ii][icoeff];
beta[ii][icoeff] += coeffi[k]*bveci;
k++;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
double bvecj = bispectrum[ii][jcoeff];
beta[ii][icoeff] += coeffi[k]*bvecj;
beta[ii][jcoeff] += coeffi[k]*bveci;
k++;
}
}
}
}
}
@ -308,8 +329,8 @@ void PairSNAP::compute_bispectrum()
snaptr->compute_zi();
snaptr->compute_bi();
for (int k = 0; k < ncoeff; k++)
bispectrum[ii][k] = snaptr->blist[k];
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
bispectrum[ii][icoeff] = snaptr->blist[icoeff];
}
}
@ -354,8 +375,6 @@ void PairSNAP::coeff(int narg, char **arg)
memory->destroy(wjelem);
memory->destroy(coeffelem);
}
memory->destroy(beta);
memory->destroy(bispectrum);
char* type1 = arg[0];
char* type2 = arg[1];
@ -425,9 +444,7 @@ void PairSNAP::coeff(int narg, char **arg)
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
snaptr = new SNA(lmp,rfac0,twojmax,
diagonalstyle,
rmin0,switchflag,bzeroflag);
snaptr->grow_rij(nmax);
if (ncoeff != snaptr->ncoeff) {
if (comm->me == 0)
@ -617,7 +634,6 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
rfac0 = 0.99363;
rmin0 = 0.0;
diagonalstyle = 3;
switchflag = 1;
bzeroflag = 1;
quadraticflag = 0;
@ -678,8 +694,6 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
rfac0 = atof(keyval);
else if (strcmp(keywd,"rmin0") == 0)
rmin0 = atof(keyval);
else if (strcmp(keywd,"diagonalstyle") == 0)
diagonalstyle = atoi(keyval);
else if (strcmp(keywd,"switchflag") == 0)
switchflag = atoi(keyval);
else if (strcmp(keywd,"bzeroflag") == 0)
@ -702,13 +716,16 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
double PairSNAP::memory_usage()
{
double bytes = Pair::memory_usage();
int n = atom->ntypes+1;
bytes += n*n*sizeof(int);
bytes += n*n*sizeof(double);
bytes += 3*nmax*sizeof(double);
bytes += nmax*sizeof(int);
bytes += (2*ncoeffall)*sizeof(double);
bytes += (ncoeff*3)*sizeof(double);
bytes += n*n*sizeof(int); // setflag
bytes += n*n*sizeof(double); // cutsq
bytes += n*sizeof(int); // map
bytes += beta_max*ncoeff*sizeof(double); // bispectrum
bytes += beta_max*ncoeff*sizeof(double); // beta
bytes += snaptr->memory_usage(); // SNA object
return bytes;
}

4
src/SNAP/pair_snap.h
View File

@ -40,9 +40,7 @@ public:
protected:
int ncoeffq, ncoeffall;
double **bvec, ***dbvec;
class SNA* snaptr;
int nmax;
virtual void allocate();
void read_files(char *, char *);
inline int equal(double* x,double* y);
@ -60,7 +58,7 @@ protected:
double** beta; // betas for all atoms in list
double** bispectrum; // bispectrum components for all atoms in list
int *map; // mapping from atom types to elements
int twojmax, diagonalstyle, switchflag, bzeroflag;
int twojmax, switchflag, bzeroflag;
double rfac0, rmin0, wj1, wj2;
int rcutfacflag, twojmaxflag; // flags for required parameters
int beta_max; // length of beta

80
src/SNAP/sna.cpp
View File

@ -113,7 +113,7 @@ using namespace MathConst;
------------------------------------------------------------------------- */
SNA::SNA(LAMMPS* lmp, double rfac0_in,
int twojmax_in, int diagonalstyle_in,
int twojmax_in,
double rmin0_in, int switch_flag_in, int bzero_flag_in) : Pointers(lmp)
{
wself = 1.0;
@ -124,21 +124,16 @@ SNA::SNA(LAMMPS* lmp, double rfac0_in,
bzero_flag = bzero_flag_in;
twojmax = twojmax_in;
diagonalstyle = diagonalstyle_in;
ncoeff = compute_ncoeff();
bvec = NULL;
dbvec = NULL;
memory->create(bvec, ncoeff, "pair:bvec");
memory->create(dbvec, ncoeff, 3, "pair:dbvec");
rij = NULL;
inside = NULL;
wj = NULL;
rcutij = NULL;
nmax = 0;
idxz = NULL;
idxb= NULL;
idxb = NULL;
build_indexlist();
create_twojmax_arrays();
@ -159,8 +154,6 @@ SNA::~SNA()
memory->destroy(inside);
memory->destroy(wj);
memory->destroy(rcutij);
memory->destroy(bvec);
memory->destroy(dbvec);
delete[] idxz;
delete[] idxb;
destroy_twojmax_arrays();
@ -168,8 +161,6 @@ SNA::~SNA()
void SNA::build_indexlist()
{
if(diagonalstyle != 3)
error->all(FLERR, "diagonal_style must be 3\n");
// index list for cglist
@ -180,7 +171,7 @@ void SNA::build_indexlist()
int idxcg_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2) {
idxcg_block[j1][j2][j] = idxcg_count;
for (int m1 = 0; m1 <= j1; m1++)
for (int m2 = 0; m2 <= j2; m2++)
@ -209,7 +200,7 @@ void SNA::build_indexlist()
int idxb_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) idxb_count++;
idxb_max = idxb_count;
@ -218,7 +209,7 @@ void SNA::build_indexlist()
idxb_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) {
idxb[idxb_count].j1 = j1;
idxb[idxb_count].j2 = j2;
@ -233,7 +224,7 @@ void SNA::build_indexlist()
idxb_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2) {
if (j >= j1) {
idxb_block[j1][j2][j] = idxb_count;
idxb_count++;
@ -246,7 +237,7 @@ void SNA::build_indexlist()
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2)
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2)
for (int mb = 0; 2*mb <= j; mb++)
for (int ma = 0; ma <= j; ma++)
idxz_count++;
@ -260,7 +251,7 @@ void SNA::build_indexlist()
idxz_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2) {
idxz_block[j1][j2][j] = idxz_count;
// find right beta[jjb] entry
@ -288,6 +279,7 @@ void SNA::build_indexlist()
}
}
}
/* ---------------------------------------------------------------------- */
void SNA::init()
@ -312,6 +304,7 @@ void SNA::grow_rij(int newnmax)
memory->create(wj, nmax, "pair:wj");
memory->create(rcutij, nmax, "pair:rcutij");
}
/* ----------------------------------------------------------------------
compute Ui by summing over neighbors j
------------------------------------------------------------------------- */
@ -401,12 +394,6 @@ void SNA::compute_zi()
icgb += j2;
} // end loop over ib
// // apply symmetry factor
// const double jfac = 1.0/(j+1);
// zlist_r[jjz] *= jfac;
// zlist_i[jjz] *= jfac;
} // end loop over jjz
}
@ -631,6 +618,11 @@ void SNA::compute_bi()
} // end if jeven
blist[jjb] = 2.0*sumzu;
// apply bzero shift
if (bzero_flag)
blist[jjb] -= bzero[j];
}
}
@ -1186,25 +1178,38 @@ double SNA::memory_usage()
int jdimpq = twojmax + 2;
int jdim = twojmax + 1;
double bytes;
bytes = ncoeff * sizeof(double); // coeff
bytes = 0;
bytes += jdimpq*jdimpq * sizeof(double); // pqarray
bytes += idxcg_max * sizeof(double); // cglist
bytes += jdim * jdim * jdim * sizeof(int); // idxcg_block
bytes += idxu_max * sizeof(double) * 2; // ulist
bytes += idxu_max * sizeof(double) * 2; // ulisttot
bytes += idxu_max * 3 * sizeof(double) * 2; // dulist
bytes += jdim * sizeof(int); // idxu_block
bytes += idxz_max * 9 * sizeof(int); // idxz
bytes += idxz_max * sizeof(double) * 2; // zlist
bytes += jdim * jdim * jdim * sizeof(int); // idxz_block
bytes += idxb_max * sizeof(double); // blist
bytes += idxb_max * 3 * sizeof(double); // dblist
bytes += idxu_max * sizeof(double) * 2; // ylist
bytes += idxb_max * 3 * sizeof(int); // idxb
bytes += jdim * jdim * jdim * sizeof(int); // idxb_block
bytes += jdim * jdim * jdim * sizeof(int); // idxcg_block
bytes += jdim * sizeof(int); // idxu_block
bytes += jdim * jdim * jdim * sizeof(int); // idxz_block
bytes += jdim * jdim * jdim * sizeof(int); // idxb_block
bytes += idxz_max * sizeof(SNA_ZINDICES); // idxz
bytes += idxb_max * sizeof(SNA_BINDICES); // idxb
bytes += jdim * sizeof(double); // bzero
bytes += nmax * 3 * sizeof(double); // rij
bytes += nmax * sizeof(int); // inside
bytes += nmax * sizeof(double); // wj
bytes += nmax * sizeof(double); // rcutij
printf("SNAP Z list Memory Usage %d\n",idxz_max * sizeof(double) * 2);
printf("SNAP CG list Memory Usage %d\n",idxcg_max * sizeof(double));
return bytes;
}
@ -1242,25 +1247,22 @@ void SNA::destroy_twojmax_arrays()
{
memory->destroy(rootpqarray);
memory->destroy(cglist);
memory->destroy(idxcg_block);
memory->destroy(ulist_r);
memory->destroy(ulist_i);
memory->destroy(ulisttot_r);
memory->destroy(ulisttot_i);
memory->destroy(dulist_r);
memory->destroy(dulist_i);
memory->destroy(idxu_block);
memory->destroy(zlist_r);
memory->destroy(zlist_i);
memory->destroy(blist);
memory->destroy(dblist);
memory->destroy(idxz_block);
memory->destroy(ylist_r);
memory->destroy(ylist_i);
memory->destroy(idxcg_block);
memory->destroy(idxu_block);
memory->destroy(idxz_block);
memory->destroy(idxb_block);
if (bzero_flag)
@ -1484,7 +1486,7 @@ void SNA::init_clebsch_gordan()
int idxcg_count = 0;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
for(int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2) {
for (int m1 = 0; m1 <= j1; m1++) {
aa2 = 2 * m1 - j1;
@ -1557,7 +1559,7 @@ int SNA::compute_ncoeff()
for (int j1 = 0; j1 <= twojmax; j1++)
for (int j2 = 0; j2 <= j1; j2++)
for (int j = abs(j1 - j2);
for (int j = j1 - j2;
j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) ncount++;

9
src/SNAP/sna.h
View File

@ -18,9 +18,7 @@
#ifndef LMP_SNA_H
#define LMP_SNA_H
#include <complex>
#include "pointers.h"
#include <ctime>
namespace LAMMPS_NS {
@ -35,7 +33,7 @@ struct SNA_BINDICES {
class SNA : protected Pointers {
public:
SNA(LAMMPS*, double, int, int, double, int, int);
SNA(LAMMPS*, double, int, double, int, int);
SNA(LAMMPS* lmp) : Pointers(lmp) {};
~SNA();
@ -61,7 +59,6 @@ public:
double compute_sfac(double, double);
double compute_dsfac(double, double);
double* bvec, ** dbvec;
double* blist;
double** dblist;
double** rij;
@ -72,7 +69,7 @@ public:
void grow_rij(int);
int twojmax, diagonalstyle;
int twojmax;
private:
double rmin0, rfac0;
@ -126,7 +123,7 @@ private:
double wself;
int bzero_flag; // 1 if bzero subtracted from barray
double *bzero; // array of B values for isolated atoms
double* bzero; // array of B values for isolated atoms
};
}