diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp
index 258aad64b8..0ffab48332 100644
--- a/src/compute_orientorder_atom.cpp
+++ b/src/compute_orientorder_atom.cpp
@@ -56,6 +56,7 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
 
   nnn = 12;
   cutsq = 0.0;
+  wlflag = 0;
   qlcompflag = 0;
 
   // specify which orders to request
@@ -103,6 +104,13 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
         if (qlist[il] > qmax) qmax = qlist[il];
       }
       iarg += nqlist;
+    } else if (strcmp(arg[iarg],"wlparams") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute orientorder/atom command");
+      if (strcmp(arg[iarg+1],"yes") == 0) wlflag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) wlflag = 0;
+      else error->all(FLERR,"Illegal compute orientorder/atom command");
+      iarg += 2;
     } else if (strcmp(arg[iarg],"components") == 0) {
       qlcompflag = 1;
       if (iarg+2 > narg)
@@ -130,8 +138,9 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
     } else error->all(FLERR,"Illegal compute orientorder/atom command");
   }
 
-  if (qlcompflag) ncol = nqlist + 2*(2*qlcomp+1);
-  else ncol = nqlist;
+  ncol = nqlist;
+  if (wlflag) ncol += nqlist;
+  if (qlcompflag) ncol += nqlist + 2*(2*qlcomp+1);
 
   peratom_flag = 1;
   size_peratom_cols = ncol;
@@ -151,7 +160,7 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
   memory->destroy(qlist);
   memory->destroy(qnm_r);
   memory->destroy(qnm_i);
-
+  memory->destroy(cglist);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -183,6 +192,8 @@ void ComputeOrientOrderAtom::init()
     if (strcmp(modify->compute[i]->style,"orientorder/atom") == 0) count++;
   if (count > 1 && comm->me == 0)
     error->warning(FLERR,"More than one compute orientorder/atom");
+
+  if (wlflag) init_clebsch_gordan();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -274,8 +285,8 @@ void ComputeOrientOrderAtom::compute_peratom()
       // if not nnn neighbors, order parameter = 0;
 
       if ((ncount == 0) || (ncount < nnn)) {
-        for (int il = 0; il < nqlist; il++)
-          qn[il] = 0.0;
+        for (int jj = 0; jj < ncol; jj++)
+          qn[jj] = 0.0;
         continue;
       }
 
@@ -461,52 +472,71 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
     }
   }
 
+  // convert sums to averages
+
+  double facn = 1.0 / ncount;
+  for (int il = 0; il < nqlist; il++) {
+    int l = qlist[il];
+    for(int m = 0; m < 2*l+1; m++) {
+      qnm_r[il][m] *= facn;
+      qnm_i[il][m] *= facn;
+    }
+  }
+
   // calculate Q_l
 
-  double fac = sqrt(MY_4PI) / ncount;
+  double facpi = sqrt(MY_4PI);
   double normfac = 0.0;
-  int jcount = 0;
+  int jj = 0;
   for (int il = 0; il < nqlist; il++) {
     int l = qlist[il];
     double qm_sum = 0.0;
     for(int m = 0; m < 2*l+1; m++)
       qm_sum += qnm_r[il][m]*qnm_r[il][m] + qnm_i[il][m]*qnm_i[il][m];
 
-    qn[jcount++] = fac * sqrt(qm_sum / (2*l+1));
+    for(int m = 0; m < 2*l+1; m++)
+      printf("Q_l = %d %g %g\n",l, qnm_r[il][m], qnm_i[il][m]);
+
+    qn[jj++] = facpi * sqrt(qm_sum / (2*l+1) );
     if (qlcompflag && iqlcomp == il) normfac = 1.0/sqrt(qm_sum);
 
   }
 
   // TODO: 
-  // 1. Need to allocate extra memory in qn[] for this option
-  // 2. Need to add keyword option
+  // 1. [done]Need to allocate extra memory in qnarray[] for this option
+  // 2. [done]Need to add keyword option
   // 3. Need to caclulate Clebsch-Gordan/Wigner 3j coefficients
-  // 4. Compate to bcc values in /Users/athomps/netapp/codes/MatMiner/matminer/matminer/featurizers/boop
+  //     (Can try getting them from boop.py first)
+  // 5. Compate to bcc values in /Users/athomps/netapp/codes/MatMiner/matminer/matminer/featurizers/boop.py
 
