diff --git a/doc/src/compute_coord_atom.txt b/doc/src/compute_coord_atom.txt index 012a87a9a7..f1a6bf7ffb 100644 --- a/doc/src/compute_coord_atom.txt +++ b/doc/src/compute_coord_atom.txt @@ -10,22 +10,34 @@ compute coord/atom command :h3 [Syntax:] -compute ID group-ID coord/atom cutoff type1 type2 ... :pre +compute ID group-ID coord/atom cstyle args ... :pre ID, group-ID are documented in "compute"_compute.html command coord/atom = style name of this compute command -cutoff = distance within which to count coordination neighbors (distance units) -typeN = atom type for Nth coordination count (see asterisk form below) :ul +one cstyle must be appended :ul + +cstyle = {cutoff} or {orientorder} + +{cutoff} args = cutoff typeN + cutoff = distance within which to count coordination neighbors (distance units) + typeN = atom type for Nth coordination count (see asterisk form below) :pre + +{orientorder} args = orientorderID threshold + orientorderID = ID of a previously defined orientorder/atom compute + threshold = minimum value of the scalar product between two 'connected' atoms (see text for explanation) :pre [Examples:] -compute 1 all coord/atom 2.0 -compute 1 all coord/atom 6.0 1 2 -compute 1 all coord/atom 6.0 2*4 5*8 * :pre +compute 1 all coord/atom cutoff 2.0 +compute 1 all coord/atom cutoff 6.0 1 2 +compute 1 all coord/atom cutoff 6.0 2*4 5*8 * +compute 1 all coord/atom orientorder 2 0.5 :pre [Description:] -Define a computation that calculates one or more coordination numbers +This compute performs generic calculations between neighboring atoms. So far, +there are two cstyles implemented: {cutoff} and {orientorder}. +The {cutoff} cstyle calculates one or more coordination numbers for each atom in a group. A coordination number is defined as the number of neighbor atoms with @@ -49,6 +61,14 @@ from 1 to N. A leading asterisk means all types from 1 to n (inclusive). A middle asterisk means all types from m to n (inclusive). +The {orientorder} cstyle calculates the number of 'connected' atoms j +around each atom i. The atom j is connected to i if the scalar product +({Ybar_lm(i)},{Ybar_lm(j)}) is larger than {threshold}. Thus, this cstyle +will work only if a "compute orientorder/atom"_compute_orientorder_atom.html +has been previously defined. This cstyle allows one to apply the +ten Wolde's criterion to identify cristal-like atoms in a system +(see "ten Wolde et al."_#tenWolde). + The value of all coordination numbers will be 0.0 for atoms not in the specified compute group. @@ -83,10 +103,19 @@ options. The per-atom vector or array values will be a number >= 0.0, as explained above. -[Restrictions:] none +[Restrictions:] +The cstyle {orientorder} can only be used if a +"compute orientorder/atom"_compute_orientorder_atom.html command +was previously defined. Otherwise, an error message will be issued. [Related commands:] "compute cluster/atom"_compute_cluster_atom.html +"compute orientorder/atom"_compute_orientorder_atom.html [Default:] none + +:line + +:link(tenWolde) +[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, J. Chem. Phys. 104, 9932 (1996). diff --git a/doc/src/compute_orientorder_atom.txt b/doc/src/compute_orientorder_atom.txt index c5ecef3cb2..74426dd33b 100644 --- a/doc/src/compute_orientorder_atom.txt +++ b/doc/src/compute_orientorder_atom.txt @@ -15,17 +15,19 @@ compute ID group-ID orientorder/atom keyword values ... :pre ID, group-ID are documented in "compute"_compute.html command :ulb,l orientorder/atom = style name of this compute command :l one or more keyword/value pairs may be appended :l -keyword = {cutoff} or {nnn} or {degrees} +keyword = {cutoff} or {nnn} or {degrees} or {components} {cutoff} value = distance cutoff {nnn} value = number of nearest neighbors - {degrees} values = nlvalues, l1, l2,... :pre + {degrees} values = nlvalues, l1, l2,... + {components} value = l :pre :ule [Examples:] compute 1 all orientorder/atom -compute 1 all orientorder/atom degrees 5 4 6 8 10 12 nnn NULL cutoff 1.5 :pre +compute 1 all orientorder/atom degrees 5 4 6 8 10 12 nnn NULL cutoff 1.5 +compute 1 all orientorder/atom degrees 4 6 components 6 nnn NULL cutoff 3.0 :pre [Description:] @@ -71,6 +73,13 @@ The numerical values of all order parameters up to {Q}12 for a range of commonly encountered high-symmetry structures are given in Table I of "Mickel et al."_#Mickel. +The optional keyword {components} will output the components of +the normalized complex vector {Ybar_lm} of degree {l}, which must be +explicitly included in the keyword {degrees}. This option can be used +in conjunction with "compute coord_atom"_compute_coord_atom.html to +calculate the ten Wolde's criterion to identify crystal-like particles +(see "ten Wolde et al."_#tenWolde96). + The value of {Ql} is set to zero for atoms not in the specified compute group, as well as for atoms that have less than {nnn} neighbors within the distance cutoff. @@ -98,6 +107,12 @@ the neighbor list. This compute calculates a per-atom array with {nlvalues} columns, giving the {Ql} values for each atom, which are real numbers on the range 0 <= {Ql} <= 1. +If the keyword {components} is set, then the real and imaginary parts of each +component of (normalized) {Ybar_lm} will be added to the output array in the +following order: +Re({Ybar_-m}) Im({Ybar_-m}) Re({Ybar_-m+1}) Im({Ybar_-m+1}) ... Re({Ybar_m}) Im({Ybar_m}). +This way, the per-atom array will have a total of {nlvalues}+2*(2{l}+1) columns. + These values can be accessed by any command that uses per-atom values from a compute as input. See "Section 6.15"_Section_howto.html#howto_15 for an overview of LAMMPS output @@ -117,5 +132,9 @@ The option defaults are {cutoff} = pair style cutoff, {nnn} = 12, {degrees} = 5 :link(Steinhardt) [(Steinhardt)] P. Steinhardt, D. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983). + :link(Mickel) [(Mickel)] W. Mickel, S. C. Kapfer, G. E. Schroeder-Turkand, K. Mecke, J. Chem. Phys. 138, 044501 (2013). + +:link(tenWolde96) +[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, J. Chem. Phys. 104, 9932 (1996). diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp index 21744dcc93..1ad553ee37 100644 --- a/src/compute_coord_atom.cpp +++ b/src/compute_coord_atom.cpp @@ -15,6 +15,7 @@ #include #include #include "compute_coord_atom.h" +#include "compute_orientorder_atom.h" #include "atom.h" #include "update.h" #include "modify.h" @@ -29,37 +30,72 @@ using namespace LAMMPS_NS; +#define INVOKED_PERATOM 8 + /* ---------------------------------------------------------------------- */ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), - typelo(NULL), typehi(NULL), cvec(NULL), carray(NULL) + cstyle(NULL), id_orientorder(NULL), typelo(NULL), typehi(NULL), cvec(NULL), carray(NULL) { - if (narg < 4) error->all(FLERR,"Illegal compute coord/atom command"); + if (narg < 5) error->all(FLERR,"Illegal compute coord/atom command"); - double cutoff = force->numeric(FLERR,arg[3]); - cutsq = cutoff*cutoff; + int n = strlen(arg[3]) + 1; + cstyle = new char[n]; + strcpy(cstyle,arg[3]); - ncol = narg-4 + 1; - int ntypes = atom->ntypes; - typelo = new int[ncol]; - typehi = new int[ncol]; + if (strcmp(cstyle,"cutoff") == 0) { + + double cutoff = force->numeric(FLERR,arg[4]); + cutsq = cutoff*cutoff; + + ncol = narg-5 + 1; + int ntypes = atom->ntypes; + typelo = new int[ncol]; + typehi = new int[ncol]; + + if (narg == 5) { + ncol = 1; + typelo[0] = 1; + typehi[0] = ntypes; + } else { + ncol = 0; + int iarg = 5; + while (iarg < narg) { + force->bounds(FLERR,arg[iarg],ntypes,typelo[ncol],typehi[ncol]); + if (typelo[ncol] > typehi[ncol]) + error->all(FLERR,"Illegal compute coord/atom command"); + ncol++; + iarg++; + } - if (narg == 4) { - ncol = 1; - typelo[0] = 1; - typehi[0] = ntypes; - } else { - ncol = 0; - int iarg = 4; - while (iarg < narg) { - force->bounds(FLERR,arg[iarg],ntypes,typelo[ncol],typehi[ncol]); - if (typelo[ncol] > typehi[ncol]) - error->all(FLERR,"Illegal compute coord/atom command"); - ncol++; - iarg++; } - } + + } else if (strcmp(cstyle,"orientorder") == 0) { + + if (narg != 6) error->all(FLERR,"Illegal compute coord/atom command"); + + n = strlen(arg[4]) + 1; + id_orientorder = new char[n]; + strcpy(id_orientorder,arg[4]); + + int iorientorder = modify->find_compute(id_orientorder); + if (iorientorder < 0) + error->all(FLERR,"Could not find compute coord/atom compute ID"); + if (strcmp(modify->compute[iorientorder]->style,"orientorder/atom") != 0) + error->all(FLERR,"Compute coord/atom compute ID does not compute orientorder/atom"); + + threshold = force->numeric(FLERR,arg[5]); + if (threshold <= -1.0 || threshold >= 1.0) + error->all(FLERR,"Compute coord/atom threshold value must lie between -1 and 1"); + + ncol = 1; + typelo = new int[ncol]; + typehi = new int[ncol]; + typelo[0] = 1; + typehi[0] = atom->ntypes; + + } else error->all(FLERR,"Invalid cstyle in compute coord/atom"); peratom_flag = 1; if (ncol == 1) size_peratom_cols = 0; @@ -82,6 +118,17 @@ ComputeCoordAtom::~ComputeCoordAtom() void ComputeCoordAtom::init() { + if (strcmp(cstyle,"orientorder") == 0) { + int iorientorder = modify->find_compute(id_orientorder); + c_orientorder = (ComputeOrientOrderAtom*)(modify->compute[iorientorder]); + cutsq = c_orientorder->cutsq; + l = c_orientorder->qlcomp; + // communicate real and imaginary 2*l+1 components of the normalized vector + comm_forward = 2*(2*l+1); + if (c_orientorder->iqlcomp < 0) + error->all(FLERR,"Compute coord/atom requires components option in compute orientorder/atom be defined"); + } + if (force->pair == NULL) error->all(FLERR,"Compute coord/atom requires a pair style be defined"); if (sqrt(cutsq) > force->pair->cutforce) @@ -122,6 +169,9 @@ void ComputeCoordAtom::compute_peratom() invoked_peratom = update->ntimestep; +// printf("Number of degrees %i components degree %i",nqlist,l); +// printf("Particle \t %i \t Norm \t %g \n",0,norm[0][0]); + // grow coordination array if necessary if (atom->nmax > nmax) { @@ -138,6 +188,19 @@ void ComputeCoordAtom::compute_peratom() } } + if (strcmp(cstyle,"orientorder") == 0) { + if (!(c_orientorder->invoked_flag & INVOKED_PERATOM)) { + c_orientorder->compute_peratom(); + c_orientorder->invoked_flag |= INVOKED_PERATOM; + } + nqlist = c_orientorder->nqlist; + int ltmp = l; +// l = c_orientorder->qlcomp; + if (ltmp != l) error->all(FLERR,"Debug error, ltmp != l\n"); + normv = c_orientorder->array_atom; + comm->forward_comm_compute(this); + } + // invoke full neighbor list (will copy or build if necessary) neighbor->build_one(list); @@ -154,65 +217,131 @@ void ComputeCoordAtom::compute_peratom() int *type = atom->type; int *mask = atom->mask; - if (ncol == 1) { - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - if (mask[i] & groupbit) { - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - jlist = firstneigh[i]; - jnum = numneigh[i]; + if (strcmp(cstyle,"cutoff") == 0) { - n = 0; - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; + if (ncol == 1) { - 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]) n++; - } + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (mask[i] & groupbit) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - cvec[i] = n; - } else cvec[i] = 0.0; - } + n = 0; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; - } else { - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - count = carray[i]; - for (m = 0; m < ncol; m++) count[m] = 0.0; + 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]) + n++; + } - if (mask[i] & groupbit) { - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - jlist = firstneigh[i]; - jnum = numneigh[i]; + cvec[i] = n; + } else cvec[i] = 0.0; + } + + } else { + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + count = carray[i]; + for (m = 0; m < ncol; m++) count[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 (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; + 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) { - for (m = 0; m < ncol; m++) - if (jtype >= typelo[m] && jtype <= typehi[m]) - count[m] += 1.0; + 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) { + for (m = 0; m < ncol; m++) + if (jtype >= typelo[m] && jtype <= typehi[m]) + count[m] += 1.0; + } } } } } + + } else if (strcmp(cstyle,"orientorder") == 0) { + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (mask[i] & groupbit) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + n = 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 dot_product = 0.0; + for (int m=0; m < 2*(2*l+1); m++) { + dot_product += normv[i][nqlist+m]*normv[j][nqlist+m]; + } + if (dot_product > threshold) n++; + } + } + cvec[i] = n; + } else cvec[i] = 0.0; } + } +} + +/* ---------------------------------------------------------------------- */ + +int ComputeCoordAtom::pack_forward_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,m=0,j; + for (i = 0; i < n; ++i) { + for (j = nqlist; j < nqlist + 2*(2*l+1); ++j) { + buf[m++] = normv[list[i]][j]; + } + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeCoordAtom::unpack_forward_comm(int n, int first, double *buf) +{ + int i,last,m=0,j; + last = first + n; + for (i = first; i < last; ++i) { + for (j = nqlist; j < nqlist + 2*(2*l+1); ++j) { + normv[i][j] = buf[m++]; + } + } + } /* ---------------------------------------------------------------------- diff --git a/src/compute_coord_atom.h b/src/compute_coord_atom.h index 0ff373f133..942882179f 100644 --- a/src/compute_coord_atom.h +++ b/src/compute_coord_atom.h @@ -31,6 +31,8 @@ class ComputeCoordAtom : public Compute { void init(); void init_list(int, class NeighList *); void compute_peratom(); + int pack_forward_comm(int, int *, double *, int, int *); + void unpack_forward_comm(int, int, double *); double memory_usage(); private: @@ -41,6 +43,12 @@ class ComputeCoordAtom : public Compute { int *typelo,*typehi; double *cvec; double **carray; + + class ComputeOrientOrderAtom *c_orientorder; + char *cstyle,*id_orientorder; + double threshold; + double **normv; + int nqlist,l; }; } diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp index 6c5a2c0c0e..42186f03b8 100644 --- a/src/compute_orientorder_atom.cpp +++ b/src/compute_orientorder_atom.cpp @@ -54,6 +54,7 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg nnn = 12; cutsq = 0.0; + qlcompflag = 0; // specify which orders to request @@ -96,6 +97,20 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg if (qlist[iw] > qmax) qmax = qlist[iw]; } iarg += nqlist; + if (strcmp(arg[iarg],"components") == 0) { + qlcompflag = 1; + if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); + qlcomp = force->numeric(FLERR,arg[iarg+1]); + if (qlcomp <= 0) error->all(FLERR,"Illegal compute orientorder/atom command"); + iqlcomp = -1; + for (int iw = 0; iw < nqlist; iw++) + if (qlcomp == qlist[iw]) { + iqlcomp = iw; + break; + } + if (iqlcomp < 0) error->all(FLERR,"Illegal compute orientorder/atom command"); + iarg += 2; + } } else if (strcmp(arg[iarg],"cutoff") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); double cutoff = force->numeric(FLERR,arg[iarg+1]); @@ -105,7 +120,9 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg } else error->all(FLERR,"Illegal compute orientorder/atom command"); } - ncol = nqlist; + if (qlcompflag) ncol = nqlist + 2*(2*qlcomp+1); + else ncol = nqlist; + peratom_flag = 1; size_peratom_cols = ncol; @@ -434,6 +451,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, } double fac = sqrt(MY_4PI) / ncount; + double normfac = 0.0; for (int iw = 0; iw < nqlist; iw++) { int n = qlist[iw]; double qm_sum = 0.0; @@ -443,6 +461,18 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, // qnm_r[iw][m]*qnm_r[iw][m] + qnm_i[iw][m]*qnm_i[iw][m]); } qn[iw] = fac * sqrt(qm_sum / (2*n+1)); + if (qlcompflag && iqlcomp == iw) normfac = 1.0/sqrt(qm_sum); + + } + + // output of the complex vector + + if (qlcompflag) { + int j = nqlist; + for(int m = 0; m < 2*qlcomp+1; m++) { + qn[j++] = qnm_r[iqlcomp][m] * normfac; + qn[j++] = qnm_i[iqlcomp][m] * normfac; + } } } diff --git a/src/compute_orientorder_atom.h b/src/compute_orientorder_atom.h index 9c5ec14f56..d6f32b7727 100644 --- a/src/compute_orientorder_atom.h +++ b/src/compute_orientorder_atom.h @@ -32,6 +32,10 @@ class ComputeOrientOrderAtom : public Compute { void init_list(int, class NeighList *); void compute_peratom(); double memory_usage(); + double cutsq; + int iqlcomp, qlcomp, qlcompflag; + int *qlist; + int nqlist; private: int nmax,maxneigh,ncol,nnn; @@ -39,11 +43,8 @@ class ComputeOrientOrderAtom : public Compute { double *distsq; int *nearest; double **rlist; - int *qlist; - int nqlist; int qmax; double **qnarray; - double cutsq; double **qnm_r; double **qnm_i;