Precompute some sqrt factors

src/KOKKOS/compute_orientorder_atom_kokkos.cpp

@@ -62,6 +62,20 @@ ComputeOrientOrderAtomKokkos<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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 = d_qnormfac(il)/d_qnarray(i,il);
-        const double qnfac = qnormfac/d_qnarray(i,il);
+        d_qnarray(i,jj++) = d_qnarray(i,jptr+il) * (qnfac*qnfac*qnfac) * d_qnormfac2(il);
-        d_qnarray(i,jj++) = d_qnarray(i,jptr+il) * (qnfac*qnfac*qnfac) * sqrt(2.0*l+1.0);
       }
     }
   }
@@ -586,8 +597,7 @@ void ComputeOrientOrderAtomKokkos<DeviceType>::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 = d_qnormfac(il)/d_qnarray(i,il);
-      const double qnfac = qnormfac/d_qnarray(i,il);
       for (int m = -l; m < 0; m++) {
         // Computed only qnm for m>=0.
         // qnm[-m] = (-1)^m * conjg(qnm[m])

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<class TagStyle>
   void check_team_size_for(int, int&, int);

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[il]/qn[il];
-        double qnfac = qnormfac/qn[il];
+        qn[jj++] = qn[jptr+il] * (qnfac*qnfac*qnfac) * qnormfac2[il];
-        qn[jj++] = qn[jptr+il] * (qnfac*qnfac*qnfac) * sqrt(2*l+1);
       }
     }
   }
@@ -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[il]/qn[il];
-      double qnfac = qnormfac/qn[il];
       for (int m = -l; m < 0; m++) {
         // Computed only qnm for m>=0.
         // qnm[-m] = (-1)^m * conjg(qnm[m])

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;