diff --git a/doc/src/compute_orientorder_atom.txt b/doc/src/compute_orientorder_atom.txt
index 7327a7b1d3..da14c866bd 100644
--- a/doc/src/compute_orientorder_atom.txt
+++ b/doc/src/compute_orientorder_atom.txt
@@ -19,6 +19,8 @@ keyword = {cutoff} or {nnn} or {degrees} or {components}
   {cutoff} value = distance cutoff
   {nnn} value = number of nearest neighbors
   {degrees} values = nlvalues, l1, l2,...
+  {wl} value = yes or no
+  {wl/hat} value = yes or no
   {components} value = ldegree  :pre
 
 :ule
@@ -27,7 +29,8 @@ keyword = {cutoff} or {nnn} or {degrees} or {components}
 
 compute 1 all orientorder/atom
 compute 1 all orientorder/atom degrees 5 4 6 8 10 12 nnn NULL cutoff 1.5
-compute 1 all orientorder/atom degrees 4 6 components 6 nnn NULL cutoff 3.0 :pre
+compute 1 all orientorder/atom wl/hat yes
+compute 1 all orientorder/atom components 6 :pre
 
 [Description:]
 
@@ -48,7 +51,7 @@ neighbors of the central atom.
 The angles theta and phi are the standard spherical polar angles
 defining the direction of the bond vector {rij}.
 The second equation defines {Ql}, which is a
-rotationally invariant scalar quantity obtained by summing
+rotationally invariant non-negative amplitude obtained by summing
 over all the components of degree {l}.
 
 The optional keyword {cutoff} defines the distance cutoff
@@ -63,7 +66,7 @@ specified distance cutoff are used.
 
 The optional keyword {degrees} defines the list of order parameters to
 be computed.  The first argument {nlvalues} is the number of order
