From a91bbaf7f27615d6c6e20e7fd1f1d29fd83f944e Mon Sep 17 00:00:00 2001 From: athomps Date: Wed, 4 Nov 2015 23:56:47 +0000 Subject: [PATCH] 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 | 44 +++------- src/compute_hexorder_atom.cpp | 146 +++++++++------------------------- src/compute_hexorder_atom.h | 1 - 3 files changed, 47 insertions(+), 144 deletions(-) diff --git a/doc/compute_hexorder_atom.txt b/doc/compute_hexorder_atom.txt index 3f0d44150c..4b7e6c4109 100644 --- a/doc/compute_hexorder_atom.txt +++ b/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. diff --git a/src/compute_hexorder_atom.cpp b/src/compute_hexorder_atom.cpp index 1fedc33c23..9555bfc4ed 100644 --- a/src/compute_hexorder_atom.cpp +++ b/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; } } } diff --git a/src/compute_hexorder_atom.h b/src/compute_hexorder_atom.h index b9430addc2..99cd987161 100644 --- a/src/compute_hexorder_atom.h +++ b/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&);