Added hexatic bond orientational order parameter

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

doc/Eqs/hexorder.tex

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

doc/compute_hexorder_atom.txt

@@ -10,28 +10,27 @@ compute hexorder/atom command :h3
 
 [Syntax:]
 
-compute ID group-ID hexorder/atom cutoff type1 type2 ... :pre
+compute ID group-ID hexorder/atom :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)
 
 [Examples:]
 
-compute 1 all hexorder/atom 2.0 :pre
+compute 1 all hexorder/atom :pre
 
 [Description:]
 
 Define a computation that calculates {q}6 the hexatic bond-orientational 
 order parameter for each atom in a group. This order 
 parameter was introduced by "Nelson and Halperin"_#Nelson as a way to detect
-hexagonal symmetry in two-dimensional systems. For a each atoms, {q}6 
+hexagonal symmetry in two-dimensional systems. For each atom, {q}6 
 is a complex number (stored as two real numbers) defined as follows:
 
 :c,image(Eqs/hexorder.jpg)
 
-where the sum is over all atoms that are within 
+where the sum is over the six nearest neighbors 
-the specified cutoff distance from the central atom. The angle theta
+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
 central atom is calculated using all three  
@@ -46,9 +45,9 @@ 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
-specified compute group. An order parameter for atoms that have no 
+specified compute group. If the atom does not have 6 neighbors (within
-neighbors of the specified atom type within the cutoff distance will
+the potential cutoff), then its centro-symmetry parameter is set to
-be 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

src/compute_hexorder_atom.cpp

@@ -38,10 +38,7 @@ 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");
+  if (narg != 3) error->all(FLERR,"Illegal compute hexorder/atom command");
-
-  double cutoff = force->numeric(FLERR,arg[3]);
-  cutsq = cutoff*cutoff;
 
   ncol = 2;
 
@@ -50,6 +47,10 @@ ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) :
 
   nmax = 0;
   q6array = NULL;
+  maxneigh = 0;
+  distsq = NULL;
+  nearest = NULL;
+  nnn = 6;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -57,6 +58,8 @@ ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) :
 ComputeHexOrderAtom::~ComputeHexOrderAtom()
 {
   memory->destroy(q6array);
+  memory->destroy(distsq);
+  memory->destroy(nearest);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -65,9 +68,6 @@ 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
 
@@ -102,7 +102,7 @@ void ComputeHexOrderAtom::compute_peratom()
 
   invoked_peratom = update->ntimestep;
 
-  // grow coordination array if necessary
+  // grow order parameter array if necessary
 
   if (atom->nlocal > nmax) {
     memory->destroy(q6array);
@@ -120,16 +120,16 @@ void ComputeHexOrderAtom::compute_peratom()
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  // compute order parameter(s) for each atom in group
+  // compute order parameter for each atom in group
   // use full neighbor list to count atoms less than cutoff
 
   double **x = atom->x;
   int *mask = atom->mask;
+  double cutsq = force->pair->cutforce * force->pair->cutforce;
 
   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];
@@ -137,31 +137,62 @@ void ComputeHexOrderAtom::compute_peratom()
       jlist = firstneigh[i];
       jnum = numneigh[i];
       
+      // insure distsq and nearest arrays are long enough
+
+      if (jnum > maxneigh) {
+        memory->destroy(distsq);
+        memory->destroy(nearest);
+        maxneigh = jnum;
+        memory->create(distsq,maxneigh,"hexcoord/atom:distsq");
+        memory->create(nearest,maxneigh,"hexcoord/atom:nearest");
+      }
+
+      // loop over list of all neighbors within force cutoff
+      // distsq[] = distance sq to each
+      // nearest[] = atom indices of neighbors
+
+      int ncount = 0;
+      for (jj = 0; jj < jnum; jj++) {
+        j = jlist[jj];
+        j &= NEIGHMASK;
+
+        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) {
+          distsq[ncount] = rsq;
+          nearest[ncount++] = j;
+        }
+      }
+
+      // if not nnn neighbors, order parameter = 0;
+
+      if (ncount < nnn) {
+	q6[0] = q6[1] = 0.0;
+        continue;
+      }
+
+      // store nnn nearest neighs in 1st nnn locations of distsq and nearest
+
+      select2(nnn,ncount,distsq,nearest);
+
       double usum = 0.0;
       double vsum = 0.0;
