diff --git a/doc/src/compute_sna_atom.txt b/doc/src/compute_sna_atom.txt
index f82df0d81601eb16d6f0ba578ab6783cc735f7a5..defd98d520d06489d8e93e93369cb33bdfb4b2d5 100644
--- a/doc/src/compute_sna_atom.txt
+++ b/doc/src/compute_sna_atom.txt
@@ -231,11 +231,12 @@ 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
-a matrix of quantities that are products of two bispectrum components. If the
-number of bispectrum components is {K}, then the number of matrix elements
-is {K}^2. These are output in subblocks of {K}^2 columns, using the same
-ordering of columns and sub-blocks as was used for the bispectrum
-components.
+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})
These values can be accessed by any command that uses per-atom values
from a compute as input. See "Section
diff --git a/doc/src/pair_snap.txt b/doc/src/pair_snap.txt
index 4e3546481dd02baa796f6b6657da7d404b245626..fa90dc34e92a2b30390317dcce77f5c0d55f5547 100644
--- a/doc/src/pair_snap.txt
+++ b/doc/src/pair_snap.txt
@@ -140,10 +140,15 @@ The default values for these keywords are
{rmin0} = 0.0
{diagonalstyle} = 3
{switchflag} = 0
-{bzeroflag} = 1 :ul
+{bzeroflag} = 1
+{quadraticflag} = 1 :ul
-Detailed definitions of these keywords are given on the "compute
+Detailed definitions for all the keywords are given on the "compute
sna/atom"_compute_sna_atom.html doc page.
+If {quadraticflag} is set to 1, then the SNAP energy expression includes the quadratic term,
+0.5*B^t.alpha.B, where alpha is a symmetric {K} by {K} matrix.
+The SNAP element file should contain {K}({K}+1)/2 additional coefficients
+for each element, the upper-triangular elements of alpha.
:line
diff --git a/examples/COUPLE/simple/simple.cpp b/examples/COUPLE/simple/simple.cpp
index b279ae704cb9febd3401fe420ec0a7db292c5189..894c708978176747a095e681245f325991161f3f 100644
--- a/examples/COUPLE/simple/simple.cpp
+++ b/examples/COUPLE/simple/simple.cpp
@@ -153,7 +153,7 @@ int main(int narg, char **arg)
for (int i = 0; i < natoms; i++) type[i] = 1;
lmp->input->one("delete_atoms group all");
- lammps_create_atoms(lmp,natoms,NULL,type,x,v);
+ lammps_create_atoms(lmp,natoms,NULL,type,x,v,NULL,0);
lmp->input->one("run 10");
}
diff --git a/examples/snap/Ta06A.snapparam b/examples/snap/Ta06A.snapparam
index 7b30312f59abbbb1759580e35f0279fdf1cc1c98..283629d6581651d070f84e785534c87ede8391aa 100644
--- a/examples/snap/Ta06A.snapparam
+++ b/examples/snap/Ta06A.snapparam
@@ -8,8 +8,8 @@ twojmax 6
# optional
-gamma 1
rfac0 0.99363
rmin0 0
diagonalstyle 3
bzeroflag 0
+quadraticflag 0
diff --git a/examples/snap/W_2940_2017_2.snap b/examples/snap/W_2940_2017_2.snap
index 51eee41a0a00be9c1f34ec77c5b8e25be4c70924..04b8d58094ea16e58668740f3b436ddbb813cc18 100644
--- a/examples/snap/W_2940_2017_2.snap
+++ b/examples/snap/W_2940_2017_2.snap
@@ -5,7 +5,7 @@ variable zblcutinner equal 4
variable zblcutouter equal 4.8
variable zblz equal 74
-# Specify hybrid with SNAP, ZBL, and long-range Coulomb
+# Specify hybrid with SNAP and ZBL
pair_style hybrid/overlay &
zbl ${zblcutinner} ${zblcutouter} &
diff --git a/examples/snap/W_2940_2017_2.snapparam b/examples/snap/W_2940_2017_2.snapparam
index f17961bdd737f379302feca3b87a1f1c60690226..27ab61a266922fd37687c73ea6f212385f9e1a79 100644
--- a/examples/snap/W_2940_2017_2.snapparam
+++ b/examples/snap/W_2940_2017_2.snapparam
@@ -6,8 +6,8 @@ twojmax 8
# optional
-gamma 1
rfac0 0.99363
rmin0 0
diagonalstyle 3
bzeroflag 0
+quadraticflag 0
diff --git a/examples/snap/W_2940_2017_2_He_JW2013.snap b/examples/snap/W_2940_2017_2_He_JW2013.snap
index 45a31955b3a1c21138bcb6a62dbd877804909fdc..cb70916ec48dc6b9ce6d4bd8f91302a199814889 100644
--- a/examples/snap/W_2940_2017_2_He_JW2013.snap
+++ b/examples/snap/W_2940_2017_2_He_JW2013.snap
@@ -5,7 +5,7 @@ variable zblcutinner equal 4
variable zblcutouter equal 4.8
variable zblz equal 74
-# Specify hybrid with SNAP, ZBL, and long-range Coulomb
+# Specify hybrid with SNAP and ZBL
pair_style hybrid/overlay zbl ${zblcutinner} ${zblcutouter} snap table spline 10000 table spline 10000
pair_coeff 1 1 zbl ${zblz} ${zblz}
diff --git a/src/SNAP/compute_sna_atom.cpp b/src/SNAP/compute_sna_atom.cpp
index cba6fae9b75fcd51c3124ef39ff78d4e4b0849ce..5341d16efa879e5a770dd0568cb634bf843be194 100644
--- a/src/SNAP/compute_sna_atom.cpp
+++ b/src/SNAP/compute_sna_atom.cpp
@@ -129,7 +129,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
ncoeff = snaptr[0]->ncoeff;
size_peratom_cols = ncoeff;
- if (quadraticflag) size_peratom_cols += ncoeff*ncoeff;
+ if (quadraticflag) size_peratom_cols += (ncoeff*(ncoeff+1))/2;
peratom_flag = 1;
nmax = 0;
@@ -275,7 +275,10 @@ void ComputeSNAAtom::compute_peratom()
int ncount = ncoeff;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bi = snaptr[tid]->bvec[icoeff];
- for (int jcoeff = 0; jcoeff < ncoeff; jcoeff++)
+
+ // upper-triangular elements of quadratic matrix
+
+ for (int jcoeff = icoeff; 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 39f34dd8cd50e8b15edace89accb10d8026e1178..c2ac5f5afd3f9ec4230fcc37cfc52cfd3e81ac0c 100644
--- a/src/SNAP/compute_snad_atom.cpp
+++ b/src/SNAP/compute_snad_atom.cpp
@@ -125,11 +125,11 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
threencoeff = 3*ncoeff;
size_peratom_cols = threencoeff*atom->ntypes;
if (quadraticflag) {
- ncoeffsq = ncoeff*ncoeff;
- twoncoeffsq = 2*ncoeffsq;
- threencoeffsq = 3*ncoeffsq;
+ ncoeffq = (ncoeff*(ncoeff+1))/2;
+ twoncoeffq = 2*ncoeffq;
+ threencoeffq = 3*ncoeffq;
size_peratom_cols +=
- threencoeffsq*atom->ntypes;
+ threencoeffq*atom->ntypes;
}
comm_reverse = size_peratom_cols;
peratom_flag = 1;
@@ -250,7 +250,7 @@ void ComputeSNADAtom::compute_peratom()
const int typeoffset = threencoeff*(atom->type[i]-1);
const int quadraticoffset = threencoeff*atom->ntypes +
- threencoeffsq*(atom->type[i]-1);
+ threencoeffq*(atom->type[i]-1);
// insure rij, inside, and typej are of size jnum
@@ -320,7 +320,10 @@ void ComputeSNADAtom::compute_peratom()
double bix = snaptr[tid]->dbvec[icoeff][0];
double biy = snaptr[tid]->dbvec[icoeff][1];
double biz = snaptr[tid]->dbvec[icoeff][2];
- for (int jcoeff = 0; jcoeff < ncoeff; jcoeff++) {
+
+ // upper-triangular elements of quadratic matrix
+
+ for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++) {
double dbxtmp = bi*snaptr[tid]->dbvec[jcoeff][0]
+ bix*snaptr[tid]->bvec[jcoeff];
double dbytmp = bi*snaptr[tid]->dbvec[jcoeff][1]
@@ -328,11 +331,11 @@ void ComputeSNADAtom::compute_peratom()
double dbztmp = bi*snaptr[tid]->dbvec[jcoeff][2]
+ biz*snaptr[tid]->bvec[jcoeff];
snadi[ncount] += dbxtmp;
- snadi[ncount+ncoeffsq] += dbytmp;
- snadi[ncount+twoncoeffsq] += dbztmp;
+ snadi[ncount+ncoeffq] += dbytmp;
+ snadi[ncount+twoncoeffq] += dbztmp;
snadj[ncount] -= dbxtmp;
- snadj[ncount+ncoeffsq] -= dbytmp;
- snadj[ncount+twoncoeffsq] -= dbztmp;
+ snadj[ncount+ncoeffq] -= dbytmp;
+ snadj[ncount+twoncoeffq] -= dbztmp;
ncount++;
}
}
@@ -385,7 +388,7 @@ double ComputeSNADAtom::memory_usage()
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
bytes += threencoeff*atom->ntypes;
- if (quadraticflag) bytes += threencoeffsq*atom->ntypes;
+ if (quadraticflag) bytes += threencoeffq*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 0d5a369ab6b9050bea1104c2d928bc55b388e886..8eaae2a67f9b5a2f4a2898604eb9be9f056542a6 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, ncoeffsq, twoncoeffsq, threencoeffsq;
+ int ncoeff, twoncoeff, threencoeff, ncoeffq, twoncoeffq, threencoeffq;
double **cutsq;
class NeighList *list;
double **snad;
diff --git a/src/SNAP/compute_snav_atom.cpp b/src/SNAP/compute_snav_atom.cpp
index 0d21d16561956cd1424c5b267067a471f93bd7bb..3b5383ddf45f68f8af7186d7c5509c40015f61e1 100644
--- a/src/SNAP/compute_snav_atom.cpp
+++ b/src/SNAP/compute_snav_atom.cpp
@@ -124,14 +124,14 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
sixncoeff = 6*ncoeff;
size_peratom_cols = sixncoeff*atom->ntypes;
if (quadraticflag) {
- ncoeffsq = ncoeff*ncoeff;
- twoncoeffsq = 2*ncoeffsq;
- threencoeffsq = 3*ncoeffsq;
- fourncoeffsq = 4*ncoeffsq;
- fivencoeffsq = 5*ncoeffsq;
- sixncoeffsq = 6*ncoeffsq;
+ ncoeffq = ncoeff*ncoeff;
+ twoncoeffq = 2*ncoeffq;
+ threencoeffq = 3*ncoeffq;
+ fourncoeffq = 4*ncoeffq;
+ fivencoeffq = 5*ncoeffq;
+ sixncoeffq = 6*ncoeffq;
size_peratom_cols +=
- sixncoeffsq*atom->ntypes;
+ sixncoeffq*atom->ntypes;
}
comm_reverse = size_peratom_cols;
peratom_flag = 1;
@@ -253,7 +253,7 @@ void ComputeSNAVAtom::compute_peratom()
const int typeoffset = sixncoeff*(atom->type[i]-1);
const int quadraticoffset = sixncoeff*atom->ntypes +
- sixncoeffsq*(atom->type[i]-1);
+ sixncoeffq*(atom->type[i]-1);
// insure rij, inside, and typej are of size jnum
@@ -330,7 +330,10 @@ void ComputeSNAVAtom::compute_peratom()
double bix = snaptr[tid]->dbvec[icoeff][0];
double biy = snaptr[tid]->dbvec[icoeff][1];
double biz = snaptr[tid]->dbvec[icoeff][2];
- for (int jcoeff = 0; jcoeff < ncoeff; jcoeff++) {
+
+ // upper-triangular elements of quadratic matrix
+
+ for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++) {
double dbxtmp = bi*snaptr[tid]->dbvec[jcoeff][0]
+ bix*snaptr[tid]->bvec[jcoeff];
double dbytmp = bi*snaptr[tid]->dbvec[jcoeff][1]
@@ -338,17 +341,17 @@ void ComputeSNAVAtom::compute_peratom()
double dbztmp = bi*snaptr[tid]->dbvec[jcoeff][2]
+ biz*snaptr[tid]->bvec[jcoeff];
snavi[ncount] += dbxtmp*xtmp;
- snavi[ncount+ncoeffsq] += dbytmp*ytmp;
- snavi[ncount+twoncoeffsq] += dbztmp*ztmp;
- snavi[ncount+threencoeffsq] += dbytmp*ztmp;
- snavi[ncount+fourncoeffsq] += dbxtmp*ztmp;
- snavi[ncount+fivencoeffsq] += dbxtmp*ytmp;
+ 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+ncoeffsq] -= dbytmp*x[j][1];
- snavj[ncount+twoncoeffsq] -= dbztmp*x[j][2];
- snavj[ncount+threencoeffsq] -= dbytmp*x[j][2];
- snavj[ncount+fourncoeffsq] -= dbxtmp*x[j][2];
- snavj[ncount+fivencoeffsq] -= dbxtmp*x[j][1];
+ 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];
ncount++;
}
}
@@ -401,7 +404,7 @@ double ComputeSNAVAtom::memory_usage()
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
bytes += sixncoeff*atom->ntypes;
- if (quadraticflag) bytes += sixncoeffsq*atom->ntypes;
+ if (quadraticflag) bytes += sixncoeffq*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 33ae4f92173945c9754b0ef17b5a5b9e730d4a79..35f1478393520dc9f2a687a0db86201b0ca5e52c 100644
--- a/src/SNAP/compute_snav_atom.h
+++ b/src/SNAP/compute_snav_atom.h
@@ -38,7 +38,7 @@ class ComputeSNAVAtom : public Compute {
private:
int nmax, njmax, diagonalstyle;
int ncoeff, twoncoeff, threencoeff, fourncoeff, fivencoeff, sixncoeff;
- int ncoeffsq, twoncoeffsq, threencoeffsq, fourncoeffsq, fivencoeffsq, sixncoeffsq;
+ int ncoeffq, twoncoeffq, threencoeffq, fourncoeffq, fivencoeffq, sixncoeffq;
double **cutsq;
class NeighList *list;
double **snav;
diff --git a/src/SNAP/pair_snap.cpp b/src/SNAP/pair_snap.cpp
index e4ed57b933cba2d65e310ff62238b6d9a92b1660..492e9ad6436ecb8e18ae794127528695ec85c5c4 100644
--- a/src/SNAP/pair_snap.cpp
+++ b/src/SNAP/pair_snap.cpp
@@ -35,6 +35,10 @@ 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)
@@ -232,10 +236,6 @@ void PairSNAP::compute_regular(int eflag, int vflag)
snaptr->compute_ui(ninside);
snaptr->compute_zi();
- if (!gammaoneflag) {
- snaptr->compute_bi();
- snaptr->copy_bi2bvec();
- }
// for neighbors of I within cutoff:
// compute dUi/drj and dBi/drj
@@ -255,17 +255,41 @@ void PairSNAP::compute_regular(int eflag, int vflag)
fij[1] = 0.0;
fij[2] = 0.0;
+ // linear contributions
+
for (int k = 1; k <= ncoeff; k++) {
- double bgb;
- if (gammaoneflag)
- bgb = coeffi[k];
- else bgb = coeffi[k]*
- gamma*pow(snaptr->bvec[k-1],gamma-1.0);
+ double bgb = coeffi[k];
fij[0] += bgb*snaptr->dbvec[k-1][0];
fij[1] += bgb*snaptr->dbvec[k-1][1];
fij[2] += bgb*snaptr->dbvec[k-1][2];
}
+ // quadratic contributions
+
+ if (quadraticflag) {
+ snaptr->compute_bi();
+ snaptr->copy_bi2bvec();
+ int k = ncoeff+1;
+ for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+ 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];
+ 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];
+ k++;
+ }
+ }
+ }
+
f[i][0] += fij[0];
f[i][1] += fij[1];
f[i][2] += fij[2];
@@ -285,14 +309,33 @@ void PairSNAP::compute_regular(int eflag, int vflag)
// evdwl = energy of atom I, sum over coeffs_k * Bi_k
evdwl = coeffi[0];
- if (gammaoneflag) {
- snaptr->compute_bi();
- snaptr->copy_bi2bvec();
- for (int k = 1; k <= ncoeff; k++)
- evdwl += coeffi[k]*snaptr->bvec[k-1];
- } else
- for (int k = 1; k <= ncoeff; k++)
- evdwl += coeffi[k]*pow(snaptr->bvec[k-1],gamma);
+ if (!quadraticflag) {
+ snaptr->compute_bi();
+ snaptr->copy_bi2bvec();
+ }
+
+ // 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 = 1; k <= ncoeff; k++)
+ evdwl += coeffi[k]*snaptr->bvec[k-1];
+
+ // quadratic contributions
+
+ if (quadraticflag) {
+ int k = ncoeff+1;
+ for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+ double bveci = snaptr->bvec[icoeff];
+ evdwl += 0.5*coeffi[k++]*bveci*bveci;
+ for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+ evdwl += coeffi[k++]*bveci*snaptr->bvec[jcoeff];
+ }
+ }
+ }
ev_tally_full(i,2.0*evdwl,0.0,0.0,delx,dely,delz);
}
@@ -562,26 +605,46 @@ void PairSNAP::compute_optimized(int eflag, int vflag)
sna[tid]->compute_dbidrj();
sna[tid]->copy_dbi2dbvec();
- if (!gammaoneflag) {
- sna[tid]->compute_bi();
- sna[tid]->copy_bi2bvec();
- }
fij[0] = 0.0;
fij[1] = 0.0;
fij[2] = 0.0;
+ // linear contributions
+
for (k = 1; k <= ncoeff; k++) {
- double bgb;
- if (gammaoneflag)
- bgb = coeffi[k];
- else bgb = coeffi[k]*
- gamma*pow(sna[tid]->bvec[k-1],gamma-1.0);
+ double bgb = coeffi[k];
fij[0] += bgb*sna[tid]->dbvec[k-1][0];
fij[1] += bgb*sna[tid]->dbvec[k-1][1];
fij[2] += bgb*sna[tid]->dbvec[k-1][2];
}
+ // quadratic contributions
+
+ if (quadraticflag) {
+ sna[tid]->compute_bi();
+ sna[tid]->copy_bi2bvec();
+ int k = ncoeff+1;
+ for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+ double bveci = sna[tid]->bvec[icoeff];
+ double fack = coeffi[k]*bveci;
+ double* dbveci = sna[tid]->dbvec[icoeff];
+ fij[0] += fack*sna[tid]->dbvec[icoeff][0];
+ fij[1] += fack*sna[tid]->dbvec[icoeff][1];
+ fij[2] += fack*sna[tid]->dbvec[icoeff][2];
+ k++;
+ for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+ double facki = coeffi[k]*bveci;
+ double fackj = coeffi[k]*sna[tid]->bvec[jcoeff];
+ double* dbvecj = sna[tid]->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];
+ k++;
+ }
+ }
+ }
+
#if defined(_OPENMP)
#pragma omp critical
#endif
@@ -606,15 +669,33 @@ void PairSNAP::compute_optimized(int eflag, int vflag)
if (eflag&&pairs[iijj][1] == 0) {
evdwl = coeffi[0];
- if (gammaoneflag) {
- sna[tid]->compute_bi();
- sna[tid]->copy_bi2bvec();
- for (int k = 1; k <= ncoeff; k++)
- evdwl += coeffi[k]*sna[tid]->bvec[k-1];
- } else
- for (int k = 1; k <= ncoeff; k++)
- evdwl += coeffi[k]*pow(sna[tid]->bvec[k-1],gamma);
+ sna[tid]->compute_bi();
+ sna[tid]->copy_bi2bvec();
+
+ // 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 = 1; k <= ncoeff; k++)
+ evdwl += coeffi[k]*sna[tid]->bvec[k-1];
+
+ // quadratic contributions
+
+ if (quadraticflag) {
+ int k = ncoeff+1;
+ for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+ double bveci = sna[tid]->bvec[icoeff];
+ evdwl += 0.5*coeffi[k++]*bveci*bveci;
+ for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+ evdwl += coeffi[k++]*bveci*sna[tid]->bvec[jcoeff];
+ }
+ }
+ }
+
#if defined(_OPENMP)
#pragma omp critical
#endif
@@ -1363,6 +1444,22 @@ void PairSNAP::coeff(int narg, char **arg)
read_files(coefffilename,paramfilename);
+ if (!quadraticflag)
+ ncoeff = ncoeffall - 1;
+ else {
+
+ // ncoeffall should be (ncoeff+2)*(ncoeff+1)/2
+ // so, ncoeff = floor(sqrt(2*ncoeffall))-1
+
+ ncoeff = sqrt(2*ncoeffall)-1;
+ ncoeffq = (ncoeff*(ncoeff+1))/2;
+ int ntmp = 1+ncoeff+ncoeffq;
+ if (ntmp != ncoeffall) {
+ printf("ncoeffall = %d ntmp = %d ncoeff = %d \n",ncoeffall,ntmp,ncoeff);
+ error->all(FLERR,"Incorrect SNAP coeff file");
+ }
+ }
+
// read args that map atom types to SNAP elements
// map[i] = which element the Ith atom type is, -1 if not mapped
// map[0] is not used
@@ -1416,6 +1513,7 @@ void PairSNAP::coeff(int narg, char **arg)
sna[tid]->grow_rij(nmax);
}
+ printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff);
if (ncoeff != sna[0]->ncoeff) {
printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff);
error->all(FLERR,"Incorrect SNAP parameter file");
@@ -1470,7 +1568,7 @@ double PairSNAP::init_one(int i, int j)
void PairSNAP::read_files(char *coefffilename, char *paramfilename)
{
- // open SNAP ceofficient file on proc 0
+ // open SNAP coefficient file on proc 0
FILE *fpcoeff;
if (comm->me == 0) {
@@ -1518,13 +1616,13 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
words[iword] = strtok(NULL,"' \t\n\r\f");
int nelemfile = atoi(words[0]);
- ncoeff = atoi(words[1])-1;
-
+ ncoeffall = atoi(words[1]);
+
// Set up element lists
memory->create(radelem,nelements,"pair:radelem");
memory->create(wjelem,nelements,"pair:wjelem");
- memory->create(coeffelem,nelements,ncoeff+1,"pair:coeffelem");
+ memory->create(coeffelem,nelements,ncoeffall,"pair:coeffelem");
int *found = new int[nelements];
for (int ielem = 0; ielem < nelements; ielem++)
@@ -1569,7 +1667,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
if (strcmp(elemtmp,elements[ielem]) == 0) break;
if (ielem == nelements) {
if (comm->me == 0)
- for (int icoeff = 0; icoeff <= ncoeff; icoeff++)
+ for (int icoeff = 0; icoeff < ncoeffall; icoeff++)
ptr = fgets(line,MAXLINE,fpcoeff);
continue;
}
@@ -1578,7 +1676,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
if (found[ielem]) {
if (comm->me == 0)
- for (int icoeff = 0; icoeff <= ncoeff; icoeff++)
+ for (int icoeff = 0; icoeff < ncoeffall; icoeff++)
ptr = fgets(line,MAXLINE,fpcoeff);
continue;
}
@@ -1595,7 +1693,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
elements[ielem], radelem[ielem], wjelem[ielem]);
}
- for (int icoeff = 0; icoeff <= ncoeff; icoeff++) {
+ for (int icoeff = 0; icoeff < ncoeffall; icoeff++) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fpcoeff);
if (ptr == NULL) {
@@ -1629,13 +1727,12 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
// Set defaults for optional keywords
- gamma = 1.0;
- gammaoneflag = 1;
rfac0 = 0.99363;
rmin0 = 0.0;
diagonalstyle = 3;
switchflag = 1;
bzeroflag = 1;
+ quadraticflag = 0;
// open SNAP parameter file on proc 0
@@ -1689,9 +1786,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
} else if (strcmp(keywd,"twojmax") == 0) {
twojmax = atoi(keyval);
twojmaxflag = 1;
- } else if (strcmp(keywd,"gamma") == 0)
- gamma = atof(keyval);
- else if (strcmp(keywd,"rfac0") == 0)
+ } else if (strcmp(keywd,"rfac0") == 0)
rfac0 = atof(keyval);
else if (strcmp(keywd,"rmin0") == 0)
rmin0 = atof(keyval);
@@ -1701,6 +1796,8 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
switchflag = atoi(keyval);
else if (strcmp(keywd,"bzeroflag") == 0)
bzeroflag = atoi(keyval);
+ else if (strcmp(keywd,"quadraticflag") == 0)
+ quadraticflag = atoi(keyval);
else
error->all(FLERR,"Incorrect SNAP parameter file");
}
@@ -1708,9 +1805,6 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
if (rcutfacflag == 0 || twojmaxflag == 0)
error->all(FLERR,"Incorrect SNAP parameter file");
- if (gamma == 1.0) gammaoneflag = 1;
- else gammaoneflag = 0;
-
delete[] found;
}
@@ -1726,7 +1820,7 @@ double PairSNAP::memory_usage()
bytes += n*n*sizeof(double);
bytes += 3*nmax*sizeof(double);
bytes += nmax*sizeof(int);
- bytes += (2*ncoeff+1)*sizeof(double);
+ bytes += (2*ncoeffall)*sizeof(double);
bytes += (ncoeff*3)*sizeof(double);
bytes += sna[0]->memory_usage()*nthreads;
return bytes;
diff --git a/src/SNAP/pair_snap.h b/src/SNAP/pair_snap.h
index 559d3ef5717f46d0f75bb0518010b46dd30700c8..9dec211e8eb4adaf7b682e240adc7be548f8e301 100644
--- a/src/SNAP/pair_snap.h
+++ b/src/SNAP/pair_snap.h
@@ -38,7 +38,7 @@ public:
double memory_usage();
protected:
- int ncoeff;
+ int ncoeff, ncoeffq, ncoeffall;
double **bvec, ***dbvec;
class SNA** sna;
int nmax;
@@ -89,7 +89,6 @@ protected:
// timespec starttime, endtime;
double timers[4];
#endif
- double gamma;
double rcutmax; // max cutoff for all elements
int nelements; // # of unique elements
@@ -98,10 +97,9 @@ protected:
double *wjelem; // elements weights
double **coeffelem; // element bispectrum coefficients
int *map; // mapping from atom types to elements
- int twojmax, diagonalstyle, switchflag, bzeroflag;
+ int twojmax, diagonalstyle, switchflag, bzeroflag, quadraticflag;
double rcutfac, rfac0, rmin0, wj1, wj2;
int rcutfacflag, twojmaxflag; // flags for required parameters
- int gammaoneflag; // 1 if parameter gamma is 1
};
}
diff --git a/src/npair_full_bin_atomonly.h b/src/npair_full_bin_atomonly.h
index 7ca5d667c7105b12dc01e75e4143e977adb2cc4d..f8c33d95588c3e814ab4b96190c113c7505085b0 100644
--- a/src/npair_full_bin_atomonly.h
+++ b/src/npair_full_bin_atomonly.h
@@ -14,7 +14,7 @@
#ifdef NPAIR_CLASS
NPairStyle(full/bin/atomonly,
- NPairFullBin,
+ NPairFullBinAtomonly,
NP_FULL | NP_BIN | NP_ATOMONLY |
NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI)