Merge pull request #2171 from athomps/orientorder-components-parity

   \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*\ .

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

explicitly included in the keyword *degrees*\ . This option can be used

*normalized* complex vector :math:`\hat{Y}_{lm} = \bar{Y}_{lm}/|\bar{Y}_{lm}|` of degree *ldegree*\,

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>`.

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

output array in the following order: :math:`Re(\bar{Y}_{-m}) Im(\bar{Y}_{-m})

Re(\bar{Y}_{-m+1}) Im(\bar{Y}_{-m+1}) ... Re(\bar{Y}_m) Im(\bar{Y}_m)`.  This

way, the per-atom array will have a total of *nlvalues*\ +2\*(2\ *l*\ +1)

columns.

of each component of *normalized* :math:`\hat{Y}_{lm}` will be added to the

output array in the following order: :math:`{\rm Re}(\hat{Y}_{-m}), {\rm Im}(\hat{Y}_{-m}), 

{\rm Re}(\hat{Y}_{-m+1}), {\rm Im}(\hat{Y}_{-m+1}), \dots , {\rm Re}(\hat{Y}_m), {\rm Im}(\hat{Y}_m)`.  

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

  for (int il = 0; il < nqlist; il++) {

    const int l = d_qlist[il];

    // calculate spherical harmonics

    // Ylm, -l <= m <= l

    // sign convention: sign(Yll(0,0)) = (-1)^l

    //d_qnm(ii,il,l).re += polar_prefactor(l, 0, costheta);

    const double polar_pf = polar_prefactor(l, 0, costheta);

    Kokkos::atomic_add(&(d_qnm(ii,il,l).re), polar_pf);

    SNAcomplex expphim = {expphi.re,expphi.im};

    for(int m = 1; m <= +l; m++) {

      const double prefactor = polar_prefactor(l, m, costheta);

      SNAcomplex c = {prefactor * expphim.re, prefactor * expphim.im};

      //d_qnm(ii,il,m+l).re += c.re;

      //d_qnm(ii,il,m+l).im += c.im;

      Kokkos::atomic_add(&(d_qnm(ii,il,m+l).re), c.re);

      Kokkos::atomic_add(&(d_qnm(ii,il,m+l).im), c.im);

      SNAcomplex ylm = {prefactor * expphim.re, prefactor * expphim.im};

      //d_qnm(ii,il,m+l).re += ylm.re;

      //d_qnm(ii,il,m+l).im += ylm.im;

      Kokkos::atomic_add(&(d_qnm(ii,il,m+l).re), ylm.re);

      Kokkos::atomic_add(&(d_qnm(ii,il,m+l).im), ylm.im);

      if(m & 1) {

        //d_qnm(ii,il,-m+l).re -= c.re;

        //d_qnm(ii,il,-m+l).im += c.im;

        Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).re), -c.re);

        Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).im), c.im);

        //d_qnm(ii,il,-m+l).re -= ylm.re;

        //d_qnm(ii,il,-m+l).im += ylm.im;

        Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).re), -ylm.re);

        Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).im), ylm.im);

      } else {

        //d_qnm(ii,il,-m+l).re += c.re;

        //d_qnm(ii,il,-m+l).im -= c.im;

        Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).re), c.re);

        Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).im), -c.im);

        //d_qnm(ii,il,-m+l).re += ylm.re;

        //d_qnm(ii,il,-m+l).im -= ylm.im;

        Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).re), ylm.re);

        Kokkos::atomic_add(&(d_qnm(ii,il,-m+l).im), -ylm.im);

      }

      SNAcomplex tmp;

      tmp.re = expphim.re*expphi.re - expphim.im*expphi.im;

          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 (d_qnarray(i,il) < QEPSILON)

        d_qnarray(i,jj++) = 0.0;

      else {

    }

  }

  // Calculate components of Q_l, for l=qlcomp

  // Calculate components of Q_l/|Q_l|, for l=qlcomp

  if (qlcompflag) {

    const int il = iqlcomp;

/* ----------------------------------------------------------------------

   associated legendre polynomial

   sign convention: P(l,l) = (2l-1)!!(-sqrt(1-x^2))^l

------------------------------------------------------------------------- */

template<class DeviceType>

  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) {

    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 c_i = prefactor * expphim_i;

        qnm_r[il][m+l] += c_r;

        qnm_i[il][m+l] += c_i;

        double ylm_r = prefactor * expphim_r;

        double ylm_i = prefactor * expphim_i;

        qnm_r[il][m+l] += ylm_r;

        qnm_i[il][m+l] += ylm_i;

        if(m & 1) {

          qnm_r[il][-m+l] -= c_r;

          qnm_i[il][-m+l] += c_i;

          qnm_r[il][-m+l] -= ylm_r;

          qnm_i[il][-m+l] += ylm_i;

        } else {

          qnm_r[il][-m+l] += c_r;

          qnm_i[il][-m+l] -= c_i;

          qnm_r[il][-m+l] += ylm_r;

          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;

    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) {

          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 {

    }

  }

  // Calculate components of Q_l, for l=qlcomp

  // Calculate components of Q_l/|Q_l|, for l=qlcomp

  if (qlcompflag) {

    int il = iqlcomp;

/* ----------------------------------------------------------------------

   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)

  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) {