You could also ask the authors of that paper via e-mail how they did it, they might even share their input or post-processing script(s) if you’re lucky.

As for Axel’s suggestion, it will probably become clearer if you read the documentation for fix ave/correlate and the list of computes LAMMPS offers.

In order to calculate orientational correlation, you first need q6,
the hexatic orientational order parameter. This is now a fairly
standard quantity in condensed matter physics. It has probably been
implemented at least once in LAMMPS, as it is a straight forward
per-atom quantity similar to the centro-symmetry parameter, but
complex not real. Unfortunately, no one has yet shared their code with
us, so you will have to code it yourself. Starting with a compute like
coord/atom or centro/atom, modify it to calculate exp(6*i*theta) =
(x+i y)^6, where x and y are coordinates of the unit vector along the
bond from central atom i to neighbor atom j. Average that over all
neighbors of the central atom and you are done.