-      int ncount = 0;
       
-      for (jj = 0; jj < jnum; jj++) {
+      for (jj = 0; jj < nnn; jj++) {
-	j = jlist[jj];
+	j = nearest[jj];
 	j &= NEIGHMASK;
 	
 	delx = xtmp - x[j][0];
 	dely = ytmp - x[j][1];
-	delz = ztmp - x[j][2];
+	double u, v;
-	rsq = delx*delx + dely*dely + delz*delz;
+	calc_q6(delx, dely, u, v);
-	if (rsq < cutsq) {
+	usum += u;
-	  double u, v;
+	vsum += 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;
       }
+      q6[0] = usum/nnn;
+      q6[1] = vsum/nnn;
     }
   }
 }
@@ -178,6 +209,70 @@ inline void ComputeHexOrderAtom::calc_q6(double delx, double dely, double &u, do
   v = ((6*a - 20*b1)*a +  6*b2)*x*y;
 }
 
+/* ----------------------------------------------------------------------
+   select2 routine from Numerical Recipes (slightly modified)
+   find k smallest values in array of length n
+   sort auxiliary array at same time
+------------------------------------------------------------------------- */
+
+#define SWAP(a,b)   tmp = a; a = b; b = tmp;
+#define ISWAP(a,b) itmp = a; a = b; b = itmp;
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHexOrderAtom::select2(int k, int n, double *arr, int *iarr)
+{
+  int i,ir,j,l,mid,ia,itmp;
+  double a,tmp;
+
+  arr--;
+  iarr--;
+  l = 1;
+  ir = n;
+  for (;;) {
+    if (ir <= l+1) {
+      if (ir == l+1 && arr[ir] < arr[l]) {
+        SWAP(arr[l],arr[ir])
+        ISWAP(iarr[l],iarr[ir])
+      }
+      return;
+    } else {
+      mid=(l+ir) >> 1;
+      SWAP(arr[mid],arr[l+1])
+      ISWAP(iarr[mid],iarr[l+1])
+      if (arr[l] > arr[ir]) {
+        SWAP(arr[l],arr[ir])
+        ISWAP(iarr[l],iarr[ir])
+      }
+      if (arr[l+1] > arr[ir]) {
+        SWAP(arr[l+1],arr[ir])
+        ISWAP(iarr[l+1],iarr[ir])
+      }
+      if (arr[l] > arr[l+1]) {
+        SWAP(arr[l],arr[l+1])
+        ISWAP(iarr[l],iarr[l+1])
+      }
+      i = l+1;
+      j = ir;
+      a = arr[l+1];
+      ia = iarr[l+1];
+      for (;;) {
+        do i++; while (arr[i] < a);
+        do j--; while (arr[j] > a);
+        if (j < i) break;
+        SWAP(arr[i],arr[j])
+        ISWAP(iarr[i],iarr[j])
+      }
+      arr[l+1] = arr[j];
+      arr[j] = a;
+      iarr[l+1] = iarr[j];
+      iarr[j] = ia;
+      if (j >= k) ir = j-1;
+      if (j <= k) l = i;
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    memory usage of local atom-based array
 ------------------------------------------------------------------------- */

src/compute_hexorder_atom.h

@@ -34,13 +34,15 @@ class ComputeHexOrderAtom : public Compute {
   double memory_usage();
 
  private:
-  int nmax,ncol;
+  int nmax,maxneigh,ncol,nnn;
-  double cutsq;
   class NeighList *list;
+  double *distsq;
+  int *nearest;
 
   double **q6array;
 
   void calc_q6(double, double, double&, double&);
+  void select2(int, int, double *, int *);
 };
 
 }
@@ -60,11 +62,6 @@ 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.