-//   // calculate W_l
+  // calculate W_l
 
-//   if (wlflag) {
-//     for (int il = 0; il < nqlist; il++) {
-//       int l = qlist[il];
-//       double wlsum = 0.0;
-//       for(int m1 = 0; m1 < 2*l+1; m1++) {
-//         for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1; m2++)) {
-//           int m = m1 + m2 - l;
-//           qm1qm2_r = qnm_r[il][m1]*qnm_r[il][m2] - qnm_i[il][m1]*qnm_i[il][m2];
-//           qm1qm2_i = qnm_r[il][m1]*qnm_i[il][m2] + qnm_i[il][m1]*qnm_r[il][m2];
-//           wlsum += (qm1qm2_r*qnm_r[il][m] + qm1qm2_i*qnm_i[il][m])*cg; 
-//         }
-//       }
-//       qn[jcount++] = wlsum;
-//     }
-//   }
+  if (wlflag) {
+    int idxcg_count = 0;
+    for (int il = 0; il < nqlist; il++) {
+      int l = qlist[il];
+      double wlsum = 0.0;
+      for(int m1 = 0; m1 < 2*l+1; m1++) {
+        for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {
+          int m = m1 + m2 - l;
+          double qm1qm2_r = qnm_r[il][m1]*qnm_r[il][m2] - qnm_i[il][m1]*qnm_i[il][m2];
+          double qm1qm2_i = qnm_r[il][m1]*qnm_i[il][m2] + qnm_i[il][m1]*qnm_r[il][m2];
+          wlsum += (qm1qm2_r*qnm_r[il][m] + qm1qm2_i*qnm_i[il][m])*cglist[idxcg_count];
+          idxcg_count++;
+        }
+      }
+      // This matches boop.py output, with cg==1
+      printf("W_l = %d %g\n",l,wlsum);
+      qn[jj++] = wlsum/sqrt(2*l+1);
+    }
+  }
 
   // output of the complex vector
 
   if (qlcompflag) {
     for(int m = 0; m < 2*qlcomp+1; m++) {
-      qn[jcount++] = qnm_r[iqlcomp][m] * normfac;
-      qn[jcount++] = qnm_i[iqlcomp][m] * normfac;
+      qn[jj++] = qnm_r[iqlcomp][m] * normfac;
+      qn[jj++] = qnm_i[iqlcomp][m] * normfac;
     }
   }
 
@@ -567,3 +597,272 @@ double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x)
 
   return p;
 }
