diff --git a/src/SNAP/compute_sna_atom.cpp b/src/SNAP/compute_sna_atom.cpp
index 9a0d8a0a6b25ccadb5ed8edd7c384ab95c241357..b77e01c5e2e44be7b7d088ecbb03fa82affcf7cf 100644
--- a/src/SNAP/compute_sna_atom.cpp
+++ b/src/SNAP/compute_sna_atom.cpp
@@ -102,17 +102,16 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
} else error->all(FLERR,"Illegal compute sna/atom command");
}
- // always unset use_shared_arrays
-
- int use_shared_arrays = 0;
- int nthreads = omp_get_max_threads();
-
- snaptr = new SNA*[nthreads];
-#pragma omp parallel shared(rfac0,twojmax,rmin0)
+ snaptr = new SNA*[comm->nthreads];
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag)
+#endif
{
int tid = omp_get_thread_num();
- snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,use_shared_arrays,
- rmin0,switchflag);
+
+ // always unset use_shared_arrays since it does not work with computes
+ snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,
+ 0 /*use_shared_arrays*/, rmin0,switchflag);
}
ncoeff = snaptr[0]->ncoeff;
@@ -160,7 +159,9 @@ void ComputeSNAAtom::init()
if (strcmp(modify->compute[i]->style,"sna/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute sna/atom");
- #pragma omp parallel
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
{
int tid = omp_get_thread_num();
snaptr[tid]->init();
@@ -197,7 +198,7 @@ void ComputeSNAAtom::compute_peratom()
const int* const ilist = list->ilist;
const int* const numneigh = list->numneigh;
int** const firstneigh = list->firstneigh;
- int *type = atom->type;
+ int * const type = atom->type;
// compute sna for each atom in group
// use full neighbor list to count atoms less than cutoff
@@ -205,7 +206,9 @@ void ComputeSNAAtom::compute_peratom()
double** const x = atom->x;
const int* const mask = atom->mask;
- #pragma omp parallel for
+#if defined(_OPENMP)
+#pragma omp parallel for default(none)
+#endif
for (int ii = 0; ii < inum; ii++) {
const int tid = omp_get_thread_num();
const int i = ilist[ii];
@@ -270,7 +273,7 @@ double ComputeSNAAtom::memory_usage()
double bytes = nmax*size_peratom_cols * sizeof(double);
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
- bytes += snaptr[0]->memory_usage()*omp_get_max_threads();
+ bytes += snaptr[0]->memory_usage()*comm->nthreads;
return bytes;
}
diff --git a/src/SNAP/compute_snad_atom.cpp b/src/SNAP/compute_snad_atom.cpp
index 79eb3838ca110efaf4335cdcc1d3d21d5bce9cf3..7039d0dc4b6e50bffc39c5b68d9947a4274845a4 100644
--- a/src/SNAP/compute_snad_atom.cpp
+++ b/src/SNAP/compute_snad_atom.cpp
@@ -45,7 +45,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
// default values
diagonalstyle = 0;
- rmin0 = 0;
+ rmin0 = 0.0;
switchflag = 1;
// process required arguments
@@ -94,15 +94,17 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else error->all(FLERR,"Illegal compute snad/atom command");
}
- // use_shared_arrays doesn't work with computes and should be disabled.
- int use_shared_arrays = 0;
- int nthreads = omp_get_max_threads();
- snaptr = new SNA*[nthreads];
-#pragma omp parallel shared(rfac0,twojmax,rmin0)
+
+ snaptr = new SNA*[comm->nthreads];
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag)
+#endif
{
int tid = omp_get_thread_num();
- snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,use_shared_arrays,
- rmin0,switchflag);
+
+ // always unset use_shared_arrays since it does not work with computes
+ snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,
+ 0 /*use_shared_arrays*/, rmin0,switchflag);
}
ncoeff = snaptr[0]->ncoeff;
@@ -151,7 +153,9 @@ void ComputeSNADAtom::init()
if (strcmp(modify->compute[i]->style,"snad/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute snad/atom");
- #pragma omp parallel
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
{
int tid = omp_get_thread_num();
snaptr[tid]->init();
@@ -198,7 +202,7 @@ void ComputeSNADAtom::compute_peratom()
const int* const ilist = list->ilist;
const int* const numneigh = list->numneigh;
int** const firstneigh = list->firstneigh;
- int *type = atom->type;
+ int * const type = atom->type;
// compute sna derivatives for each atom in group
// use full neighbor list to count atoms less than cutoff
@@ -206,7 +210,9 @@ void ComputeSNADAtom::compute_peratom()
double** const x = atom->x;
const int* const mask = atom->mask;
- #pragma omp parallel for
+#if defined(_OPENMP)
+#pragma omp parallel for default(none)
+#endif
for (int ii = 0; ii < inum; ii++) {
const int tid = omp_get_thread_num();
const int i = ilist[ii];
@@ -324,6 +330,6 @@ double ComputeSNADAtom::memory_usage()
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
bytes += ncoeff*3;
- bytes += snaptr[0]->memory_usage()*omp_get_max_threads();
+ bytes += snaptr[0]->memory_usage()*comm->nthreads;
return bytes;
}
diff --git a/src/SNAP/compute_snav_atom.cpp b/src/SNAP/compute_snav_atom.cpp
index 4b366e649bdfc06e15369ac7415d4fef08fa46b7..76f3d3aa620a7f897ae0fc339656333b41a69311 100644
--- a/src/SNAP/compute_snav_atom.cpp
+++ b/src/SNAP/compute_snav_atom.cpp
@@ -47,7 +47,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
// default values
diagonalstyle = 0;
- rmin0 = 0;
+ rmin0 = 0.0;
switchflag = 1;
// process required arguments
@@ -96,15 +96,17 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else error->all(FLERR,"Illegal compute snav/atom command");
}
- // use_shared_arrays doesn't work with computes and should be disabled.
- int use_shared_arrays = 0;
- int nthreads = omp_get_max_threads();
- snaptr = new SNA*[nthreads];
-#pragma omp parallel shared(rfac0,twojmax,rmin0,switchflag)
+
+ snaptr = new SNA*[comm->nthreads];
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag)
+#endif
{
int tid = omp_get_thread_num();
- snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,use_shared_arrays,
- rmin0,switchflag);
+
+ // always unset use_shared_arrays since it does not work with computes
+ snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,
+ 0 /*use_shared_arrays*/, rmin0,switchflag);
}
ncoeff = snaptr[0]->ncoeff;
@@ -154,7 +156,9 @@ void ComputeSNAVAtom::init()
if (strcmp(modify->compute[i]->style,"snav/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute snav/atom");
- #pragma omp parallel
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
{
int tid = omp_get_thread_num();
snaptr[tid]->init();
@@ -201,14 +205,16 @@ void ComputeSNAVAtom::compute_peratom()
const int* const ilist = list->ilist;
const int* const numneigh = list->numneigh;
int** const firstneigh = list->firstneigh;
- int *type = atom->type;
+ int * const type = atom->type;
// compute sna derivatives for each atom in group
// use full neighbor list to count atoms less than cutoff
double** const x = atom->x;
const int* const mask = atom->mask;
- #pragma omp parallel for
+#if defined(_OPENMP)
+#pragma omp parallel for default(none)
+#endif
for (int ii = 0; ii < inum; ii++) {
const int tid = omp_get_thread_num();
const int i = ilist[ii];
@@ -334,6 +340,6 @@ double ComputeSNAVAtom::memory_usage()
bytes += 3*njmax*sizeof(double);
bytes += njmax*sizeof(int);
bytes += ncoeff*nvirial;
- bytes += snaptr[0]->memory_usage()*omp_get_max_threads();
+ bytes += snaptr[0]->memory_usage()*comm->nthreads;
return bytes;
}
diff --git a/src/SNAP/openmp_snap.h b/src/SNAP/openmp_snap.h
index 28becaea70e9b68461598f190bb5a550197e99d8..60a3138c9c751b3d567841ad404f83b9834d7ebb 100644
--- a/src/SNAP/openmp_snap.h
+++ b/src/SNAP/openmp_snap.h
@@ -1,11 +1,16 @@
-#ifdef _OPENMP
+
+#ifndef LMP_OPENMP_SNAP_H
+#define LMP_OPENMP_SNAP_H
+
+#if defined(_OPENMP)
#include <omp.h>
#else
enum omp_sched_t {omp_sched_static, omp_sched_dynamic, omp_sched_guided, omp_sched_auto};
inline int omp_get_thread_num() { return 0;}
-inline int omp_get_max_threads() { return 1;}
inline int omp_set_num_threads(int num_threads) {return 1;}
/* inline int __sync_fetch_and_add(int* ptr, int value) {int tmp = *ptr; ptr[0]+=value; return tmp;} */
inline void omp_set_schedule(omp_sched_t schedule,int modifier=1) {}
inline int omp_in_parallel() {return 0;}
#endif
+
+#endif
diff --git a/src/SNAP/pair_snap.cpp b/src/SNAP/pair_snap.cpp
index 9060dbbc54e840f1d8fd1a0fcc98ea4d3bb64520..1c56bba07cc35dbda33a419a064b197dc4771759 100644
--- a/src/SNAP/pair_snap.cpp
+++ b/src/SNAP/pair_snap.cpp
@@ -170,7 +170,6 @@ void PairSNAP::compute_regular(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
- int *tag = atom->tag;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
class SNA* snaptr = sna[0];
@@ -377,7 +376,9 @@ void PairSNAP::compute_optimized(int eflag, int vflag)
if (use_shared_arrays)
build_per_atom_arrays();
- #pragma omp parallel shared(time_dynamic,time_guided) firstprivate(numpairs)
+#if defined(_OPENMP)
+#pragma omp parallel shared(eflag,vflag,time_dynamic,time_guided) firstprivate(numpairs) default(none)
+#endif
{
// begin of pragma omp parallel
@@ -419,24 +420,28 @@ void PairSNAP::compute_optimized(int eflag, int vflag)
}
int ielem;
- int i,j,ii,jj,k,inum,jnum,itype,jtype,ninside;
+ int jj,k,inum,jnum,jtype,ninside;
double delx,dely,delz,evdwl,rsq;
double fij[3];
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
- #pragma omp master
+#if defined(_OPENMP)
+#pragma omp master
+#endif
{
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
}
- #pragma omp barrier
+#if defined(_OPENMP)
+#pragma omp barrier
+ { ; }
+#endif
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
- int *tag = atom->tag;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
@@ -463,7 +468,9 @@ void PairSNAP::compute_optimized(int eflag, int vflag)
if (time_dynamic || time_guided)
starttime = MPI_Wtime();
- #pragma omp for schedule(runtime)
+#if defined(_OPENMP)
+#pragma omp for schedule(runtime)
+#endif
for (int iijj = 0; iijj < numpairs; iijj++) {
int i = 0;
if (use_shared_arrays) {
@@ -571,7 +578,9 @@ void PairSNAP::compute_optimized(int eflag, int vflag)
fij[2] += bgb*sna[tid]->dbvec[k-1][2];
}
- #pragma omp critical
+#if defined(_OPENMP)
+#pragma omp critical
+#endif
{
f[i][0] += fij[0];
f[i][1] += fij[1];
@@ -602,7 +611,9 @@ void PairSNAP::compute_optimized(int eflag, int vflag)
for (int k = 1; k <= ncoeff; k++)
evdwl += coeffi[k]*pow(sna[tid]->bvec[k-1],gamma);
- #pragma omp critical
+#if defined(_OPENMP)
+#pragma omp critical
+#endif
ev_tally_full(i,2.0*evdwl,0.0,0.0,delx,dely,delz);
}
@@ -724,17 +735,12 @@ void PairSNAP::load_balance()
for (int i=0; i < list->inum; i++)
if (ilistmask[i]) nlocal++;
int ***grid2proc = comm->grid2proc;
- int* myloc = comm->myloc;
int* procgrid = comm->procgrid;
int nlocal_up,nlocal_down;
MPI_Status status;
MPI_Request request;
- double sendbuf[100];
- double recvbuf[100];
-
- double** x=atom->x;
double sub_mid[3];
for (int dim=0; dim<3; dim++)
sub_mid[dim] = (subhi[dim] + sublo[dim])/2;
@@ -754,7 +760,6 @@ void PairSNAP::load_balance()
int nrecv = ghostinum;
int totalsend = 0;
int nsend = 0;
- int recvatoms[27][4];
int depth = 1;
for (int dx = -depth; dx < depth+1; dx++)
@@ -1163,10 +1168,11 @@ void PairSNAP::build_per_atom_arrays()
clock_gettime(CLOCK_REALTIME,&starttime);
#endif
- #pragma omp parallel for shared(count)
+#if defined(_OPENMP)
+#pragma omp parallel for shared(count) default(none)
+#endif
for (int ii=0; ii < count; ii++) {
int tid = omp_get_thread_num();
- int i=i_list[ii];
set_sna_to_shared(tid,ii);
//sna[tid]->compute_ui(i_ninside[ii]);
#ifdef TIMING_INFO
@@ -1185,7 +1191,6 @@ void PairSNAP::build_per_atom_arrays()
#endif
for (int ii=0; ii < count; ii++) {
int tid = 0;//omp_get_thread_num();
- int i=i_list[ii];
set_sna_to_shared(tid,ii);
sna[tid]->compute_zi_omp(MAX(int(nthreads/count),1));
}
@@ -1237,7 +1242,12 @@ void PairSNAP::settings(int narg, char **arg)
" Too few arguments.");
if (strcmp(arg[i],"nthreads")==0) {
nthreads=atoi(arg[++i]);
+#if defined(LMP_USER_OMP)
+ error->all(FLERR,"Please set number of threads via package omp command");
+#else
omp_set_num_threads(nthreads);
+ comm->nthreads=nthreads;
+#endif
continue;
}
if (strcmp(arg[i],"optimized")==0) {
@@ -1295,12 +1305,13 @@ void PairSNAP::settings(int narg, char **arg)
}
if (nthreads < 0)
- nthreads = omp_get_max_threads();
+ nthreads = comm->nthreads;
- if (use_shared_arrays < 0)
- if (nthreads > 1 && atom->nlocal <= 2*nthreads)
+ if (use_shared_arrays < 0) {
+ if (nthreads > 1 && atom->nlocal <= 2*nthreads)
use_shared_arrays = 1;
else use_shared_arrays = 0;
+ }
// check if running non-optimized code with
// optimization flags set
@@ -1405,7 +1416,9 @@ void PairSNAP::coeff(int narg, char **arg)
// allocate memory for per OpenMP thread data which
// is wrapped into the sna class
-#pragma omp parallel shared(rfac0,twojmax,diagonalstyle,rmin0,switchflag)
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
{
int tid = omp_get_thread_num();
sna[tid] = new SNA(lmp,rfac0,twojmax,
@@ -1443,7 +1456,9 @@ void PairSNAP::init_style()
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
- #pragma omp parallel
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
{
int tid = omp_get_thread_num();
sna[tid]->init();
diff --git a/src/SNAP/sna.cpp b/src/SNAP/sna.cpp
index 970453d8d522d4c93b260dc332748c3ae69b00d3..10e9dac4afb8f84034c5347302eaf11e55999975 100644
--- a/src/SNAP/sna.cpp
+++ b/src/SNAP/sna.cpp
@@ -160,7 +160,6 @@ void SNA::build_indexlist()
{
if(diagonalstyle == 0) {
int idxj_count = 0;
- int ma2, mb2;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
@@ -186,7 +185,6 @@ void SNA::build_indexlist()
if(diagonalstyle == 1) {
int idxj_count = 0;
- int ma2, mb2;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j = 0; j <= MIN(twojmax, 2 * j1); j += 2) {
@@ -211,7 +209,6 @@ void SNA::build_indexlist()
if(diagonalstyle == 2) {
int idxj_count = 0;
- int ma2, mb2;
for(int j1 = 0; j1 <= twojmax; j1++) {
idxj_count++;
@@ -234,7 +231,6 @@ void SNA::build_indexlist()
if(diagonalstyle == 3) {
int idxj_count = 0;
- int ma2, mb2;
for(int j1 = 0; j1 <= twojmax; j1++)
for(int j2 = 0; j2 <= j1; j2++)
@@ -358,7 +354,9 @@ void SNA::compute_ui_omp(int jnum, int sub_threads)
z0 = r / tan(theta0);
omp_set_num_threads(sub_threads);
- #pragma omp parallel shared(x,y,z,z0,r,sub_threads)
+#if defined(_OPENMP)
+#pragma omp parallel shared(x,y,z,z0,r,sub_threads) default(none)
+#endif
{
compute_uarray_omp(x, y, z, z0, r, sub_threads);
}
@@ -464,7 +462,9 @@ void SNA::compute_zi_omp(int sub_threads)
// compute_dbidrj() requires full j1/j2/j chunk of z elements
// use zarray j1/j2 symmetry
- #pragma omp parallel for schedule(auto)
+#if defined(_OPENMP)
+#pragma omp parallel for schedule(auto) default(none)
+#endif
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) {
@@ -670,7 +670,7 @@ void SNA::compute_dbidrj_nonsymm()
double* dbdr;
double* dudr_r, *dudr_i;
double sumb1_r[3], sumb1_i[3], dzdr_r[3], dzdr_i[3];
- int ma2, mb2;
+ int ma2;
#ifdef TIMING_INFO
clock_gettime(CLOCK_REALTIME, &starttime);
@@ -1124,7 +1124,9 @@ void SNA::add_uarraytot_omp(double r, double wj, double rcut)
sfac *= wj;
- #pragma omp for
+#if defined(_OPENMP)
+#pragma omp for
+#endif
for (int j = 0; j <= twojmax; j++)
for (int ma = 0; ma <= j; ma++)
for (int mb = 0; mb <= j; mb++) {
@@ -1232,7 +1234,9 @@ void SNA::compute_uarray_omp(double x, double y, double z,
uarray_i[0][0][0] = 0.0;
for (int j = 1; j <= twojmax; j++) {
- #pragma omp for
+#if defined(_OPENMP)
+#pragma omp for
+#endif
for (int mb = 0; mb < j; mb++) {
uarray_r[j][0][mb] = 0.0;
uarray_i[j][0][mb] = 0.0;
@@ -1264,7 +1268,9 @@ void SNA::compute_uarray_omp(double x, double y, double z,
uarray_r[j][0][mb] = 0.0;
uarray_i[j][0][mb] = 0.0;
- #pragma omp for
+#if defined(_OPENMP)
+#pragma omp for
+#endif
for (int ma = 0; ma < j; ma++) {
rootpq = rootpqarray[j - ma][mb];
uarray_r[j][ma][mb] +=
@@ -1917,6 +1923,7 @@ double SNA::compute_sfac(double r, double rcut)
return 0.5 * (cos((r - rmin0) * rcutfac) + 1.0);
}
}
+ return 0.0;
}
/* ---------------------------------------------------------------------- */
@@ -1932,151 +1939,6 @@ double SNA::compute_dsfac(double r, double rcut)
return -0.5 * sin((r - rmin0) * rcutfac) * rcutfac;
}
}
-}
-
-/* ----------------------------------------------------------------------
-
- This compares SNA results with Mathematica output for BCC crystal
-
- To build this test with LAMMPS,
- include sna.h in lammps.cpp
- and add the line:
-
- SNA *sna = new SNA(this);
- sna->test();
-
- to the bottom of LAMMPS::create()
-
- To run the test, just run LAMMPS with no input script
-------------------------------------------------------------------------- */
-
-void SNA::test()
-{
-
- // rotated nearest-neighbor bond vectors for BCC crystal
-
- int natoms = 8;
- double rijpolar[] = {
- 1.959, 0.9730636300319065, 0.825999311173714, 1.959,
- 0.8980297487727392, 2.423647971942037, 1.959, 0.93989780135463,
- -2.2344890150374734, 1.959, 1.0127992115193267,
- -0.6943594028137655, 1.959, 2.201694852235163, 0.9071036385523198,
- 1.959, 2.1287934420704664, 2.4472332507760277, 1.959,
- 2.1685290235578867, -2.3155933424160793, 1.959,
- 2.2435629048170536, -0.7179446816477564
- };
-
- // results from Mathematica
- // the quantities bcbrt0 are the cube roots of the
- // bispectrum coefficients, as plotted in Bartok's
- // thesis, Fig 2.5.
-
- twojmax = 4;
- double bcbrt0[] = {
- 9., 4.199196953160132, 0.8023880044758301, 7.111020498078093,
- 6.698177303873683, 4.199196953160132, 3.3329048304340256,
- -1.1951270630399309, -1.0440396551493507, 1.525005828080911,
- 1.3855597471974264, -2.289695002886823, -2.125564551555691,
- 0.8023880044758301, -1.0440396551493507, 1.5250058280809107,
- 0.5563447693070245, -0.19949645249634843, 0.32967616972020086,
- 1.2103979278766366, -1.7680017087520263, 0.2780596505135809,
- 1.4288971890382558, 7.111020498078093, 1.3855597471974261,
- -2.289695002886823, 1.2103979278766366, -1.7680017087520263,
- 4.479662205881257, -1.606336156567859, 2.654537184765498,
- -1.687061702582909, 2.464254031043488, 6.698177303873683,
- -2.125564551555691, 0.2780596505135809, 1.428897189038256,
- -1.687061702582909, 2.464254031043488, 3.917117850044296,
- 1.2051785645927136, -3.627660818200106
- };
-
- // assume fixed theta0, all neighbors are equidistant
-
- double theta0 = rijpolar[0];
- diagonalstyle = 0;
- ncoeff = compute_ncoeff();
- create_twojmax_arrays();
-
- // theta0 = rfac0*MY_PI*r
- // r = 1 by construction
-
- double rcut = 0.0;
- rfac0 = theta0 / MY_PI;
- double rscale0 = MY_PI * rfac0;
-
- printf("SNA::test() validating bispectrum coefficients "
- "against Mathamatica calculation for BCC crystal\n");
- printf("\n ncoeff = %d\n rcut = %g\n rfac0 = %g\n"
- "twojmax = %d\n diagonalstyle = %d\n",
- ncoeff, rcut, rfac0, twojmax, diagonalstyle);
-
- //double *bvec = (double *) malloc(ncoeff*sizeof(double));
- memory->destroy(bvec);
- memory->destroy(dbvec);
- memory->create(bvec, ncoeff, "pair:bvec");
- memory->create(dbvec, ncoeff, 3, "pair:dbvec");
- double* data = (double*) malloc(3 * natoms * sizeof(double));
- grow_rij(natoms);
-
- double theta, phi;
- int n = 0;
-
- for (int iatom = 0; iatom < natoms; iatom++) {
- rij[iatom] = &data[n];
- n += 3;
- }
-
- n = 0;
-
- for (int iatom = 0; iatom < natoms; iatom++) {
- theta = rijpolar[n + 1];
- phi = rijpolar[n + 2];
- rij[iatom][0] = sin(theta) * cos(phi);
- rij[iatom][1] = sin(theta) * sin(phi);
- rij[iatom][2] = cos(theta);
- n += 3;
- }
-
- init();
-
- // compute Ui and Zi for atom I
-
- compute_ui(natoms);
- compute_zi();
-
- compute_bi();
- copy_bi2bvec();
-
- double bcbrt, pcterr;
-
- printf("\nBCC Bispectrum Coefficients (Bartok, Figure 2.5) \n");
- printf("2x[J1,J2,J] Mathematica SNA Delta (%%) \n");
-
- // use explicit nested loop because
- // j2 > j1 coefficients skipped in bvec.
-
- int k = 0;
- int kk = 0;
-
- for (int j1 = 0; j1 <= twojmax; j1++)
- for (int j2 = 0; j2 <= twojmax; j2++)
- for (int j = abs(j1 - j2);
- j <= MIN(twojmax, j1 + j2); j += 2) {
- if(j2 <= j1) {
- bcbrt = cbrt(bvec[kk]);
- pcterr = 100.0 * (bcbrt / bcbrt0[k] - 1);
- kk++;
- } else {
- bcbrt = 0.0;
- pcterr = 0.0;
- }
-
- printf("%3d %1d%1d%1d %g %g %g\n", k + 1, j1, j2, j, bcbrt0[k],
- bcbrt, pcterr);
- k++;
- }
-
- MPI_Finalize();
- exit(1);
-
+ return 0.0;
}
diff --git a/src/SNAP/sna.h b/src/SNAP/sna.h
index e391c0e7c146e8c3609eb1371bb8b1c4f3832eb4..1c38bcb920fb1c5a71c609dd1a776e3f6943fae7 100644
--- a/src/SNAP/sna.h
+++ b/src/SNAP/sna.h
@@ -28,7 +28,6 @@ public:
~SNA();
void build_indexlist();
void init();
- void test();
double memory_usage();
int ncoeff;