Switched the sign of spherical harmonics for m odd

doc/src/compute_orientorder_atom.rst

@@ -48,14 +48,17 @@ For each atom, :math:`Q_l` is a real number defined as follows:
    \bar{Y}_{lm} = & \frac{1}{nnn}\sum_{j = 1}^{nnn} Y_{lm}( \theta( {\bf r}_{ij} ), \phi( {\bf r}_{ij} ) ) \\
    Q_l = & \sqrt{\frac{4 \pi}{2 l + 1} \sum_{m = -l}^{m = l} \bar{Y}_{lm} \bar{Y}^*_{lm}}
 
-The first equation defines the spherical harmonic order parameters.
+The first equation defines the local order parameters as averages
+of the spherical harmonics :math:`Y_{lm}` for each neighbor.
 These are complex number components of the 3D analog of the 2D order
 parameter :math:`q_n`, which is implemented as LAMMPS compute
 :doc:`hexorder/atom <compute_hexorder_atom>`.
 The summation is over the *nnn* nearest
 neighbors of the central atom.
-The angles theta and phi are the standard spherical polar angles
+The angles :math:`theta` and :math:`phi` are the standard spherical polar angles
 defining the direction of the bond vector :math:`r_{ij}`.
+The phase and sign of :math:`Y_{lm}` follow the standard conventions, 
+so that :math:`{\rm sign}(Y_{ll}(0,0)) = (-1)^l`.   
 The second equation defines :math:`Q_l`, which is a
 rotationally invariant non-negative amplitude obtained by summing
 over all the components of degree *l*\ .
@@ -98,8 +101,8 @@ structures are given in Table I of :ref:`Steinhardt <Steinhardt>`, and these
 can be reproduced with this keyword.
 
 The optional keyword *components* will output the components of the
-normalized complex vector :math:`\bar{Y}_{lm}` of degree *ldegree*\ , which must be
+*normalized* complex vector :math:`\hat{Y}_{lm} = \bar{Y}_{lm}/|\bar{Y}_{lm}|` of degree *ldegree*\,
-explicitly included in the keyword *degrees*\ . This option can be used
+which must be included in the list of order parameters to be computed. This option can be used
 in conjunction with :doc:`compute coord_atom <compute_coord_atom>` to
 calculate the ten Wolde's criterion to identify crystal-like
 particles, as discussed in :ref:`ten Wolde <tenWolde2>`.
@@ -141,11 +144,15 @@ If the keyword *wl/hat* is set to yes, then the :math:`\hat{W}_l`
 values for each atom will be added to the output array, which are real numbers.
 
 If the keyword *components* is set, then the real and imaginary parts
-of each component of (normalized) :math:`\bar{Y}_{lm}` will be added to the
+of each component of *normalized* :math:`\hat{Y}_{lm}` will be added to the
-output array in the following order: :math:`Re(\bar{Y}_{-m}) Im(\bar{Y}_{-m})
+output array in the following order: :math:`{\rm Re}(\hat{Y}_{-m}), {\rm Im}(\hat{Y}_{-m}), 
-Re(\bar{Y}_{-m+1}) Im(\bar{Y}_{-m+1}) ... Re(\bar{Y}_m) Im(\bar{Y}_m)`.  This
+{\rm Re}(\hat{Y}_{-m+1}), {\rm Im}(\hat{Y}_{-m+1}), \dots , {\rm Re}(\hat{Y}_m), {\rm Im}(\hat{Y}_m)`.  
-way, the per-atom array will have a total of *nlvalues*\ +2\*(2\ *l*\ +1)
+
-columns.
+In summary, the per-atom array will contain *nlvalues* columns, followed by
+an additional *nlvalues* columns if *wl* is set to yes, followed by
+an additional *nlvalues* columns if *wl/hat* is set to yes, followed
+by an additional 2\*(2\* *ldegree*\ +1) columns if the *components* 
+keyword is set.
 
 These values can be accessed by any command that uses per-atom values
 from a compute as input.  See the :doc:`Howto output <Howto_output>` doc