+
+/* ----------------------------------------------------------------------
+   assign Clebsch-Gordan coefficients with l1=l2=l
+   using the quasi-binomial formula VMK 8.2.1(3)
+------------------------------------------------------------------------- */
+
+void ComputeOrientOrderAtom::init_clebsch_gordan()
+{
+  double sum,dcg,sfaccg;
+  int m, aa2, bb2, cc2, j1, j2, j;
+  int ifac, idxcg_count;
+
+  idxcg_count = 0;
+  for (int il = 0; il < nqlist; il++) {
+    int l = qlist[il];
+    for(int m1 = 0; m1 < 2*l+1; m1++) {
+      for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {
+        idxcg_count++;
+      }
+    }
+  }
+  idxcg_max = idxcg_count;
+  memory->create(cglist, idxcg_max, "computeorientorderatom:cglist");
+
+  idxcg_count = 0;
+  for (int il = 0; il < nqlist; il++) {
+    int l = qlist[il];
+    j1 = j2 = j = l;
+    for(int m1 = 0; m1 < 2*l+1; m1++) {
+        aa2 = m1 - j1;
+        for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {
+          bb2 = m2 - j2;
+          m = (aa2 + bb2 + j) / 2;
+          
+          sum = 0.0;
+          for (int z = MAX(0, MAX(-(j - j2 + aa2)
+                                  , -(j - j1 - bb2) ));
+               z <= MIN((j1 + j2 - j),
+                        MIN((j1 - aa2), (j2 + bb2)));
+               z++) {
+            ifac = z % 2 ? -1 : 1;
+            sum += ifac /
+              (factorial(z) *
+               factorial((j1 + j2 - j) - z) *
+               factorial((j1 - aa2) - z) *
+               factorial((j2 + bb2) - z) *
+               factorial((j - j2 + aa2) + z) *
+               factorial((j - j1 - bb2) + z));
+          }
+          
+          cc2 = m - j;
+          dcg = deltacg(j1, j2, j);
+          sfaccg = sqrt(factorial((j1 + aa2)) *
+                        factorial((j1 - aa2)) *
+                        factorial((j2 + bb2)) *
+                        factorial((j2 - bb2)) *
+                        factorial((j  + cc2)) *
+                        factorial((j  - cc2)) *
+                        (j + 1));
+          
+          cglist[idxcg_count] = sum * dcg * sfaccg;
+          idxcg_count++;
+        }
+      }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   factorial n, wrapper for precomputed table
+------------------------------------------------------------------------- */
+
+double ComputeOrientOrderAtom::factorial(int n)
+{
+  if (n < 0 || n > nmaxfactorial) {
+    char str[128];
+    sprintf(str, "Invalid argument to factorial %d", n);
+    error->all(FLERR, str);
+  }
+
+  return nfac_table[n];
+}
+
+/* ----------------------------------------------------------------------
+   factorial n table, size SNA::nmaxfactorial+1
+------------------------------------------------------------------------- */
+
+const double ComputeOrientOrderAtom::nfac_table[] = {
+  1,
+  1,
+  2,
+  6,
+  24,
+  120,
+  720,
+  5040,
+  40320,
+  362880,
+  3628800,
+  39916800,
+  479001600,
+  6227020800,
+  87178291200,
+  1307674368000,
+  20922789888000,
+  355687428096000,
+  6.402373705728e+15,
+  1.21645100408832e+17,
+  2.43290200817664e+18,
+  5.10909421717094e+19,
+  1.12400072777761e+21,
+  2.5852016738885e+22,
+  6.20448401733239e+23,
+  1.5511210043331e+25,
+  4.03291461126606e+26,
+  1.08888694504184e+28,
+  3.04888344611714e+29,
+  8.8417619937397e+30,
+  2.65252859812191e+32,
+  8.22283865417792e+33,
+  2.63130836933694e+35,
+  8.68331761881189e+36,
+  2.95232799039604e+38,
+  1.03331479663861e+40,
+  3.71993326789901e+41,
+  1.37637530912263e+43,
+  5.23022617466601e+44,
+  2.03978820811974e+46,
+  8.15915283247898e+47,
+  3.34525266131638e+49,
+  1.40500611775288e+51,
+  6.04152630633738e+52,
+  2.65827157478845e+54,
+  1.1962222086548e+56,
+  5.50262215981209e+57,
+  2.58623241511168e+59,
+  1.24139155925361e+61,
+  6.08281864034268e+62,
+  3.04140932017134e+64,
+  1.55111875328738e+66,
+  8.06581751709439e+67,
+  4.27488328406003e+69,
+  2.30843697339241e+71,
+  1.26964033536583e+73,
+  7.10998587804863e+74,
+  4.05269195048772e+76,
+  2.35056133128288e+78,
+  1.3868311854569e+80,
+  8.32098711274139e+81,
+  5.07580213877225e+83,
+  3.14699732603879e+85,
+  1.98260831540444e+87,
+  1.26886932185884e+89,
+  8.24765059208247e+90,
+  5.44344939077443e+92,
+  3.64711109181887e+94,
+  2.48003554243683e+96,
+  1.71122452428141e+98,
+  1.19785716699699e+100,
+  8.50478588567862e+101,
+  6.12344583768861e+103,
+  4.47011546151268e+105,
+  3.30788544151939e+107,
+  2.48091408113954e+109,
+  1.88549470166605e+111,
+  1.45183092028286e+113,
+  1.13242811782063e+115,
+  8.94618213078297e+116,
+  7.15694570462638e+118,
+  5.79712602074737e+120,
+  4.75364333701284e+122,
+  3.94552396972066e+124,
+  3.31424013456535e+126,
+  2.81710411438055e+128,
+  2.42270953836727e+130,
+  2.10775729837953e+132,
+  1.85482642257398e+134,
+  1.65079551609085e+136,
+  1.48571596448176e+138,
+  1.3520015276784e+140,
+  1.24384140546413e+142,
+  1.15677250708164e+144,
+  1.08736615665674e+146,
+  1.03299784882391e+148,
+  9.91677934870949e+149,
+  9.61927596824821e+151,
+  9.42689044888324e+153,
+  9.33262154439441e+155,
+  9.33262154439441e+157,
+  9.42594775983835e+159,
+  9.61446671503512e+161,
+  9.90290071648618e+163,
+  1.02990167451456e+166,
+  1.08139675824029e+168,
+  1.14628056373471e+170,
+  1.22652020319614e+172,
+  1.32464181945183e+174,
+  1.44385958320249e+176,
+  1.58824554152274e+178,
+  1.76295255109024e+180,
+  1.97450685722107e+182,
+  2.23119274865981e+184,
+  2.54355973347219e+186,
+  2.92509369349301e+188,
+  3.3931086844519e+190,
+  3.96993716080872e+192,
+  4.68452584975429e+194,
+  5.5745857612076e+196,
+  6.68950291344912e+198,
+  8.09429852527344e+200,
+  9.8750442008336e+202,
+  1.21463043670253e+205,
+  1.50614174151114e+207,
+  1.88267717688893e+209,
+  2.37217324288005e+211,
+  3.01266001845766e+213,
+  3.8562048236258e+215,
+  4.97450422247729e+217,
+  6.46685548922047e+219,
+  8.47158069087882e+221,
+  1.118248651196e+224,
+  1.48727070609069e+226,
+  1.99294274616152e+228,
+  2.69047270731805e+230,
+  3.65904288195255e+232,
+  5.01288874827499e+234,
+  6.91778647261949e+236,
+  9.61572319694109e+238,
+  1.34620124757175e+241,
+  1.89814375907617e+243,
+  2.69536413788816e+245,
+  3.85437071718007e+247,
+  5.5502938327393e+249,
+  8.04792605747199e+251,
+  1.17499720439091e+254,
+  1.72724589045464e+256,
+  2.55632391787286e+258,
+  3.80892263763057e+260,
+  5.71338395644585e+262,
+  8.62720977423323e+264,
+  1.31133588568345e+267,
+  2.00634390509568e+269,
+  3.08976961384735e+271,
+  4.78914290146339e+273,
+  7.47106292628289e+275,
+  1.17295687942641e+278,
+  1.85327186949373e+280,
+  2.94670227249504e+282,
+  4.71472363599206e+284,
+  7.59070505394721e+286,
+  1.22969421873945e+289,
+  2.0044015765453e+291,
+  3.28721858553429e+293,
+  5.42391066613159e+295,
+  9.00369170577843e+297,
+  1.503616514865e+300, // nmaxfactorial = 167
+};
+
+/* ----------------------------------------------------------------------
+   the function delta given by VMK Eq. 8.2(1)
+------------------------------------------------------------------------- */
+
+double ComputeOrientOrderAtom::deltacg(int j1, int j2, int j)
+{
+  double sfaccg = factorial((j1 + j2 + j) + 1);
+  return sqrt(factorial((j1 + j2 - j)) *
+              factorial((j1 - j2 + j)) *
+              factorial((-j1 + j2 + j)) / sfaccg);
+}
+
diff --git a/src/compute_orientorder_atom.h b/src/compute_orientorder_atom.h
index d5905df63b..b824e736f9 100644
--- a/src/compute_orientorder_atom.h
+++ b/src/compute_orientorder_atom.h
@@ -33,7 +33,7 @@ class ComputeOrientOrderAtom : public Compute {
   void compute_peratom();
   double memory_usage();
   double cutsq;
-  int iqlcomp, qlcomp, qlcompflag;
+  int iqlcomp, qlcomp, qlcompflag, wlflag;
   int *qlist;
   int nqlist;
 
@@ -55,6 +55,14 @@ class ComputeOrientOrderAtom : public Compute {
 
   double polar_prefactor(int, int, double);
   double associated_legendre(int, int, double);
+
+  static const int nmaxfactorial = 167;
+  static const double nfac_table[];
+  double factorial(int);
+  void init_clebsch_gordan();
+  double deltacg(int, int, int);
+  double *cglist;
+  int idxcg_max;
 };
 
 }