Added hexatic bond orientational order parameter

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

doc/compute_hexorder_atom.txt

@@ -15,53 +15,35 @@ 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
+compute 1 all hexorder/atom 2.0 :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:
+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 
+is a complex number (stored as two real numbers) defined as follows:
 
 :c,image(Eqs/hexorder.jpg)
 
-where the sum is over atoms of the specified atom type(s) that are within 
+where the sum is over all atoms 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 
+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, 
 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).
+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 
@@ -88,19 +70,15 @@ 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
+This compute calculates a per-atom array with 2 columns, giving the
 real and imaginary parts of {q}6, 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 {q}6.  
 
 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 pairs of numbers representing the
+The per-atom array contain pairs of numbers representing the
 real and imaginary parts of {q}6, a complex number subject to the
 constraint |{q}6| <= 1.
 

src/compute_hexorder_atom.cpp

@@ -38,33 +38,12 @@ 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 != 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++;
-    }
-  }
+  ncol = 2;
 
   peratom_flag = 1;
   size_peratom_cols = ncol;
@@ -77,8 +56,6 @@ ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) :
 
 ComputeHexOrderAtom::~ComputeHexOrderAtom()
 {
-  delete [] typelo;
-  delete [] typehi;
   memory->destroy(q6array);
 }
 
@@ -119,10 +96,9 @@ void ComputeHexOrderAtom::init_list(int id, NeighList *ptr)
 
 void ComputeHexOrderAtom::compute_peratom()
 {
-  int i,j,m,ii,jj,inum,jnum,jtype;
+  int i,j,m,ii,jj,inum,jnum;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  double *count;
 
   invoked_peratom = update->ntimestep;
 
@@ -148,93 +124,43 @@ void ComputeHexOrderAtom::compute_peratom()
   // 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;
+  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;
+	
+	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) {
+	  double u, v;
+	  calc_q6(delx, dely, u, v);
+	  usum += u;
+	  vsum += v;
+	  ncount++;
 	}
       }
-    }
-  } 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;
-	    }
-	  }
-        }
+      if (ncount > 0) {
+	double ninv = 1.0/ncount ;
+	q6[0] = usum*ninv;
+	q6[1] = vsum*ninv;
       }
     }
   }

src/compute_hexorder_atom.h

@@ -38,7 +38,6 @@ class ComputeHexOrderAtom : public Compute {
   double cutsq;
   class NeighList *list;
 
-  int *typelo,*typehi;
   double **q6array;
 
   void calc_q6(double, double, double&, double&);