diff --git a/doc/src/compute_sna_atom.txt b/doc/src/compute_sna_atom.txt index efbf2e9ea3..10e68f5698 100644 --- a/doc/src/compute_sna_atom.txt +++ b/doc/src/compute_sna_atom.txt @@ -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 diff --git a/doc/src/pair_snap.txt b/doc/src/pair_snap.txt index a796cfdeba..1fba74a188 100644 --- a/doc/src/pair_snap.txt +++ b/doc/src/pair_snap.txt @@ -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]: diff --git a/potentials/Ta06A.snapparam b/potentials/Ta06A.snapparam index 283629d658..629d96d708 100644 --- a/potentials/Ta06A.snapparam +++ b/potentials/Ta06A.snapparam @@ -10,6 +10,5 @@ twojmax 6 rfac0 0.99363 rmin0 0 -diagonalstyle 3 bzeroflag 0 quadraticflag 0 diff --git a/potentials/W_2940_2017_2.snapparam b/potentials/W_2940_2017_2.snapparam index 27ab61a266..49f3094d08 100644 --- a/potentials/W_2940_2017_2.snapparam +++ b/potentials/W_2940_2017_2.snapparam @@ -8,6 +8,5 @@ twojmax 8 rfac0 0.99363 rmin0 0 -diagonalstyle 3 bzeroflag 0 quadraticflag 0 diff --git a/src/SNAP/compute_sna_atom.cpp b/src/SNAP/compute_sna_atom.cpp index fea37faca0..cc7a84281e 100644 --- a/src/SNAP/compute_sna_atom.cpp +++ b/src/SNAP/compute_sna_atom.cpp @@ -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; } diff --git a/src/SNAP/compute_sna_atom.h b/src/SNAP/compute_sna_atom.h index 56ffccfa7e..105a62a37a 100644 --- a/src/SNAP/compute_sna_atom.h +++ b/src/SNAP/compute_sna_atom.h @@ -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; diff --git a/src/SNAP/compute_snad_atom.cpp b/src/SNAP/compute_snad_atom.cpp index 156380eccc..37587a0aae 100644 --- a/src/SNAP/compute_snad_atom.cpp +++ b/src/SNAP/compute_snad_atom.cpp @@ -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; } diff --git a/src/SNAP/compute_snad_atom.h b/src/SNAP/compute_snad_atom.h index 1fcf540d7c..ac353d8553 100644 --- a/src/SNAP/compute_snad_atom.h +++ b/src/SNAP/compute_snad_atom.h @@ -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; diff --git a/src/SNAP/compute_snav_atom.cpp b/src/SNAP/compute_snav_atom.cpp index 6caff0820c..1f702496ed 100644 --- a/src/SNAP/compute_snav_atom.cpp +++ b/src/SNAP/compute_snav_atom.cpp @@ -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; } diff --git a/src/SNAP/compute_snav_atom.h b/src/SNAP/compute_snav_atom.h index 6bcce346e0..9df17cc667 100644 --- a/src/SNAP/compute_snav_atom.h +++ b/src/SNAP/compute_snav_atom.h @@ -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; diff --git a/src/SNAP/pair_snap.cpp b/src/SNAP/pair_snap.cpp index 6a65f872fd..6eb05f85a4 100644 --- a/src/SNAP/pair_snap.cpp +++ b/src/SNAP/pair_snap.cpp @@ -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; } diff --git a/src/SNAP/pair_snap.h b/src/SNAP/pair_snap.h index b5871c1527..c64eaa5d4e 100644 --- a/src/SNAP/pair_snap.h +++ b/src/SNAP/pair_snap.h @@ -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 diff --git a/src/SNAP/sna.cpp b/src/SNAP/sna.cpp index 131ac48fdb..75601b8e17 100644 --- a/src/SNAP/sna.cpp +++ b/src/SNAP/sna.cpp @@ -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++; diff --git a/src/SNAP/sna.h b/src/SNAP/sna.h index b54ad3482a..1e08ef123c 100644 --- a/src/SNAP/sna.h +++ b/src/SNAP/sna.h @@ -18,9 +18,7 @@ #ifndef LMP_SNA_H #define LMP_SNA_H -#include #include "pointers.h" -#include 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 }; }