-parameters. This is followed by that number of integers giving the
+parameters. This is followed by that number of non-negative integers giving the
 degree of each order parameter. Because {Q}2 and all odd-degree order
 parameters are zero for atoms in cubic crystals (see
 "Steinhardt"_#Steinhardt), the default order parameters are {Q}4,
@@ -71,7 +74,20 @@ parameters are zero for atoms in cubic crystals (see
 = sqrt(7/3)/8 = 0.19094....  The numerical values of all order
 parameters up to {Q}12 for a range of commonly encountered
 high-symmetry structures are given in Table I of "Mickel et
-al."_#Mickel.
+al."_#Mickel, and these can be reproduced with this compute
+
+The optional keyword {wl} will output the third-order invariants {Wl} 
+(see Eq. 1.4 in "Steinhardt"_#Steinhardt) for the same degrees as
+for the {Ql} parameters. For the FCC crystal with {nnn} =12,
+{W}4 = -sqrt(14/143).(49/4096)/Pi^1.5 = -0.0006722136...
+
+The optional keyword {wl/hat} will output the normalized third-order 
+invariants {Wlhat} (see Eq. 2.2 in "Steinhardt"_#Steinhardt) 
+for the same degrees as for the {Ql} parameters. For the FCC crystal 
+with {nnn} =12, {W}4hat = -7/3*sqrt(2/429) = -0.159317...The numerical
+values of {Wlhat} for a range of commonly encountered high-symmetry
+structures are given in Table I of "Steinhardt"_#Steinhardt, and these
+can be reproduced with this keyword.
 
 The optional keyword {components} will output the components of the
 normalized complex vector {Ybar_lm} of degree {ldegree}, which must be
@@ -82,7 +98,7 @@ particles, as discussed in "ten Wolde"_#tenWolde2.
 
 The value of {Ql} is set to zero for atoms not in the
 specified compute group, as well as for atoms that have less than
-{nnn} neighbors within the distance cutoff.
+{nnn} neighbors within the distance cutoff, unless {nnn} is NULL.
 
 The neighbor list needed to compute this quantity is constructed each
 time the calculation is performed (i.e. each time a snapshot of atoms
@@ -108,6 +124,12 @@ This compute calculates a per-atom array with {nlvalues} columns,
 giving the {Ql} values for each atom, which are real numbers on the
 range 0 <= {Ql} <= 1.
 
+If the keyword {wl} is set to yes, then the {Wl} values for each
+atom will be added to the output array, which are real numbers.
+
+If the keyword {wl/hat} is set to yes, then the {Wl_hat} 
+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) {Ybar_lm} will be added to the
 output array in the following order: Re({Ybar_-m}) Im({Ybar_-m})
@@ -130,7 +152,8 @@ hexorder/atom"_compute_hexorder_atom.html
 [Default:]
 
 The option defaults are {cutoff} = pair style cutoff, {nnn} = 12,
-{degrees} = 5 4 6 8 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12.
+{degrees} = 5 4 6 8 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12,
+{wl} = no, {wl/hat} = no, and {components} off
 
 :line
 
diff --git a/examples/steinhardt/in.bcc b/examples/steinhardt/in.bcc
new file mode 100644
index 0000000000..a0348fdcd9
--- /dev/null
+++ b/examples/steinhardt/in.bcc
@@ -0,0 +1,91 @@
+# Steinhardt-Nelson bond orientational order parameters for BCC
+
+variable  	rcut equal 3.0
+
+boundary	p p p
+
+atom_style	atomic
+neighbor	0.3 bin
+neigh_modify	delay 5
+
+# create geometry
+
+lattice		bcc 1.0
+region		box block 0 3 0 3 0 3
+create_box	1 box
+create_atoms	1 box
+
+mass		1 1.0
+
+# LJ potentials
+
+pair_style	lj/cut ${rcut}
+pair_coeff	* * 1.0 1.0 ${rcut}
+
+# 14 neighbors, perfect crystal 
+
+compute 	qlwlhat all orientorder/atom degrees 6 2 4 6 8 10 12 nnn 14 wl/hat yes
+compute 	avql all reduce ave c_qlwlhat[1] c_qlwlhat[2] c_qlwlhat[3] c_qlwlhat[4] c_qlwlhat[5] c_qlwlhat[6] 
+compute 	avwlhat all reduce ave c_qlwlhat[7] c_qlwlhat[8] c_qlwlhat[9] c_qlwlhat[10] c_qlwlhat[11] c_qlwlhat[12] 
+
+thermo_style    custom step temp epair etotal c_avql[*] c_avwlhat[*]
+
+run		0
+
+# check Q_l values
+
+print " "
+print "*******************************************************************"
+print " "
+print "Comparison with reference values of Q_l "
+print "   [Table I in W. Mickel, S. C. Kapfer," 
+print "   G. E. Schroeder-Turkand, K. Mecke, "
+print "   J. Chem. Phys. 138, 044501 (2013).]"
+print " "
+
+variable 	q2ref equal 0.0
+variable 	q4ref equal 0.036
+variable 	q6ref equal 0.511
+variable 	q8ref equal 0.429
+variable 	q10ref equal 0.195
+variable 	q12ref equal 0.405
+
+variable 	q2 equal c_avql[1]
+variable 	q4 equal c_avql[2]
+variable 	q6 equal c_avql[3]
+variable 	q8 equal c_avql[4]
+variable 	q10 equal c_avql[5]
+variable 	q12 equal c_avql[6]
+
+print "q2     = $(v_q2:%10.6f) delta = $(v_q2-v_q2ref:%10.4f)"
+print "q4     = $(v_q4:%10.6f) delta = $(v_q4-v_q4ref:%10.4f)"
+print "q6     = $(v_q6:%10.6f) delta = $(v_q6-v_q6ref:%10.4f)"
+print "q8     = $(v_q8:%10.6f) delta = $(v_q8-v_q8ref:%10.4f)"
+print "q10    = $(v_q10:%10.6f) delta = $(v_q10-v_q10ref:%10.4f)"
+print "q12    = $(v_q12:%10.6f) delta = $(v_q12-v_q12ref:%10.4f)"
+
+# check W_l_hat values
+
+print " "
+print "Comparison with reference values of W_l_hat"
+print "   [Table I in P. Steinhardt, D. Nelson, and M. Ronchetti, "
+print "   Phys. Rev. B 28, 784 (1983).]"
+print " "
+
+variable 	w4hatref equal 0.159317
+variable 	w6hatref equal 0.013161
+variable 	w8hatref equal -0.058455
+variable 	w10hatref equal -0.090130
+
+variable 	w4hat equal c_avwlhat[2]
+variable 	w6hat equal c_avwlhat[3]
+variable 	w8hat equal c_avwlhat[4]
+variable 	w10hat equal c_avwlhat[5]
+
+print "w4hat  = $(v_w4hat:%10.6f) delta = $(v_w4hat-v_w4hatref:%10.6f)"
+print "w6hat  = $(v_w6hat:%10.6f) delta = $(v_w6hat-v_w6hatref:%10.6f)"
+print "w8hat  = $(v_w8hat:%10.6f) delta = $(v_w8hat-v_w8hatref:%10.6f)"
+print "w10hat = $(v_w10hat:%10.6f) delta = $(v_w10hat-v_w10hatref:%10.6f)"
+print " "
+print "*******************************************************************"
+print " "
diff --git a/examples/steinhardt/in.fcc b/examples/steinhardt/in.fcc
new file mode 100644
index 0000000000..6d2775d0bb
--- /dev/null
+++ b/examples/steinhardt/in.fcc
@@ -0,0 +1,88 @@
+# Steinhardt-Nelson bond orientational order parameters for FCC
+
+variable  	rcut equal 3.0
+
+boundary	p p p
+
+atom_style	atomic
+neighbor	0.3 bin
+neigh_modify	delay 5
+
+# create geometry
+
+lattice		fcc 1.0
+region		box block 0 3 0 3 0 3
+create_box	1 box
+create_atoms	1 box
+
+mass		1 1.0
+
+# LJ potentials
+
+pair_style	lj/cut ${rcut}
+pair_coeff	* * 1.0 1.0 ${rcut}
+
+# 12 neighbors, perfect crystal
+
+compute 	qlwlhat all orientorder/atom wl/hat yes
+compute 	avql all reduce ave c_qlwlhat[1] c_qlwlhat[2] c_qlwlhat[3] c_qlwlhat[4] c_qlwlhat[5] 
+compute 	avwlhat all reduce ave c_qlwlhat[6] c_qlwlhat[7] c_qlwlhat[8] c_qlwlhat[9] c_qlwlhat[10]
+
+thermo_style    custom step temp epair etotal c_avql[*] c_avwlhat[*]
+
+run		0
+
+# check Q_l values
+
+print " "
+print "*******************************************************************"
+print " "
+print "Comparison with reference values of Q_l "
+print "   [Table I in W. Mickel, S. C. Kapfer," 
+print "   G. E. Schroeder-Turkand, K. Mecke, "
+print "   J. Chem. Phys. 138, 044501 (2013).]"
+print " "
+
+variable 	q4ref equal 0.190
+variable 	q6ref equal 0.575
+variable 	q8ref equal 0.404
+variable 	q10ref equal 0.013
+variable 	q12ref equal 0.600
+
+variable 	q4 equal c_avql[1]
+variable 	q6 equal c_avql[2]
+variable 	q8 equal c_avql[3]
+variable 	q10 equal c_avql[4]
+variable 	q12 equal c_avql[5]
+
+print "q4     = $(v_q4:%10.6f) delta = $(v_q4-v_q4ref:%10.4f)"
+print "q6     = $(v_q6:%10.6f) delta = $(v_q6-v_q6ref:%10.4f)"
+print "q8     = $(v_q8:%10.6f) delta = $(v_q8-v_q8ref:%10.4f)"
+print "q10    = $(v_q10:%10.6f) delta = $(v_q10-v_q10ref:%10.4f)"
+print "q12    = $(v_q12:%10.6f) delta = $(v_q12-v_q12ref:%10.4f)"
+
+# check W_l_hat values
+
+print " "
+print "Comparison with reference values of W_l_hat"
+print "   [Table I in P. Steinhardt, D. Nelson, and M. Ronchetti, "
+print "   Phys. Rev. B 28, 784 (1983).]"
+print " "
+
+variable 	w4hatref equal -0.159316
+variable 	w6hatref equal -0.013161
+variable 	w8hatref equal 0.058454
+variable 	w10hatref equal -0.090128
+
+variable 	w4hat equal c_avwlhat[1]
+variable 	w6hat equal c_avwlhat[2]
+variable 	w8hat equal c_avwlhat[3]
+variable 	w10hat equal c_avwlhat[4]
+
+print "w4hat  = $(v_w4hat:%10.6f) delta = $(v_w4hat-v_w4hatref:%10.6f)"
+print "w6hat  = $(v_w6hat:%10.6f) delta = $(v_w6hat-v_w6hatref:%10.6f)"
+print "w8hat  = $(v_w8hat:%10.6f) delta = $(v_w8hat-v_w8hatref:%10.6f)"
+print "w10hat = $(v_w10hat:%10.6f) delta = $(v_w10hat-v_w10hatref:%10.6f)"
+print " "
+print "*******************************************************************"
+print " "
diff --git a/examples/steinhardt/in.icos b/examples/steinhardt/in.icos
new file mode 100644
index 0000000000..d0d61a902d
--- /dev/null
+++ b/examples/steinhardt/in.icos
@@ -0,0 +1,106 @@
+# Steinhardt-Nelson bond orientational order parameters for icosahedral cluster
+# W_6_hat is sensitive to icosohedral order
+
+variable        rcut equal 1.2 # a bit bigger than LJ Rmin 
+variable        rcutred equal 0.75 # a bit bigger than 1/sqrt(2)
+
+# create a perfect fcc crystallite
+
+atom_style	atomic
+boundary	s s s
+lattice		fcc 1.0 # neighbors at LJ Rmin
+region		box block 0 2 0 2 0 2
+create_box	1 box
+create_atoms	1 box
+mass		1 1.0
+
+region 		centralatom sphere 1 1 1 0.0 side in
+group 		centralatom region centralatom
+
+region 		mysphere sphere 1 1 1 ${rcutred} side out
+delete_atoms 	region mysphere
+
+# LJ potential
+
+pair_style	lj/cut 100.0
+pair_coeff	* * 1.0 1.0 100.0
+
+# define output for central atom
+
+compute 	qlwlhat all orientorder/atom wl/hat yes cutoff ${rcut} nnn NULL
+compute 	avql centralatom reduce ave c_qlwlhat[1] c_qlwlhat[2] c_qlwlhat[3] c_qlwlhat[4] c_qlwlhat[5] 
+compute 	avwlhat centralatom reduce ave c_qlwlhat[6] c_qlwlhat[7] c_qlwlhat[8] c_qlwlhat[9] c_qlwlhat[10] 
+variable	q6 equal c_avql[2]
+variable	w6hat equal c_avwlhat[2]
+
+compute 	mype all pe/atom
+compute 	centralatompe centralatom reduce ave c_mype
+
+# gently equilibrate the crystallite
+
+velocity     	all create 0.001 482748
+fix	     	1 all nve
+neighbor	0.3 bin
+neigh_modify	every 1 check no delay 0
+timestep	0.003
+thermo_style    custom step temp epair etotal c_centralatompe v_q6 v_w6hat
+thermo		10
+
+run		10
+
+# quench to icosehedral cluster
+
+minimize      	1.0e-10 1.0e-6 100 1000
+
+# check Q_l values
+
+print " "
+print "*******************************************************************"
+print " "
+print "Comparison with reference values of Q_l "
+print "   [Table I in W. Mickel, S. C. Kapfer," 
+print "   G. E. Schroeder-Turkand, K. Mecke, "
+print "   J. Chem. Phys. 138, 044501 (2013).]"
+print " "
+
+variable 	q4ref equal 0.0
+variable 	q6ref equal 0.663
+variable 	q8ref equal 0.0
+variable 	q10ref equal 0.363
+variable 	q12ref equal 0.585
+
+variable 	q4 equal c_avql[1]
+variable 	q6 equal c_avql[2]
+variable 	q8 equal c_avql[3]
+variable 	q10 equal c_avql[4]
+variable 	q12 equal c_avql[5]
+
+print "q4     = $(v_q4:%10.6f) delta = $(v_q4-v_q4ref:%10.4f)"
+print "q6     = $(v_q6:%10.6f) delta = $(v_q6-v_q6ref:%10.4f)"
+print "q8     = $(v_q8:%10.6f) delta = $(v_q8-v_q8ref:%10.4f)"
+print "q10    = $(v_q10:%10.6f) delta = $(v_q10-v_q10ref:%10.4f)"
+print "q12    = $(v_q12:%10.6f) delta = $(v_q12-v_q12ref:%10.4f)"
+
+# check W_l_hat values
+
+print " "
+print "Comparison with reference values of W_l_hat"
+print "   [Table I in P. Steinhardt, D. Nelson, and M. Ronchetti, "
+print "   Phys. Rev. B 28, 784 (1983).]"
+print " "
+
+variable 	w6hatref equal -0.169754
+variable 	w10hatref equal -0.093967
+
+variable 	w4hat equal c_avwlhat[1]
+variable 	w6hat equal c_avwlhat[2]
+variable 	w8hat equal c_avwlhat[3]
+variable 	w10hat equal c_avwlhat[4]
+variable 	w12hat equal c_avwlhat[5]
+
+print "w6hat  = $(v_w6hat:%10.6f) delta = $(v_w6hat-v_w6hatref:%10.6f)"
+print "w10hat = $(v_w10hat:%10.6f) delta = $(v_w10hat-v_w10hatref:%10.6f)"
+print " "
+print "*******************************************************************"
+print " "
+
diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp
index 266df575f9..0a78356127 100644
--- a/src/compute_orientorder_atom.cpp
+++ b/src/compute_orientorder_atom.cpp
@@ -43,12 +43,14 @@ using namespace std;
   #define MY_EPSILON (10.0*2.220446049250313e-16)
 #endif
 
+#define QEPSILON 1.0e-6
+
 /* ---------------------------------------------------------------------- */
 
 ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg),
   qlist(NULL), distsq(NULL), nearest(NULL), rlist(NULL),
-  qnarray(NULL), qnm_r(NULL), qnm_i(NULL)
+  qnarray(NULL), qnm_r(NULL), qnm_i(NULL), cglist(NULL)
 {
   if (narg < 3 ) error->all(FLERR,"Illegal compute orientorder/atom command");
 
@@ -56,6 +58,8 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
 
   nnn = 12;
   cutsq = 0.0;
+  wlflag = 0;
+  wlhatflag = 0;
   qlcompflag = 0;
 
   // specify which orders to request
@@ -96,27 +100,39 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
       if (iarg+nqlist > narg)
         error->all(FLERR,"Illegal compute orientorder/atom command");
       qmax = 0;
-      for (int iw = 0; iw < nqlist; iw++) {
-        qlist[iw] = force->numeric(FLERR,arg[iarg+iw]);
-        if (qlist[iw] < 0)
+      for (int il = 0; il < nqlist; il++) {
+        qlist[il] = force->numeric(FLERR,arg[iarg+il]);
+        if (qlist[il] < 0)
           error->all(FLERR,"Illegal compute orientorder/atom command");
-        if (qlist[iw] > qmax) qmax = qlist[iw];
+        if (qlist[il] > qmax) qmax = qlist[il];
       }
       iarg += nqlist;
+    } else if (strcmp(arg[iarg],"wl") == 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],"wl/hat") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute orientorder/atom command");
+      if (strcmp(arg[iarg+1],"yes") == 0) wlhatflag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) wlhatflag = 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)
         error->all(FLERR,"Illegal compute orientorder/atom command");
       qlcomp = force->numeric(FLERR,arg[iarg+1]);
-      if (qlcomp <= 0)
-        error->all(FLERR,"Illegal compute orientorder/atom command");
       iqlcomp = -1;
-      for (int iw = 0; iw < nqlist; iw++)
-        if (qlcomp == qlist[iw]) {
-          iqlcomp = iw;
+      for (int il = 0; il < nqlist; il++)
+        if (qlcomp == qlist[il]) {
+          iqlcomp = il;
           break;
         }
-      if (iqlcomp < 0)
+      if (iqlcomp == -1)
         error->all(FLERR,"Illegal compute orientorder/atom command");
       iarg += 2;
     } else if (strcmp(arg[iarg],"cutoff") == 0) {
@@ -130,8 +146,10 @@ 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 (wlhatflag) ncol += nqlist;
+  if (qlcompflag) ncol += 2*(2*qlcomp+1);
 
   peratom_flag = 1;
   size_peratom_cols = ncol;
@@ -151,7 +169,7 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
   memory->destroy(qlist);
   memory->destroy(qnm_r);
   memory->destroy(qnm_i);
-
+  memory->destroy(cglist);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -166,8 +184,8 @@ void ComputeOrientOrderAtom::init()
     error->all(FLERR,"Compute orientorder/atom cutoff is "
                "longer than pairwise cutoff");
 
-  memory->create(qnm_r,qmax,2*qmax+1,"orientorder/atom:qnm_r");
-  memory->create(qnm_i,qmax,2*qmax+1,"orientorder/atom:qnm_i");
+  memory->create(qnm_r,nqlist,2*qmax+1,"orientorder/atom:qnm_r");
+  memory->create(qnm_i,nqlist,2*qmax+1,"orientorder/atom:qnm_i");
 
   // need an occasional full neighbor list
 
@@ -183,6 +201,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 || wlhatflag) init_clebsch_gordan();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -274,8 +294,8 @@ void ComputeOrientOrderAtom::compute_peratom()
       // if not nnn neighbors, order parameter = 0;
 
       if ((ncount == 0) || (ncount < nnn)) {
-        for (int iw = 0; iw < nqlist; iw++)
-          qn[iw] = 0.0;
+        for (int jj = 0; jj < ncol; jj++)
+          qn[jj] = 0.0;
         continue;
       }
 
@@ -287,6 +307,7 @@ void ComputeOrientOrderAtom::compute_peratom()
       }
 
       calc_boop(rlist, ncount, qn, qlist, nqlist);