src/compute_orientorder_atom.cpp

@@ -456,21 +456,26 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
     for (int il = 0; il < nqlist; il++) {
       int l = qlist[il];
 
+      // calculate spherical harmonics
+      // Ylm, -l <= m <= l
+      // sign convention: sign(Yll(0,0)) = (-1)^l
+
       qnm_r[il][l] += polar_prefactor(l, 0, costheta);
       double expphim_r = expphi_r;
       double expphim_i = expphi_i;
       for(int m = 1; m <= +l; m++) {
+
         double prefactor = polar_prefactor(l, m, costheta);
-        double c_r = prefactor * expphim_r;
+        double ylm_r = prefactor * expphim_r;
-        double c_i = prefactor * expphim_i;
+        double ylm_i = prefactor * expphim_i;
-        qnm_r[il][m+l] += c_r;
+        qnm_r[il][m+l] += ylm_r;
-        qnm_i[il][m+l] += c_i;
+        qnm_i[il][m+l] += ylm_i;
         if(m & 1) {
-          qnm_r[il][-m+l] -= c_r;
+          qnm_r[il][-m+l] -= ylm_r;
-          qnm_i[il][-m+l] += c_i;
+          qnm_i[il][-m+l] += ylm_i;
         } else {
-          qnm_r[il][-m+l] += c_r;
+          qnm_r[il][-m+l] += ylm_r;
-          qnm_i[il][-m+l] -= c_i;
+          qnm_i[il][-m+l] -= ylm_i;
         }
         double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i;
         double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r;
@@ -505,19 +510,6 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
     qn[jj++] = qnormfac * sqrt(qm_sum);
   }
 
-  // TODO:
-  // 1. [done]Need to allocate extra memory in qnarray[] for this option
-  // 2. [done]Need to add keyword option
-  // 3. [done]Need to calculate Clebsch-Gordan/Wigner 3j coefficients
-  //     (Can try getting them from boop.py first)
-  // 5. [done]Compare to bcc values in /Users/athomps/netapp/codes/MatMiner/matminer/matminer/featurizers/boop.py
-  // 6. [done]I get the right answer for W_l, but need to make sure that factor of 1/sqrt(l+1) is right for cglist
-  // 7. Add documentation
-  // 8. [done] run valgrind
-  // 9. [done] Add Wlhat
-  // 10. Update memory_usage()
-  // 11. Add exact FCC values for W_4, W_4_hat
-
   // calculate W_l
 
   if (wlflag) {
@@ -554,7 +546,6 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
           idxcg_count++;
         }
       }
-  //      Whats = [w/(q/np.sqrt(np.pi * 4 / (2 * l + 1)))**3 if abs(q) > 1.0e-6 else 0.0 for l,q,w in zip(range(1,max_l+1),Qs,Ws)]
       if (qn[il] < QEPSILON)
         qn[jj++] = 0.0;
       else {
@@ -565,7 +556,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
     }
   }
 
-  // Calculate components of Q_l, for l=qlcomp
+  // Calculate components of Q_l/|Q_l|, for l=qlcomp
 
   if (qlcompflag) {
     int il = iqlcomp;
@@ -619,6 +610,7 @@ double ComputeOrientOrderAtom::polar_prefactor(int l, int m, double costheta)
 
 /* ----------------------------------------------------------------------
    associated legendre polynomial
+   sign convention: P(l,l) = (2l-1)!!(-sqrt(1-x^2))^l
 ------------------------------------------------------------------------- */
 
 double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x)
@@ -628,9 +620,9 @@ double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x)
   double p(1.0), pm1(0.0), pm2(0.0);
 
   if (m != 0) {
-    const double sqx = sqrt(1.0-x*x);
+    const double msqx = -sqrt(1.0-x*x);
     for (int i=1; i < m+1; ++i)
-      p *= static_cast<double>(2*i-1) * sqx;
+      p *= static_cast<double>(2*i-1) * msqx;
   }
 
   for (int i=m+1; i < l+1; ++i) {