diff --git a/src/KOKKOS/compute_orientorder_atom_kokkos.cpp b/src/KOKKOS/compute_orientorder_atom_kokkos.cpp index 5cc7ef104b..6edb500870 100644 --- a/src/KOKKOS/compute_orientorder_atom_kokkos.cpp +++ b/src/KOKKOS/compute_orientorder_atom_kokkos.cpp @@ -62,6 +62,20 @@ ComputeOrientOrderAtomKokkos::ComputeOrientOrderAtomKokkos(LAMMPS *l datamask_modify = EMPTY_MASK; host_flag = (execution_space == Host); + + d_qnormfac = t_sna_1d("orientorder/atom:qnormfac",nqlist); + d_qnormfac2 = t_sna_1d("orientorder/atom:qnormfac2",nqlist); + + auto h_qnormfac = Kokkos::create_mirror_view(d_qnormfac); + auto h_qnormfac2 = Kokkos::create_mirror_view(d_qnormfac2); + + for (int il = 0; il < nqlist; il++) { + h_qnormfac[il] = qnormfac[il]; + h_qnormfac2[il] = qnormfac2[il]; + } + + Kokkos::deep_copy(d_qnormfac,h_qnormfac); + Kokkos::deep_copy(d_qnormfac2,h_qnormfac2); } /* ---------------------------------------------------------------------- */ @@ -518,11 +532,10 @@ void ComputeOrientOrderAtomKokkos::calc_boop2(int ncount, int ii) co int jj = 0; for (int il = 0; il < nqlist; il++) { int l = d_qlist[il]; - double qnormfac = sqrt(MY_4PI/(2*l+1)); double qm_sum = d_qnm(ii,il,0).re*d_qnm(ii,il,0).re; for (int m = 1; m < l+1; m++) qm_sum += 2.0*(d_qnm(ii,il,m).re*d_qnm(ii,il,m).re + d_qnm(ii,il,m).im*d_qnm(ii,il,m).im); - d_qnarray(i,jj++) = qnormfac * sqrt(qm_sum); + d_qnarray(i,jj++) = d_qnormfac(il) * sqrt(qm_sum); } // calculate W_l @@ -553,7 +566,7 @@ void ComputeOrientOrderAtomKokkos::calc_boop2(int ncount, int ii) co wlsum += Q1Q2Q3*c; } } - d_qnarray(i,jj++) = wlsum/sqrt(2.0*l+1.0); + d_qnarray(i,jj++) = wlsum/d_qnormfac2(il); nterms++; } } @@ -564,13 +577,11 @@ void ComputeOrientOrderAtomKokkos::calc_boop2(int ncount, int ii) co const int jptr = jj-nterms; if (!wlflag) jj = jptr; for (int il = 0; il < nqlist; il++) { - int l = d_qlist[il]; if (d_qnarray(i,il) < QEPSILON) d_qnarray(i,jj++) = 0.0; else { - const double qnormfac = sqrt(MY_4PI/(2*l+1)); - const double qnfac = qnormfac/d_qnarray(i,il); - d_qnarray(i,jj++) = d_qnarray(i,jptr+il) * (qnfac*qnfac*qnfac) * sqrt(2.0*l+1.0); + const double qnfac = d_qnormfac(il)/d_qnarray(i,il); + d_qnarray(i,jj++) = d_qnarray(i,jptr+il) * (qnfac*qnfac*qnfac) * d_qnormfac2(il); } } } @@ -586,8 +597,7 @@ void ComputeOrientOrderAtomKokkos::calc_boop2(int ncount, int ii) co d_qnarray(i,jj++) = 0.0; } else { - const double qnormfac = sqrt(MY_4PI/(2*l+1)); - const double qnfac = qnormfac/d_qnarray(i,il); + const double qnfac = d_qnormfac(il)/d_qnarray(i,il); for (int m = -l; m < 0; m++) { // Computed only qnm for m>=0. // qnm[-m] = (-1)^m * conjg(qnm[m]) diff --git a/src/KOKKOS/compute_orientorder_atom_kokkos.h b/src/KOKKOS/compute_orientorder_atom_kokkos.h index 8afa1dcd40..2d75bd5ab8 100644 --- a/src/KOKKOS/compute_orientorder_atom_kokkos.h +++ b/src/KOKKOS/compute_orientorder_atom_kokkos.h @@ -71,6 +71,7 @@ class ComputeOrientOrderAtomKokkos : public ComputeOrientOrderAtom { void init(); void compute_peratom(); t_sna_1i d_qlist; + t_sna_1d d_qnormfac,d_qnormfac2; template void check_team_size_for(int, int&, int); diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp index fb84bbce06..34357ae850 100644 --- a/src/compute_orientorder_atom.cpp +++ b/src/compute_orientorder_atom.cpp @@ -57,7 +57,8 @@ using namespace MathSpecial; ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), qlist(nullptr), distsq(nullptr), nearest(nullptr), rlist(nullptr), - qnarray(nullptr), qnm_r(nullptr), qnm_i(nullptr), w3jlist(nullptr) + qnarray(nullptr), qnm_r(nullptr), qnm_i(nullptr), w3jlist(nullptr), + qnormfac(nullptr),qnormfac2(nullptr) { if (narg < 3 ) error->all(FLERR,"Illegal compute orientorder/atom command"); @@ -154,9 +155,17 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg nmax = 0; maxneigh = 0; + + memory->create(qnormfac,nqlist,"orientorder/atom:qnormfac"); + memory->create(qnormfac2,nqlist,"orientorder/atom:qnormfac2"); + for (int il = 0; il < nqlist; il++) { + int l = qlist[il]; + qnormfac[il] = sqrt(MY_4PI/(2*l+1)); + qnormfac2[il] = sqrt(2*l+1); + } } -/* ---------------------------------------------------------------------- */ +/* --------------------------------------------------------------------- */ ComputeOrientOrderAtom::~ComputeOrientOrderAtom() { @@ -167,6 +176,8 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom() memory->destroy(rlist); memory->destroy(nearest); memory->destroy(qlist); + memory->destroy(qnormfac); + memory->destroy(qnormfac2); memory->destroy(qnm_r); memory->destroy(qnm_i); memory->destroy(w3jlist); @@ -497,11 +508,10 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, int jj = 0; for (int il = 0; il < nqlist; il++) { int l = qlist[il]; - double qnormfac = sqrt(MY_4PI/(2*l+1)); double qm_sum = qnm_r[il][0]*qnm_r[il][0]; for (int m = 1; m < l+1; m++) qm_sum += 2.0*(qnm_r[il][m]*qnm_r[il][m] + qnm_i[il][m]*qnm_i[il][m]); - qn[jj++] = qnormfac * sqrt(qm_sum); + qn[jj++] = qnormfac[il] * sqrt(qm_sum); } // calculate W_l @@ -531,7 +541,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, wlsum += Q1Q2Q3*c; } } - qn[jj++] = wlsum/sqrt(2*l+1); + qn[jj++] = wlsum/qnormfac2[il]; nterms++; } } @@ -542,13 +552,11 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, const int jptr = jj-nterms; if (!wlflag) jj = jptr; for (int il = 0; il < nqlist; il++) { - int l = qlist[il]; if (qn[il] < QEPSILON) qn[jj++] = 0.0; else { - double qnormfac = sqrt(MY_4PI/(2*l+1)); - double qnfac = qnormfac/qn[il]; - qn[jj++] = qn[jptr+il] * (qnfac*qnfac*qnfac) * sqrt(2*l+1); + double qnfac = qnormfac[il]/qn[il]; + qn[jj++] = qn[jptr+il] * (qnfac*qnfac*qnfac) * qnormfac2[il]; } } } @@ -564,8 +572,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, qn[jj++] = 0.0; } else { - double qnormfac = sqrt(MY_4PI/(2*l+1)); - double qnfac = qnormfac/qn[il]; + double qnfac = qnormfac[il]/qn[il]; for (int m = -l; m < 0; m++) { // Computed only qnm for m>=0. // qnm[-m] = (-1)^m * conjg(qnm[m]) diff --git a/src/compute_orientorder_atom.h b/src/compute_orientorder_atom.h index bc45333dd6..3941575c1f 100644 --- a/src/compute_orientorder_atom.h +++ b/src/compute_orientorder_atom.h @@ -36,6 +36,7 @@ class ComputeOrientOrderAtom : public Compute { int iqlcomp, qlcomp, qlcompflag, wlflag, wlhatflag; int *qlist; int nqlist; + double *qnormfac,*qnormfac2; protected: int nmax, maxneigh, ncol, nnn;