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

updates to quadratic form of SNAP potential

This commit is contained in:
Steven J. Plimpton 2018-05-03 10:23:50 -06:00
parent d5ec76290b
commit 75cffc78c5
8 changed files with 135 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
doc/src/compute_sna_atom.txt
View File

@ -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

6
src/SNAP/compute_sna_atom.cpp
View File

@ -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];
}
}

74
src/SNAP/compute_snad_atom.cpp
View File

@ -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;
}

2
src/SNAP/compute_snad_atom.h
View File

@ -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;

108
src/SNAP/compute_snav_atom.cpp
View File

@ -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;
}

3
src/SNAP/compute_snav_atom.h
View File

@ -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;

15
src/SNAP/pair_snap.cpp
View File

@ -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

9
src/SNAP/pair_snap.h
View File

@ -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
};