Added hexatic bond orientational order parameter

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14236 f3b2605a-c512-4ea7-a41b-209d697bcdaa

doc/Eqs/hexorder.tex

@@ -2,7 +2,7 @@
 \begin{document}
 
 $$
-   q_6 = \frac{1}{6}\sum_{j = 1}^{6} e^{6 i \theta({\bf r}_{ij})}
+   q_n = \frac{1}{n}\sum_{j = 1}^{n} e^{n i \theta({\bf r}_{ij})}
 $$
 
 \end{document}

doc/compute_hexorder_atom.txt

@@ -10,26 +10,31 @@ compute hexorder/atom command :h3
 
 [Syntax:]
 
-compute ID group-ID hexorder/atom :pre
+compute ID group-ID hexorder/atom keyword values ... :pre
 
-ID, group-ID are documented in "compute"_compute.html command
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
-hexorder/atom = style name of this compute command
+hexorder/atom = style name of this compute command :l
+zero or more keyword/value pairs may be appended :l
+keyword = {degree} :l
+  {n} value = degree of order parameter :pre
+:ule
 
 [Examples:]
 
-compute 1 all hexorder/atom :pre
+compute 1 all hexorder/atom 
+compute 1 all hexorder/atom n 4 :pre
 
 [Description:]
 
-Define a computation that calculates {q}6 the hexatic bond-orientational 
+Define a computation that calculates {qn} the bond-orientational 
-order parameter for each atom in a group. This order 
+order parameter for each atom in a group. The hexatic ({n} = 6) order 
 parameter was introduced by "Nelson and Halperin"_#Nelson as a way to detect
-hexagonal symmetry in two-dimensional systems. For each atom, {q}6 
+hexagonal symmetry in two-dimensional systems. For each atom, {qn} 
 is a complex number (stored as two real numbers) defined as follows:
 
 :c,image(Eqs/hexorder.jpg)
 
-where the sum is over the six nearest neighbors 
+where the sum is over the {n} nearest neighbors 
 of 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
@@ -38,16 +43,16 @@ central atom is calculated using all three
 Neighbor 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, 
+The optional keyword {n} sets the degree of the order parameter. 
-then |{q}6| = 1.  
+The default value is 6. 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 symmetry, |{q}6| << 1.
 
-The value of all order parameters will be zero for atoms not in the
+The value of {qn} will be zero for atoms not in the
-specified compute group. If the atom does not have 6 neighbors (within
+specified compute group. If the atom has less than {n} neighbors (within
-the potential cutoff), then its centro-symmetry parameter is set to
+the potential cutoff), then {qn} is set to zero.
-zero.
 
 The neighbor list needed to compute this quantity is constructed each
 time the calculation is performed (i.e. each time a snapshot of atoms
@@ -70,7 +75,7 @@ the neighbor list.
 [Output info:]
 
 This compute calculates a per-atom array with 2 columns, giving the
-real and imaginary parts of {q}6, respectively.  
+real and imaginary parts of {qn}, respectively.  
 
 These values can be accessed by any command that uses
 per-atom values from a compute as input.  See "Section_howto
@@ -78,8 +83,8 @@ per-atom values from a compute as input.  See "Section_howto
 options.
 
 The per-atom array contain pairs of numbers representing the
-real and imaginary parts of {q}6, a complex number subject to the
+real and imaginary parts of {qn}, a complex number subject to the
-constraint |{q}6| <= 1.
+constraint |{qn}| <= 1.
 
 [Restrictions:] none
 
@@ -87,7 +92,9 @@ constraint |{q}6| <= 1.
 
 "compute coord/atom"_compute_coord_atom.html, "compute centro/atom"_compute_centro_atom.html
 
-[Default:] none
+[Default:] 
+
+The option default is {n} = 6.
 
 :line
 

src/compute_hexorder_atom.cpp

@@ -16,6 +16,7 @@
 ------------------------------------------------------------------------- */
 
 #include <math.h>
+#include <complex>
 #include <string.h>
 #include <stdlib.h>
 #include "compute_hexorder_atom.h"
@@ -38,10 +39,24 @@ using namespace LAMMPS_NS;
 ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg)
 {
-  if (narg != 3) error->all(FLERR,"Illegal compute hexorder/atom command");
+  if (narg < 3 ) error->all(FLERR,"Illegal compute hexorder/atom command");
+
+  nnn = 6;
+
+  // process optional args
+
+  int iarg = 3;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"degree") == 0) {
+      if (iarg+1 > narg) error->all(FLERR,"Illegal lattice command");
+      nnn = force->numeric(FLERR,arg[iarg+1]);
+      if (nnn < 0)
+        error->all(FLERR,"Illegal lattice command");
+      iarg += 2;
+    }
+  }
 
   ncol = 2;
-
   peratom_flag = 1;
   size_peratom_cols = ncol;
 
@@ -50,7 +65,6 @@ ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) :
   maxneigh = 0;
   distsq = NULL;
   nearest = NULL;
-  nnn = 6;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -187,7 +201,7 @@ void ComputeHexOrderAtom::compute_peratom()
 	delx = xtmp - x[j][0];
 	dely = ytmp - x[j][1];
 	double u, v;
-	calc_q6(delx, dely, u, v);
+	calc_qn(delx, dely, u, v);
 	usum += u;
 	vsum += v;
       }
@@ -197,6 +211,8 @@ void ComputeHexOrderAtom::compute_peratom()
   }
 }
 
+// this might be faster than pow(std::complex) on some platforms
+
 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;
@@ -205,10 +221,23 @@ inline void ComputeHexOrderAtom::calc_q6(double delx, double dely, double &u, do
   double b1 = y*y;
   double b2 = b1*b1;
   double b3 = b2*b1;
+
+  // (x + i y)^6 coeffs: 1, 6, -15, -20, 15, 6, -1
+
   u = ((  a - 15*b1)*a + 15*b2)*a - b3;
   v = ((6*a - 20*b1)*a +  6*b2)*x*y;
 }
 
+inline void ComputeHexOrderAtom::calc_qn(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;
+  std::complex<double> z = x + y*1i;
+  std::complex<double> zn = pow(z,nnn);
+  u = real(zn);
+  v = imag(zn);
+}
+
 /* ----------------------------------------------------------------------
    select2 routine from Numerical Recipes (slightly modified)
    find k smallest values in array of length n

src/compute_hexorder_atom.h

@@ -42,6 +42,8 @@ class ComputeHexOrderAtom : public Compute {
   double **q6array;
 
   void calc_q6(double, double, double&, double&);
+  void calc_q4(double, double, double&, double&);
+  void calc_qn(double, double, double&, double&);
   void select2(int, int, double *, int *);
 };