From 95706ac846fba3fdf2a84e95a406e0b96dd0dd17 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Tue, 10 Jan 2017 12:29:22 -0500
Subject: [PATCH] import contributed code for computes coord/atom and
orientorder/atom
---
doc/src/compute_coord_atom.txt | 45 ++++-
doc/src/compute_orientorder_atom.txt | 25 ++-
src/compute_coord_atom.cpp | 283 +++++++++++++++++++--------
src/compute_coord_atom.h | 8 +
src/compute_orientorder_atom.cpp | 32 ++-
src/compute_orientorder_atom.h | 7 +-
6 files changed, 308 insertions(+), 92 deletions(-)
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 <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,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");
-
- 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 = 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++;
+ if (narg < 5) error->all(FLERR,"Illegal compute coord/atom command");
+
+ 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-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++;
+ }
+
}
- }
+
+ } 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,67 +217,133 @@ 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];
-
- n = 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]) n++;
- }
+ if (strcmp(cstyle,"cutoff") == 0) {
- cvec[i] = n;
- } else cvec[i] = 0.0;
- }
+ if (ncol == 1) {
- } 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;
-
- 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;
+ 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;
+
+ 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++;
+ }
+
+ 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;
+
+ 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++];
+ }
+ }
+
+}
+
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
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;
--
GitLab