diff --git a/doc/Eqs/hexorder.jpg b/doc/Eqs/hexorder.jpg
new file mode 100644
index 0000000000..cebb8d9d9e
Binary files /dev/null and b/doc/Eqs/hexorder.jpg differ
diff --git a/doc/Eqs/hexorder.tex b/doc/Eqs/hexorder.tex
new file mode 100644
index 0000000000..09fd4e2706
--- /dev/null
+++ b/doc/Eqs/hexorder.tex
@@ -0,0 +1,8 @@
+\documentclass[12pt]{article}
+\begin{document}
+
+$$
+   q_6 = \frac{1}{N_{neigh}}\sum_{j = 1}^{N_{neigh}} e^{6 i \theta({\bf r}_{ij})}
+$$
+
+\end{document}
diff --git a/doc/compute_hexorder_atom.txt b/doc/compute_hexorder_atom.txt
new file mode 100644
index 0000000000..5a5171dc50
--- /dev/null
+++ b/doc/compute_hexorder_atom.txt
@@ -0,0 +1,117 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+compute hexorder/atom command :h3
+
+[Syntax:]
+
+compute ID group-ID hexorder/atom cutoff type1 type2 ... :pre
+
+ID, group-ID are documented in "compute"_compute.html command
+hexorder/atom = style name of this compute command
+cutoff = distance within which to count neighbors (distance units)
+typeN = atom type for Nth order parameter (see asterisk form below) :ul
+
+[Examples:]
+
+compute 1 all hexorder/atom 2.0
+compute 1 all hexorder/atom 6.0 1 2
+compute 1 all hexorder/atom 6.0 2*4 5*8 * :pre
+
+[Description:]
+
+Define a computation that calculates one or more hexatic bond orientational 
+order parameters for each atom in a group. The hexatic bond orientational order 
+parameter {q}6 "(Nelson)"_#Nelson for an atom is a complex number (stored as two real numbers).
+It is defined as follows:
+
+:c,image(Eqs/hexorder.jpg)
+
+where the sum is over atoms of the specified atom type(s) that are within 
+the specified cutoff distance from the central atom. The angle theta
+is formed by the bond vector rij and the {x} axis. theta is calculated
+only using the {x} and {y} components, whereas the distance from the
+central atom is calculated using all three  
+{x}, {y}, and {z} components of the bond vector.
+Atoms not in the group 
+are included in the order parameter of atoms in the group. 
+
+If the neighbors of the central atom lie on a hexagonal lattice, 
+then |{q}6| = 1.  
+The complex phase of {q}6 depends on the orientation of the 
+lattice relative to the {x} axis. For a liquid in which the 
+atomic neighborhood lacks orientational symmettry, |{q}6| << 1.
+
+The {typeN} keywords allow you to specify which atom types contribute
+to each order parameter.  One order parameter is computed for
+each of the {typeN} keywords listed.  If no {typeN} keywords are
+listed, a single order parameter is calculated, which includes
+atoms of all types (same as the "*" format, see below).
+
+The {typeN} keywords can be specified in one of two ways.  An explicit
+numeric value can be used, as in the 2nd example above.  Or a
+wild-card asterisk can be used to specify a range of atom types.  This
+takes the form "*" or "*n" or "n*" or "m*n".  If N = the number of
+atom types, then an asterisk with no numeric values means all types
+from 1 to N.  A leading asterisk means all types from 1 to n
+(inclusive).  A trailing asterisk means all types from n to N
+(inclusive).  A middle asterisk means all types from m to n
+(inclusive).
+
+The value of all order parameters will be 0.0+0.0i for atoms not in the
+specified compute group. An order parameter for atoms that have no 
+neighbors of the specified atom type within the cutoff distance will
+be 0.0+0.0i.
+
+The neighbor list needed to compute this quantity is constructed each
+time the calculation is performed (i.e. each time a snapshot of atoms
+is dumped).  Thus it can be inefficient to compute/dump this quantity
+too frequently.
+
+IMPORTANT NOTE: If you have a bonded system, then the settings of
+"special_bonds"_special_bonds.html command can remove pairwise
+interactions between atoms in the same bond, angle, or dihedral.  This
+is the default setting for the "special_bonds"_special_bonds.html
+command, and means those pairwise interactions do not appear in the
+neighbor list.  Because this fix uses the neighbor list, it also means
+those pairs will not be included in the order parameter.  One way
+to get around this, is to write a dump file, and use the
+"rerun"_rerun.html command to compute the order parameter for snapshots
+in the dump file.  The rerun script can use a
+"special_bonds"_special_bonds.html command that includes all pairs in
+the neighbor list.
+
+[Output info:]
+
+If single {type1} keyword is specified (or if none are specified),
+this compute calculates a per-atom array with 2 columns, giving the
+real and imaginary parts of the order parameter, respectively.  
+If multiple {typeN} keywords are specified, this compute calculates 
+a per-atom array with 2*N columns, with each consecutive pair of 
+columns giving the real and imaginary parts of one order parameter.  
+
+These values can be accessed by any command that uses
+per-atom values from a compute as input.  See "Section_howto
+15"_Section_howto.html#howto_15 for an overview of LAMMPS output
+options.
+
+The per-atom array values will be floating point numbers, as
+explained above.
+
+[Restrictions:] none
+
+[Related commands:]
+
+"compute coord/atom"_compute_coord_atom.html
+
+[Default:] none
+
+:line
+
+:link(Nelson)
+[(Nelson)] Nelson, Halperin, Phys Rev B, 19, 2457 (1979).
diff --git a/src/compute_hexorder_atom.cpp b/src/compute_hexorder_atom.cpp
new file mode 100644
index 0000000000..1fedc33c23
--- /dev/null
+++ b/src/compute_hexorder_atom.cpp
@@ -0,0 +1,263 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author:  Aidan Thompson (SNL)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+#include "compute_hexorder_atom.h"
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "force.h"
+#include "pair.h"
+#include "comm.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg < 4) error->all(FLERR,"Illegal compute hexorder/atom command");
+
+  double cutoff = force->numeric(FLERR,arg[3]);
+  cutsq = cutoff*cutoff;
+
+  ncol = narg-4 + 1;
+  int ntypes = atom->ntypes;
+  typelo = new int[ncol];
+  typehi = new int[ncol];
+
+  if (narg == 4) {
+    ncol = 2;
+    typelo[0] = 1;
+    typehi[0] = ntypes;
+  } else {
+    ncol = 0;
+    int iarg = 4;
+    while (iarg < narg) {
+      force->bounds(arg[iarg],ntypes,typelo[ncol],typehi[ncol]);
+      if (typelo[ncol] > typehi[ncol])
+        error->all(FLERR,"Illegal compute hexorder/atom command");
+      typelo[ncol+1] = typelo[ncol];
+      typehi[ncol+1] = typehi[ncol];
+      ncol+=2;
+      iarg++;
+    }
+  }
+
+  peratom_flag = 1;
+  size_peratom_cols = ncol;
+
+  nmax = 0;
+  q6array = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeHexOrderAtom::~ComputeHexOrderAtom()
+{
+  delete [] typelo;
+  delete [] typehi;
+  memory->destroy(q6array);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHexOrderAtom::init()
+{
+  if (force->pair == NULL)
+    error->all(FLERR,"Compute hexorder/atom requires a pair style be defined");
+  if (sqrt(cutsq) > force->pair->cutforce)
+    error->all(FLERR,
+               "Compute hexorder/atom cutoff is longer than pairwise cutoff");
+
+  // need an occasional full neighbor list
+
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->compute = 1;
+  neighbor->requests[irequest]->half = 0;
+  neighbor->requests[irequest]->full = 1;
+  neighbor->requests[irequest]->occasional = 1;
+
+  int count = 0;
+  for (int i = 0; i < modify->ncompute; i++)
+    if (strcmp(modify->compute[i]->style,"hexorder/atom") == 0) count++;
+  if (count > 1 && comm->me == 0)
+    error->warning(FLERR,"More than one compute hexorder/atom");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHexOrderAtom::init_list(int id, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHexOrderAtom::compute_peratom()
+{
+  int i,j,m,ii,jj,inum,jnum,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+  double *count;
+
+  invoked_peratom = update->ntimestep;
+
+  // grow coordination array if necessary
+
+  if (atom->nlocal > nmax) {
+    memory->destroy(q6array);
+    nmax = atom->nmax;
+    memory->create(q6array,nmax,ncol,"hexorder/atom:q6array");
+    array_atom = q6array;
+  }
+
+  // invoke full neighbor list (will copy or build if necessary)
+
+  neighbor->build_one(list);
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // compute order parameter(s) for each atom in group
+  // use full neighbor list to count atoms less than cutoff
+
+  double **x = atom->x;
+  int *type = atom->type;
+  int *mask = atom->mask;
+
+  if (ncol == 2) {
+    for (ii = 0; ii < inum; ii++) {
+      i = ilist[ii];
+      double* q6 = q6array[i];
+      q6[0] = q6[1] = 0.0;
+      if (mask[i] & groupbit) {
+        xtmp = x[i][0];
+        ytmp = x[i][1];
+        ztmp = x[i][2];
+        jlist = firstneigh[i];
+        jnum = numneigh[i];
+
+	double usum = 0.0;
+        double vsum = 0.0;
+	int ncount = 0;
+
+        for (jj = 0; jj < jnum; jj++) {
+          j = jlist[jj];
+          j &= NEIGHMASK;
+
+          jtype = type[j];
+          delx = xtmp - x[j][0];
+          dely = ytmp - x[j][1];
+          delz = ztmp - x[j][2];
+          rsq = delx*delx + dely*dely + delz*delz;
+          if (rsq < cutsq && jtype >= typelo[0] && jtype <= typehi[0]) {
+	    double u, v;
+	    calc_q6(delx, dely, u, v);
+	    usum += u;
+	    vsum += v;
+	    ncount++;
+	  }
+	}
+	if (ncount > 0) {
+	  double ninv = 1.0/ncount ;
+	  q6[0] = usum*ninv;
+	  q6[1] = vsum*ninv;
+	}
+      }
+    }
+  } else {
+    for (ii = 0; ii < inum; ii++) {
+      i = ilist[ii];
+      double* q6 = q6array[i];
+      for (m = 0; m < ncol; m++) q6[m] = 0.0;
+
+      if (mask[i] & groupbit) {
+        xtmp = x[i][0];
+        ytmp = x[i][1];
+        ztmp = x[i][2];
+        jlist = firstneigh[i];
+        jnum = numneigh[i];
+
+	for (m = 0; m < ncol; m+=2) {
+
+	  double usum = 0.0;
+	  double vsum = 0.0;
+	  int ncount = 0;
+	
+	  for (jj = 0; jj < jnum; jj++) {
+	    j = jlist[jj];
+	    j &= NEIGHMASK;
+	
+	    jtype = type[j];
+	    delx = xtmp - x[j][0];
+	    dely = ytmp - x[j][1];
+	    delz = ztmp - x[j][2];
+	    rsq = delx*delx + dely*dely + delz*delz;
+	    if (rsq < cutsq) {
+	      if (jtype >= typelo[m] && jtype <= typehi[m]) {
+		double u, v;
+		calc_q6(delx, dely, u, v);
+		usum += u;
+		vsum += v;
+		ncount++;
+	      }
+	    }
+	    if (ncount > 0) {
+	      double ninv = 1.0/ncount ;
+	      q6[m] = usum*ninv;
+	      q6[m+1] = vsum*ninv;
+	    }
+	  }
+        }
+      }
+    }
+  }
+}
+
+inline void ComputeHexOrderAtom::calc_q6(double delx, double dely, double &u, double &v) {
+  double rinv = 1.0/sqrt(delx*delx+dely*dely);
+  double x = delx*rinv;
+  double y = dely*rinv;
+  double a = x*x;
+  double b1 = y*y;
+  double b2 = b1*b1;
+  double b3 = b2*b1;
+  u = ((  a - 15*b1)*a + 15*b2)*a - b3;
+  v = ((6*a - 20*b1)*a +  6*b2)*x*y;
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeHexOrderAtom::memory_usage()
+{
+  double bytes = ncol*nmax * sizeof(double);
+  return bytes;
+}
diff --git a/src/compute_hexorder_atom.h b/src/compute_hexorder_atom.h
new file mode 100644
index 0000000000..b9430addc2
--- /dev/null
+++ b/src/compute_hexorder_atom.h
@@ -0,0 +1,73 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(hexorder/atom,ComputeHexOrderAtom)
+
+#else
+
+#ifndef LMP_COMPUTE_HEXORDER_ATOM_H
+#define LMP_COMPUTE_HEXORDER_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeHexOrderAtom : public Compute {
+ public:
+  ComputeHexOrderAtom(class LAMMPS *, int, char **);
+  ~ComputeHexOrderAtom();
+  void init();
+  void init_list(int, class NeighList *);
+  void compute_peratom();
+  double memory_usage();
+
+ private:
+  int nmax,ncol;
+  double cutsq;
+  class NeighList *list;
+
+  int *typelo,*typehi;
+  double **q6array;
+
+  void calc_q6(double, double, double&, double&);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Compute hexorder/atom requires a pair style be defined
+
+Self-explantory.
+
+E: Compute hexorder/atom cutoff is longer than pairwise cutoff
+
+Cannot compute order parameter at distances longer than the pair cutoff,
+since those atoms are not in the neighbor list.
+
+W: More than one compute hexorder/atom
+
+It is not efficient to use compute hexorder/atom more than once.
+
+*/