diff --git a/doc/src/compute_orientorder_atom.txt b/doc/src/compute_orientorder_atom.txt
index 2e83f3ea07..f4e1cbb924 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
@@ -63,7 +65,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,
diff --git a/examples/steinhardt/in.bcc b/examples/steinhardt/in.bcc
index acdfca3647..a0348fdcd9 100644
--- a/examples/steinhardt/in.bcc
+++ b/examples/steinhardt/in.bcc
@@ -22,26 +22,70 @@ mass		1 1.0
 pair_style	lj/cut ${rcut}
 pair_coeff	* * 1.0 1.0 ${rcut}
 
-# initial velocities
+# 14 neighbors, perfect crystal 
 
-velocity     all create 5.0 482748
-fix	     1 all nve
+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] 
 
-# 8 neighbors, perfect crystal
-
-compute 	qn all orientorder/atom degrees 5 4 6 8 10 12 nnn 8
-compute 	avqn all reduce ave c_qn[*] 
-
-thermo_style    custom step temp epair etotal c_avqn[*]
+thermo_style    custom step temp epair etotal c_avql[*] c_avwlhat[*]
 
 run		0
 
-# 14 neighbors, perfect crystal dynamically melting
+# check Q_l values
 
-uncompute       qn
-compute 	qn all orientorder/atom degrees 5 4 6 8 10 12 nnn 14
+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 " "
 
-timestep	0.003
-thermo		1
+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
 
-run		20
+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
index 6c6a919d7b..6d2775d0bb 100644
--- a/examples/steinhardt/in.fcc
+++ b/examples/steinhardt/in.fcc
@@ -22,18 +22,67 @@ mass		1 1.0
 pair_style	lj/cut ${rcut}
 pair_coeff	* * 1.0 1.0 ${rcut}
 
-# initial velocities
+# 12 neighbors, perfect crystal
 
-velocity     all create 5.0 482748
-fix	     1 all nve
+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]
 
-# 12 neighbors, perfect crystal dynamically melting
+thermo_style    custom step temp epair etotal c_avql[*] c_avwlhat[*]
 
-compute 	qn all orientorder/atom # degrees 5 4 6 8 10 12 nnn 12
-compute 	avqn all reduce ave c_qn[*] 
+run		0
 
-timestep	0.003
-thermo_style    custom step temp epair etotal c_avqn[*]
-thermo		1
+# check Q_l values
 
-run		20
+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 " "
+