import contributed code for computes coord/atom and orientorder/atom

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
+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) :ul
+  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 cutoff 2.0
-compute 1 all coord/atom 6.0 1 2
+compute 1 all coord/atom cutoff 6.0 1 2
-compute 1 all coord/atom 6.0 2*4 5*8 * :pre
+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).

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).

src/compute_coord_atom.cpp

@@ -15,6 +15,7 @@
 #include <string.h>
 #include <stdlib.h>
 #include "compute_coord_atom.h"
+#include "compute_orientorder_atom.h"
 #include "atom.h"
 #include "update.h"
 #include "modify.h"
@@ -29,29 +30,37 @@
 
 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]);
+  int n = strlen(arg[3]) + 1;
+  cstyle = new char[n];
+  strcpy(cstyle,arg[3]);
+
+  if (strcmp(cstyle,"cutoff") == 0) {
+
+    double cutoff = force->numeric(FLERR,arg[4]);
     cutsq = cutoff*cutoff;
 
-  ncol = narg-4 + 1;
+    ncol = narg-5 + 1;
     int ntypes = atom->ntypes;
     typelo = new int[ncol];
     typehi = new int[ncol];
 
-  if (narg == 4) {
+    if (narg == 5) {
       ncol = 1;
       typelo[0] = 1;
       typehi[0] = ntypes;
     } else {
       ncol = 0;
-    int iarg = 4;
+      int iarg = 5;
       while (iarg < narg) {
         force->bounds(FLERR,arg[iarg],ntypes,typelo[ncol],typehi[ncol]);
         if (typelo[ncol] > typehi[ncol])
@@ -59,8 +68,35 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
         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;
   else size_peratom_cols = ncol;
@@ -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,7 +217,10 @@ void ComputeCoordAtom::compute_peratom()
   int *type = atom->type;
   int *mask = atom->mask;
 
+  if (strcmp(cstyle,"cutoff") == 0) {
+
     if (ncol == 1) {
+
       for (ii = 0; ii < inum; ii++) {
         i = ilist[ii];
         if (mask[i] & groupbit) {
@@ -174,7 +240,8 @@ void ComputeCoordAtom::compute_peratom()
             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 (rsq < cutsq && jtype >= typelo[0] && jtype <= typehi[0])
+	      n++;
           }
 
           cvec[i] = n;
@@ -213,6 +280,68 @@ void ComputeCoordAtom::compute_peratom()
         }
       }
     }
+
+  } 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++];
+    }
+  }
+
 }
 
 /* ----------------------------------------------------------------------

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;
 };
 
 }

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;
+    }
   }
 }
 

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;