diff --git a/doc/src/compute_sna_atom.txt b/doc/src/compute_sna_atom.txt index 1c3787e696..268e23ac28 100644 --- a/doc/src/compute_sna_atom.txt +++ b/doc/src/compute_sna_atom.txt @@ -161,9 +161,9 @@ function. The keyword {bzeroflag} determines whether or not {B0}, the bispectrum components of an atom with no neighbors, are subtracted from -the calculated bispectrum components. This optional keyword is only -available for compute {sna/atom}, as {snad/atom} and {snav/atom} -are unaffected by the removal of constant terms. +the calculated bispectrum components. This optional keyword +normally only affects compute {sna/atom}. However, when +{quadraticflag} is on, it also affects {snad/atom} and {snav/atom}. The keyword {quadraticflag} determines whether or not the quadratic analogs to the bispectrum quantities are generated. @@ -230,13 +230,18 @@ are 30, 90, and 180, respectively. With {quadratic} value=1, the numbers of columns are 930, 2790, and 5580, respectively. If the {quadratic} keyword value is set to 1, then additional -columns are appended to each per-atom array, corresponding to +columns are generated, corresponding to the products of all distinct pairs of bispectrum components. If the number of bispectrum components is {K}, then the number of distinct pairs -is {K}({K}+1)/2. These are output in subblocks of {K}({K}+1)/2 columns, using the same -ordering of sub-blocks as was used for the bispectrum -components. Within each sub-block, the ordering is upper-triangular, -(1,1),(1,2)...(1,{K}),(2,1)...({K}-1,{K}-1),({K}-1,{K}),({K},{K}) +is {K}({K}+1)/2. +For compute {sna/atom} these columns are appended to existing {K} columns. +The ordering of quadratic terms is upper-triangular, +(1,1),(1,2)...(1,{K}),(2,1)...({K}-1,{K}-1),({K}-1,{K}),({K},{K}). +For computes {snad/atom} and {snav/atom} each set of {K}({K}+1)/2 +additional columns is inserted directly after each of sub-block +of linear terms i.e. linear and quadratic terms are contiguous. +So the nesting order from inside to outside is bispectrum component, +linear then quadratic, vector/tensor component, type. These values can be accessed by any command that uses per-atom values from a compute as input. See "Section diff --git a/src/SNAP/compute_sna_atom.cpp b/src/SNAP/compute_sna_atom.cpp index 4b114d9ce7..9f56cee3f1 100644 --- a/src/SNAP/compute_sna_atom.cpp +++ b/src/SNAP/compute_sna_atom.cpp @@ -276,9 +276,13 @@ void ComputeSNAAtom::compute_peratom() for (int icoeff = 0; icoeff < ncoeff; icoeff++) { double bi = snaptr[tid]->bvec[icoeff]; + // diagonal element of quadratic matrix + + sna[i][ncount++] = 0.5*bi*bi; + // upper-triangular elements of quadratic matrix - for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++) + for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) sna[i][ncount++] = bi*snaptr[tid]->bvec[jcoeff]; } } diff --git a/src/SNAP/compute_snad_atom.cpp b/src/SNAP/compute_snad_atom.cpp index 54a6ce7612..586fb3d210 100644 --- a/src/SNAP/compute_snad_atom.cpp +++ b/src/SNAP/compute_snad_atom.cpp @@ -95,6 +95,11 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute snad/atom command"); rmin0 = atof(arg[iarg+1]); 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]); + iarg += 2; } else if (strcmp(arg[iarg],"switchflag") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute snad/atom command"); @@ -121,16 +126,11 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : } ncoeff = snaptr[0]->ncoeff; - twoncoeff = 2*ncoeff; - threencoeff = 3*ncoeff; - size_peratom_cols = threencoeff*atom->ntypes; - if (quadraticflag) { - ncoeffq = (ncoeff*(ncoeff+1))/2; - twoncoeffq = 2*ncoeffq; - threencoeffq = 3*ncoeffq; - size_peratom_cols += - threencoeffq*atom->ntypes; - } + nperdim = ncoeff; + if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2; + yoffset = nperdim; + zoffset = 2*nperdim; + size_peratom_cols = 3*nperdim*atom->ntypes; comm_reverse = size_peratom_cols; peratom_flag = 1; @@ -248,9 +248,10 @@ void ComputeSNADAtom::compute_peratom() const int* const jlist = firstneigh[i]; const int jnum = numneigh[i]; - const int typeoffset = threencoeff*(atom->type[i]-1); - const int quadraticoffset = threencoeff*atom->ntypes + - threencoeffq*(atom->type[i]-1); + // 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); // insure rij, inside, and typej are of size jnum @@ -304,16 +305,17 @@ void ComputeSNADAtom::compute_peratom() for (int icoeff = 0; icoeff < ncoeff; icoeff++) { snadi[icoeff] += snaptr[tid]->dbvec[icoeff][0]; - snadi[icoeff+ncoeff] += snaptr[tid]->dbvec[icoeff][1]; - snadi[icoeff+twoncoeff] += snaptr[tid]->dbvec[icoeff][2]; + snadi[icoeff+yoffset] += snaptr[tid]->dbvec[icoeff][1]; + snadi[icoeff+zoffset] += snaptr[tid]->dbvec[icoeff][2]; snadj[icoeff] -= snaptr[tid]->dbvec[icoeff][0]; - snadj[icoeff+ncoeff] -= snaptr[tid]->dbvec[icoeff][1]; - snadj[icoeff+twoncoeff] -= snaptr[tid]->dbvec[icoeff][2]; + snadj[icoeff+yoffset] -= snaptr[tid]->dbvec[icoeff][1]; + snadj[icoeff+zoffset] -= snaptr[tid]->dbvec[icoeff][2]; } if (quadraticflag) { - double *snadi = snad[i]+quadraticoffset; - double *snadj = snad[j]+quadraticoffset; + const int quadraticoffset = ncoeff; + snadi += quadraticoffset; + snadj += quadraticoffset; int ncount = 0; for (int icoeff = 0; icoeff < ncoeff; icoeff++) { double bi = snaptr[tid]->bvec[icoeff]; @@ -321,21 +323,36 @@ void ComputeSNADAtom::compute_peratom() double biy = snaptr[tid]->dbvec[icoeff][1]; double biz = snaptr[tid]->dbvec[icoeff][2]; + // diagonal elements of quadratic matrix + + double dbxtmp = bi*bix; + double dbytmp = bi*biy; + double dbztmp = bi*biz; + + snadi[ncount] += dbxtmp; + snadi[ncount+yoffset] += dbytmp; + snadi[ncount+zoffset] += dbztmp; + snadj[ncount] -= dbxtmp; + snadj[ncount+yoffset] -= dbytmp; + snadj[ncount+zoffset] -= dbztmp; + ncount++; + // upper-triangular elements of quadratic matrix - for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++) { + for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) { double dbxtmp = bi*snaptr[tid]->dbvec[jcoeff][0] + bix*snaptr[tid]->bvec[jcoeff]; double dbytmp = bi*snaptr[tid]->dbvec[jcoeff][1] + biy*snaptr[tid]->bvec[jcoeff]; double dbztmp = bi*snaptr[tid]->dbvec[jcoeff][2] + biz*snaptr[tid]->bvec[jcoeff]; - snadi[ncount] += dbxtmp; - snadi[ncount+ncoeffq] += dbytmp; - snadi[ncount+twoncoeffq] += dbztmp; - snadj[ncount] -= dbxtmp; - snadj[ncount+ncoeffq] -= dbytmp; - snadj[ncount+twoncoeffq] -= dbztmp; + + snadi[ncount] += dbxtmp; + snadi[ncount+yoffset] += dbytmp; + snadi[ncount+zoffset] += dbztmp; + snadj[ncount] -= dbxtmp; + snadj[ncount+yoffset] -= dbytmp; + snadj[ncount+zoffset] -= dbztmp; ncount++; } } @@ -361,7 +378,7 @@ int ComputeSNADAtom::pack_reverse_comm(int n, int first, double *buf) for (i = first; i < last; i++) for (icoeff = 0; icoeff < size_peratom_cols; icoeff++) buf[m++] = snad[i][icoeff]; - return comm_reverse; + return m; } /* ---------------------------------------------------------------------- */ @@ -387,8 +404,7 @@ double ComputeSNADAtom::memory_usage() double bytes = nmax*size_peratom_cols * sizeof(double); bytes += 3*njmax*sizeof(double); bytes += njmax*sizeof(int); - bytes += threencoeff*atom->ntypes; - if (quadraticflag) bytes += threencoeffq*atom->ntypes; + bytes += 3*nperdim*atom->ntypes; bytes += snaptr[0]->memory_usage()*comm->nthreads; return bytes; } diff --git a/src/SNAP/compute_snad_atom.h b/src/SNAP/compute_snad_atom.h index 5e40632d8c..a33e6047c2 100644 --- a/src/SNAP/compute_snad_atom.h +++ b/src/SNAP/compute_snad_atom.h @@ -37,7 +37,7 @@ class ComputeSNADAtom : public Compute { private: int nmax, njmax, diagonalstyle; - int ncoeff, twoncoeff, threencoeff, ncoeffq, twoncoeffq, threencoeffq; + int ncoeff, nperdim, yoffset, zoffset; double **cutsq; class NeighList *list; double **snad; diff --git a/src/SNAP/compute_snav_atom.cpp b/src/SNAP/compute_snav_atom.cpp index 0278be5a97..1248a223c3 100644 --- a/src/SNAP/compute_snav_atom.cpp +++ b/src/SNAP/compute_snav_atom.cpp @@ -96,6 +96,11 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute snav/atom command"); switchflag = atoi(arg[iarg+1]); 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]); + iarg += 2; } else if (strcmp(arg[iarg],"quadraticflag") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute snav/atom command"); @@ -117,22 +122,9 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : } ncoeff = snaptr[0]->ncoeff; - twoncoeff = 2*ncoeff; - threencoeff = 3*ncoeff; - fourncoeff = 4*ncoeff; - fivencoeff = 5*ncoeff; - sixncoeff = 6*ncoeff; - size_peratom_cols = sixncoeff*atom->ntypes; - if (quadraticflag) { - ncoeffq = ncoeff*ncoeff; - twoncoeffq = 2*ncoeffq; - threencoeffq = 3*ncoeffq; - fourncoeffq = 4*ncoeffq; - fivencoeffq = 5*ncoeffq; - sixncoeffq = 6*ncoeffq; - size_peratom_cols += - sixncoeffq*atom->ntypes; - } + nperdim = ncoeff; + if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2; + size_peratom_cols = 6*nperdim*atom->ntypes; comm_reverse = size_peratom_cols; peratom_flag = 1; @@ -251,9 +243,7 @@ void ComputeSNAVAtom::compute_peratom() const int* const jlist = firstneigh[i]; const int jnum = numneigh[i]; - const int typeoffset = sixncoeff*(atom->type[i]-1); - const int quadraticoffset = sixncoeff*atom->ntypes + - sixncoeffq*(atom->type[i]-1); + const int typeoffset = 6*nperdim*(atom->type[i]-1); // insure rij, inside, and typej are of size jnum @@ -307,23 +297,24 @@ void ComputeSNAVAtom::compute_peratom() double *snavj = snav[j]+typeoffset; for (int icoeff = 0; icoeff < ncoeff; icoeff++) { - snavi[icoeff] += snaptr[tid]->dbvec[icoeff][0]*xtmp; - snavi[icoeff+ncoeff] += snaptr[tid]->dbvec[icoeff][1]*ytmp; - snavi[icoeff+twoncoeff] += snaptr[tid]->dbvec[icoeff][2]*ztmp; - snavi[icoeff+threencoeff] += snaptr[tid]->dbvec[icoeff][1]*ztmp; - snavi[icoeff+fourncoeff] += snaptr[tid]->dbvec[icoeff][0]*ztmp; - snavi[icoeff+fivencoeff] += snaptr[tid]->dbvec[icoeff][0]*ytmp; - snavj[icoeff] -= snaptr[tid]->dbvec[icoeff][0]*x[j][0]; - snavj[icoeff+ncoeff] -= snaptr[tid]->dbvec[icoeff][1]*x[j][1]; - snavj[icoeff+twoncoeff] -= snaptr[tid]->dbvec[icoeff][2]*x[j][2]; - snavj[icoeff+threencoeff] -= snaptr[tid]->dbvec[icoeff][1]*x[j][2]; - snavj[icoeff+fourncoeff] -= snaptr[tid]->dbvec[icoeff][0]*x[j][2]; - snavj[icoeff+fivencoeff] -= snaptr[tid]->dbvec[icoeff][0]*x[j][1]; + snavi[icoeff] += snaptr[tid]->dbvec[icoeff][0]*xtmp; + snavi[icoeff+nperdim] += snaptr[tid]->dbvec[icoeff][1]*ytmp; + snavi[icoeff+2*nperdim] += snaptr[tid]->dbvec[icoeff][2]*ztmp; + snavi[icoeff+3*nperdim] += snaptr[tid]->dbvec[icoeff][1]*ztmp; + snavi[icoeff+4*nperdim] += snaptr[tid]->dbvec[icoeff][0]*ztmp; + snavi[icoeff+5*nperdim] += snaptr[tid]->dbvec[icoeff][0]*ytmp; + snavj[icoeff] -= snaptr[tid]->dbvec[icoeff][0]*x[j][0]; + snavj[icoeff+nperdim] -= snaptr[tid]->dbvec[icoeff][1]*x[j][1]; + snavj[icoeff+2*nperdim] -= snaptr[tid]->dbvec[icoeff][2]*x[j][2]; + snavj[icoeff+3*nperdim] -= snaptr[tid]->dbvec[icoeff][1]*x[j][2]; + snavj[icoeff+4*nperdim] -= snaptr[tid]->dbvec[icoeff][0]*x[j][2]; + snavj[icoeff+5*nperdim] -= snaptr[tid]->dbvec[icoeff][0]*x[j][1]; } if (quadraticflag) { - double *snavi = snav[i]+quadraticoffset; - double *snavj = snav[j]+quadraticoffset; + const int quadraticoffset = ncoeff; + snavi += quadraticoffset; + snavj += quadraticoffset; int ncount = 0; for (int icoeff = 0; icoeff < ncoeff; icoeff++) { double bi = snaptr[tid]->bvec[icoeff]; @@ -331,27 +322,46 @@ void ComputeSNAVAtom::compute_peratom() double biy = snaptr[tid]->dbvec[icoeff][1]; double biz = snaptr[tid]->dbvec[icoeff][2]; + // diagonal element of quadratic matrix + + double dbxtmp = bi*bix; + 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; + 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]; + ncount++; + // upper-triangular elements of quadratic matrix - for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++) { + for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) { double dbxtmp = bi*snaptr[tid]->dbvec[jcoeff][0] + bix*snaptr[tid]->bvec[jcoeff]; double dbytmp = bi*snaptr[tid]->dbvec[jcoeff][1] + biy*snaptr[tid]->bvec[jcoeff]; double dbztmp = bi*snaptr[tid]->dbvec[jcoeff][2] + biz*snaptr[tid]->bvec[jcoeff]; - snavi[ncount] += dbxtmp*xtmp; - snavi[ncount+ncoeffq] += dbytmp*ytmp; - snavi[ncount+twoncoeffq] += dbztmp*ztmp; - snavi[ncount+threencoeffq] += dbytmp*ztmp; - snavi[ncount+fourncoeffq] += dbxtmp*ztmp; - snavi[ncount+fivencoeffq] += dbxtmp*ytmp; - snavj[ncount] -= dbxtmp*x[j][0]; - snavj[ncount+ncoeffq] -= dbytmp*x[j][1]; - snavj[ncount+twoncoeffq] -= dbztmp*x[j][2]; - snavj[ncount+threencoeffq] -= dbytmp*x[j][2]; - snavj[ncount+fourncoeffq] -= dbxtmp*x[j][2]; - snavj[ncount+fivencoeffq] -= dbxtmp*x[j][1]; + 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; + 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]; ncount++; } } @@ -377,7 +387,7 @@ int ComputeSNAVAtom::pack_reverse_comm(int n, int first, double *buf) for (i = first; i < last; i++) for (icoeff = 0; icoeff < size_peratom_cols; icoeff++) buf[m++] = snav[i][icoeff]; - return comm_reverse; + return m; } /* ---------------------------------------------------------------------- */ @@ -403,8 +413,8 @@ double ComputeSNAVAtom::memory_usage() double bytes = nmax*size_peratom_cols * sizeof(double); bytes += 3*njmax*sizeof(double); bytes += njmax*sizeof(int); - bytes += sixncoeff*atom->ntypes; - if (quadraticflag) bytes += sixncoeffq*atom->ntypes; + bytes += 6*nperdim*atom->ntypes; + if (quadraticflag) bytes += 6*nperdim*atom->ntypes; bytes += snaptr[0]->memory_usage()*comm->nthreads; return bytes; } diff --git a/src/SNAP/compute_snav_atom.h b/src/SNAP/compute_snav_atom.h index 35f1478393..2eb9fb804f 100644 --- a/src/SNAP/compute_snav_atom.h +++ b/src/SNAP/compute_snav_atom.h @@ -37,8 +37,7 @@ class ComputeSNAVAtom : public Compute { private: int nmax, njmax, diagonalstyle; - int ncoeff, twoncoeff, threencoeff, fourncoeff, fivencoeff, sixncoeff; - int ncoeffq, twoncoeffq, threencoeffq, fourncoeffq, fivencoeffq, sixncoeffq; + int ncoeff, nperdim; double **cutsq; class NeighList *list; double **snav; diff --git a/src/SNAP/pair_snap.cpp b/src/SNAP/pair_snap.cpp index 377235685c..8714d55c5c 100644 --- a/src/SNAP/pair_snap.cpp +++ b/src/SNAP/pair_snap.cpp @@ -278,14 +278,15 @@ void PairSNAP::compute_regular(int eflag, int vflag) double bveci = snaptr->bvec[icoeff]; double fack = coeffi[k]*bveci; double* dbveci = snaptr->dbvec[icoeff]; - fij[0] += fack*snaptr->dbvec[icoeff][0]; - fij[1] += fack*snaptr->dbvec[icoeff][1]; - fij[2] += fack*snaptr->dbvec[icoeff][2]; + fij[0] += fack*dbveci[0]; + fij[1] += fack*dbveci[1]; + fij[2] += fack*dbveci[2]; k++; for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) { double facki = coeffi[k]*bveci; double fackj = coeffi[k]*snaptr->bvec[jcoeff]; double* dbvecj = snaptr->dbvec[jcoeff]; + fij[0] += facki*dbvecj[0]+fackj*dbveci[0]; fij[1] += facki*dbvecj[1]+fackj*dbveci[1]; fij[2] += facki*dbvecj[2]+fackj*dbveci[2]; @@ -1529,10 +1530,10 @@ void PairSNAP::coeff(int narg, char **arg) } if (comm->me == 0) - printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff); - if (ncoeff != sna[0]->ncoeff) { - error->all(FLERR,"Incorrect SNAP parameter file"); - } + if (ncoeff != sna[0]->ncoeff) { + printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff); + error->all(FLERR,"Incorrect SNAP parameter file"); + } // Calculate maximum cutoff for all elements diff --git a/src/SNAP/pair_snap.h b/src/SNAP/pair_snap.h index d39cb0f8d4..b60ab3c3e4 100644 --- a/src/SNAP/pair_snap.h +++ b/src/SNAP/pair_snap.h @@ -37,11 +37,8 @@ public: virtual double init_one(int, int); virtual double memory_usage(); - double rcutfac, quadraticflag; // declared public to workaround gcc 4.9 - int ncoeff; // compiler bug, manifest in KOKKOS package - protected: - int ncoeffq, ncoeffall; + int ncoeff, ncoeffq, ncoeffall; double **bvec, ***dbvec; class SNA** sna; int nmax; @@ -100,8 +97,8 @@ protected: double *wjelem; // elements weights double **coeffelem; // element bispectrum coefficients int *map; // mapping from atom types to elements - int twojmax, diagonalstyle, switchflag, bzeroflag; - double rfac0, rmin0, wj1, wj2; + int twojmax, diagonalstyle, switchflag, bzeroflag, quadraticflag; + double rcutfac, rfac0, rmin0, wj1, wj2; int rcutfacflag, twojmaxflag; // flags for required parameters };