diff --git a/doc/Eqs/hexorder.jpg b/doc/Eqs/hexorder.jpg index cebb8d9d9e..466d6b8090 100644 Binary files a/doc/Eqs/hexorder.jpg and b/doc/Eqs/hexorder.jpg differ diff --git a/doc/Eqs/hexorder.tex b/doc/Eqs/hexorder.tex index 09fd4e2706..fd7dac717b 100644 --- a/doc/Eqs/hexorder.tex +++ b/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} diff --git a/doc/compute_hexorder_atom.txt b/doc/compute_hexorder_atom.txt index 4b7e6c4109..90df71d749 100644 --- a/doc/compute_hexorder_atom.txt +++ b/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 -the specified cutoff distance from the central atom. The angle theta +where the sum is over the six 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 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 -neighbors of the specified atom type within the cutoff distance will -be zero. +specified compute group. If the atom does not have 6 neighbors (within +the potential cutoff), then its centro-symmetry parameter is set to +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 diff --git a/src/compute_hexorder_atom.cpp b/src/compute_hexorder_atom.cpp index 9555bfc4ed..54b018512c 100644 --- a/src/compute_hexorder_atom.cpp +++ b/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"); - - double cutoff = force->numeric(FLERR,arg[3]); - cutsq = cutoff*cutoff; + if (narg != 3) error->all(FLERR,"Illegal compute hexorder/atom command"); 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++) { - j = jlist[jj]; + for (jj = 0; jj < nnn; jj++) { + j = nearest[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++; - } - } - if (ncount > 0) { - double ninv = 1.0/ncount ; - q6[0] = usum*ninv; - q6[1] = vsum*ninv; + double u, v; + calc_q6(delx, dely, u, v); + usum += u; + vsum += v; } + 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 ------------------------------------------------------------------------- */ diff --git a/src/compute_hexorder_atom.h b/src/compute_hexorder_atom.h index 99cd987161..38f34cab89 100644 --- a/src/compute_hexorder_atom.h +++ b/src/compute_hexorder_atom.h @@ -34,13 +34,15 @@ class ComputeHexOrderAtom : public Compute { double memory_usage(); private: - int nmax,ncol; - double cutsq; + int nmax,maxneigh,ncol,nnn; 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.