+
     }
   }
 }
@@ -403,13 +424,12 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
 void ComputeOrientOrderAtom::calc_boop(double **rlist,
                                        int ncount, double qn[],
                                        int qlist[], int nqlist) {
-  for (int iw = 0; iw < nqlist; iw++) {
-    int n = qlist[iw];
 
-    qn[iw] = 0.0;
-    for(int m = 0; m < 2*n+1; m++) {
-      qnm_r[iw][m] = 0.0;
-      qnm_i[iw][m] = 0.0;
+  for (int il = 0; il < nqlist; il++) {
+    int l = qlist[il];
+    for(int m = 0; m < 2*l+1; m++) {
+      qnm_r[il][m] = 0.0;
+      qnm_i[il][m] = 0.0;
     }
   }
 
@@ -433,24 +453,24 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
       expphi_i *= rxymaginv;
     }
 
-    for (int iw = 0; iw < nqlist; iw++) {
-      int n = qlist[iw];
+    for (int il = 0; il < nqlist; il++) {
+      int l = qlist[il];
 
-      qnm_r[iw][n] += polar_prefactor(n, 0, costheta);
+      qnm_r[il][l] += polar_prefactor(l, 0, costheta);
       double expphim_r = expphi_r;
       double expphim_i = expphi_i;
-      for(int m = 1; m <= +n; m++) {
-        double prefactor = polar_prefactor(n, m, costheta);
+      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[iw][m+n] += c_r;
-        qnm_i[iw][m+n] += c_i;
+        qnm_r[il][m+l] += c_r;
+        qnm_i[il][m+l] += c_i;
         if(m & 1) {
-          qnm_r[iw][-m+n] -= c_r;
-          qnm_i[iw][-m+n] += c_i;
+          qnm_r[il][-m+l] -= c_r;
+          qnm_i[il][-m+l] += c_i;
         } else {
-          qnm_r[iw][-m+n] += c_r;
-          qnm_i[iw][-m+n] -= c_i;
+          qnm_r[il][-m+l] += c_r;
+          qnm_i[il][-m+l] -= c_i;
         }
         double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i;
         double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r;
@@ -461,30 +481,110 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
     }
   }
 
-  double fac = sqrt(MY_4PI) / ncount;
-  double normfac = 0.0;
-  for (int iw = 0; iw < nqlist; iw++) {
-    int n = qlist[iw];
-    double qm_sum = 0.0;
-    for(int m = 0; m < 2*n+1; m++) {
-      qm_sum += qnm_r[iw][m]*qnm_r[iw][m] + qnm_i[iw][m]*qnm_i[iw][m];
-      //      printf("Ylm^2 = %d %d %g\n",n,m,
-      //     qnm_r[iw][m]*qnm_r[iw][m] + qnm_i[iw][m]*qnm_i[iw][m]);
-    }
-    qn[iw] = fac * sqrt(qm_sum / (2*n+1));
-    if (qlcompflag && iqlcomp == iw) normfac = 1.0/sqrt(qm_sum);
+  // 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;
+    }
   }
 
-  // output of the complex vector
+  // calculate Q_l
+  // NOTE: optional W_l_hat and components of Q_qlcomp use these stored Q_l values 
+
+  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 = 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[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 caclulate 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) {
+    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++;
+        }
+      }
+      qn[jj++] = wlsum/sqrt(2*l+1);
+    }
+  }
+
+  // calculate W_l_hat
+
+  if (wlhatflag) {
+    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++;
+        }
+      }
+  //      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 {
+        double qnormfac = sqrt(MY_4PI/(2*l+1));
+        double qnfac = qnormfac/qn[il];
+        qn[jj++] = wlsum/sqrt(2*l+1)*(qnfac*qnfac*qnfac);
+      }
+    }
+  }
+
+  // Calculate components of Q_l, for l=qlcomp
 
   if (qlcompflag) {
-    int j = nqlist;
-    for(int m = 0; m < 2*qlcomp+1; m++) {
-      qn[j++] = qnm_r[iqlcomp][m] * normfac;
-      qn[j++] = qnm_i[iqlcomp][m] * normfac;
+    int il = iqlcomp;
+    int l = qlcomp;
+    if (qn[il] < QEPSILON)
+      for(int m = 0; m < 2*l+1; m++) {
+        qn[jj++] = 0.0;
+        qn[jj++] = 0.0;
+      }
+    else {
+      double qnormfac = sqrt(MY_4PI/(2*l+1));
+      double qnfac = qnormfac/qn[il];
+      for(int m = 0; m < 2*l+1; m++) {
+        qn[jj++] = qnm_r[il][m] * qnfac;
+        qn[jj++] = qnm_i[il][m] * qnfac;
+      }
     }
   }
+
 }
 
 /* ----------------------------------------------------------------------
@@ -542,3 +642,258 @@ double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x)
 
   return p;
 }
+
+/* ----------------------------------------------------------------------
+   assign Clebsch-Gordan coefficients
+   using the quasi-binomial formula VMK 8.2.1(3)
+   specialized for case j1=j2=j=l
+------------------------------------------------------------------------- */
+
+void ComputeOrientOrderAtom::init_clebsch_gordan()
+{
+  double sum,dcg,sfaccg, sfac1, sfac2;
+  int m, aa2, bb2, cc2;
+  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];
+    for(int m1 = 0; m1 < 2*l+1; m1++) {
+        aa2 = m1 - l;
+        for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {
+          bb2 = m2 - l;
+          m = aa2 + bb2 + l;
+          
+          sum = 0.0;
+          for (int z = MAX(0, MAX(-aa2, bb2));
+               z <= MIN(l, MIN(l - aa2, l + bb2)); z++) {
+            ifac = z % 2 ? -1 : 1;
+            sum += ifac /
+              (factorial(z) *
+               factorial(l - z) *
+               factorial(l - aa2 - z) *
+               factorial(l + bb2 - z) *
+               factorial(aa2 + z) *
+               factorial(-bb2 + z));
+          }
+          
+          cc2 = m - l;
+          sfaccg = sqrt(factorial(l + aa2) *
+                        factorial(l - aa2) *
+                        factorial(l + bb2) *
+                        factorial(l - bb2) *
+                        factorial(l + cc2) *
+                        factorial(l - cc2) *
+                        (2*l + 1));
+          
+          sfac1 = factorial(3*l + 1);
+          sfac2 = factorial(l);
+          dcg = sqrt(sfac2*sfac2*sfac2 / sfac1);
+
+          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
+};
+
diff --git a/src/compute_orientorder_atom.h b/src/compute_orientorder_atom.h
index d5905df63b..643875ccd0 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, wlhatflag;
   int *qlist;
   int nqlist;
 
@@ -55,6 +55,13 @@ 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 *cglist;                      // Clebsch-Gordan coeffs
+  int idxcg_max;
 };
 
 }