diff --git a/doc/src/compute_coord_atom.txt b/doc/src/compute_coord_atom.txt
index 012a87a9a7b17bdbae1f38e4c57cf3cbb2d7bf84..f1a6bf7ffb14085cebd25775f050703456e0b766 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 c5ecef3cb29ee3656a79f01965ddcf62a0687baa..74426dd33b2a5d74dc683b8d9bf1e5375ffc5714 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/examples/prd/in.prd b/examples/prd/in.prd
index ea5220ab445516e6128f128fe9cbb007be2a4281..be3454bcb969ab0fb74f7a70a4ef43cb08ad35dc 100644
--- a/examples/prd/in.prd
+++ b/examples/prd/in.prd
@@ -78,7 +78,7 @@ run 100
# only output atoms near vacancy
-compute coord all coord/atom $r
+compute coord all coord/atom cutoff $r
#dump events all custom 1 dump.prd id type x y z
#dump_modify events thresh c_coord != 4
diff --git a/examples/tad/in.tad b/examples/tad/in.tad
index da3d2175aa26def703b7c03f795d503a3f325cc8..687e1dde0cb260ab06c401cc96e76c57e801e8d5 100644
--- a/examples/tad/in.tad
+++ b/examples/tad/in.tad
@@ -80,7 +80,7 @@ velocity all zero linear
# only output atoms near vacancy
-compute coord all coord/atom $r
+compute coord all coord/atom cutoff $r
#dump events all custom 1 dump.prd id type x y z
#dump_modify events thresh c_coord != 4
diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp
index 21744dcc93839b4e1baf264c554b77de8dfe13b2..36f0b635047fb4aadaa8331325694b0d628648da 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,70 @@
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)
+ typelo(NULL), typehi(NULL), cvec(NULL), carray(NULL),
+ id_orientorder(NULL), normv(NULL)
{
- if (narg < 4) error->all(FLERR,"Illegal compute coord/atom command");
+ if (narg < 5) error->all(FLERR,"Illegal compute coord/atom command");
+
+ cstyle = NONE;
+
+ if (strcmp(arg[3],"cutoff") == 0) {
+ cstyle = CUTOFF;
+ 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(arg[3],"orientorder") == 0) {
+ cstyle = ORIENT;
+ if (narg != 6) error->all(FLERR,"Illegal compute coord/atom command");
- double cutoff = force->numeric(FLERR,arg[3]);
- cutsq = cutoff*cutoff;
+ int n = strlen(arg[4]) + 1;
+ id_orientorder = new char[n];
+ strcpy(id_orientorder,arg[4]);
- ncol = narg-4 + 1;
- int ntypes = atom->ntypes;
- typelo = new int[ncol];
- typehi = new int[ncol];
+ 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");
- if (narg == 4) {
ncol = 1;
+ typelo = new int[ncol];
+ typehi = new int[ncol];
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++;
- }
- }
+ typehi[0] = atom->ntypes;
+
+ } else error->all(FLERR,"Invalid cstyle in compute coord/atom");
peratom_flag = 1;
if (ncol == 1) size_peratom_cols = 0;
@@ -76,12 +110,25 @@ ComputeCoordAtom::~ComputeCoordAtom()
delete [] typehi;
memory->destroy(cvec);
memory->destroy(carray);
+ delete [] id_orientorder;
}
/* ---------------------------------------------------------------------- */
void ComputeCoordAtom::init()
{
+ if (cstyle == ORIENT) {
+ 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 (cstyle == ORIENT) {
+ 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,72 @@ void ComputeCoordAtom::compute_peratom()
int *type = atom->type;
int *mask = atom->mask;
- if (ncol == 1) {
+ if (cstyle == CUTOFF) {
+
+ 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++;
+ }
+
+ 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 (cstyle == ORIENT) {
+
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
@@ -168,51 +296,51 @@ void ComputeCoordAtom::compute_peratom()
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 (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;
}
+ }
+}
- } 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];
+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;
+}
- 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;
- }
- }
- }
+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 0ff373f133695b4c90bca9508950479311a3bb3a..2ad46fa854b233e461cd0bc0d53fc7b017a8006e 100644
--- a/src/compute_coord_atom.h
+++ b/src/compute_coord_atom.h
@@ -31,7 +31,10 @@ 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();
+ enum {NONE,CUTOFF,ORIENT};
private:
int nmax,ncol;
@@ -41,6 +44,12 @@ class ComputeCoordAtom : public Compute {
int *typelo,*typehi;
double *cvec;
double **carray;
+
+ class ComputeOrientOrderAtom *c_orientorder;
+ char *id_orientorder;
+ double threshold;
+ double **normv;
+ int cstyle,nqlist,l;
};
}
diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp
index 6c5a2c0c0e10b0df9e0d32521a7e0230ca1c2eb8..5f78b33b61d02e08411dbcfcb91e2492c3625ac1 100644
--- a/src/compute_orientorder_atom.cpp
+++ b/src/compute_orientorder_atom.cpp
@@ -46,7 +46,8 @@ using namespace std;
ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
- distsq(NULL), nearest(NULL), rlist(NULL), qlist(NULL), qnarray(NULL), qnm_r(NULL), qnm_i(NULL)
+ qlist(NULL), distsq(NULL), nearest(NULL), rlist(NULL),
+ qnarray(NULL), qnm_r(NULL), qnm_i(NULL)
{
if (narg < 3 ) error->all(FLERR,"Illegal compute orientorder/atom command");
@@ -54,9 +55,10 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
nnn = 12;
cutsq = 0.0;
+ qlcompflag = 0;
// specify which orders to request
-
+
nqlist = 5;
memory->create(qlist,nqlist,"orientorder/atom:qlist");
qlist[0] = 4;
@@ -72,40 +74,63 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg
while (iarg < narg) {
if (strcmp(arg[iarg],"nnn") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
- if (strcmp(arg[iarg+1],"NULL") == 0)
- nnn = 0;
- else {
- nnn = force->numeric(FLERR,arg[iarg+1]);
- if (nnn <= 0)
- error->all(FLERR,"Illegal compute orientorder/atom command");
+ if (strcmp(arg[iarg+1],"NULL") == 0) {
+ nnn = 0;
+ } else {
+ nnn = force->numeric(FLERR,arg[iarg+1]);
+ if (nnn <= 0)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
}
iarg += 2;
} else if (strcmp(arg[iarg],"degrees") == 0) {
- if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
+ if (iarg+2 > narg)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
nqlist = force->numeric(FLERR,arg[iarg+1]);
- if (nqlist <= 0) error->all(FLERR,"Illegal compute orientorder/atom command");
+ if (nqlist <= 0)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
memory->destroy(qlist);
memory->create(qlist,nqlist,"orientorder/atom:qlist");
iarg += 2;
if (iarg+nqlist > narg) error->all(FLERR,"Illegal compute orientorder/atom command");
qmax = 0;
for (int iw = 0; iw < nqlist; iw++) {
- qlist[iw] = force->numeric(FLERR,arg[iarg+iw]);
- if (qlist[iw] < 0)
- error->all(FLERR,"Illegal compute orientorder/atom command");
- if (qlist[iw] > qmax) qmax = qlist[iw];
+ qlist[iw] = force->numeric(FLERR,arg[iarg+iw]);
+ if (qlist[iw] < 0)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
+ 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");
+ if (iarg+2 > narg)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
double cutoff = force->numeric(FLERR,arg[iarg+1]);
- if (cutoff <= 0.0) error->all(FLERR,"Illegal compute orientorder/atom command");
+ if (cutoff <= 0.0)
+ error->all(FLERR,"Illegal compute orientorder/atom command");
cutsq = cutoff*cutoff;
iarg += 2;
} 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;
@@ -124,7 +149,7 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
memory->destroy(qlist);
memory->destroy(qnm_r);
memory->destroy(qnm_i);
-
+
}
/* ---------------------------------------------------------------------- */
@@ -207,7 +232,7 @@ void ComputeOrientOrderAtom::compute_peratom()
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
-
+
// insure distsq and nearest arrays are long enough
if (jnum > maxneigh) {
@@ -236,9 +261,9 @@ void ComputeOrientOrderAtom::compute_peratom()
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutsq) {
distsq[ncount] = rsq;
- rlist[ncount][0] = delx;
- rlist[ncount][1] = dely;
- rlist[ncount][2] = delz;
+ rlist[ncount][0] = delx;
+ rlist[ncount][1] = dely;
+ rlist[ncount][2] = delz;
nearest[ncount++] = j;
}
}
@@ -246,16 +271,16 @@ void ComputeOrientOrderAtom::compute_peratom()
// if not nnn neighbors, order parameter = 0;
if ((ncount == 0) || (ncount < nnn)) {
- for (int iw = 0; iw < nqlist; iw++)
- qn[iw] = 0.0;
+ for (int iw = 0; iw < nqlist; iw++)
+ qn[iw] = 0.0;
continue;
}
// if nnn > 0, use only nearest nnn neighbors
if (nnn > 0) {
- select3(nnn,ncount,distsq,nearest,rlist);
- ncount = nnn;
+ select3(nnn,ncount,distsq,nearest,rlist);
+ ncount = nnn;
}
calc_boop(rlist, ncount, qn, qlist, nqlist);
@@ -270,8 +295,8 @@ void ComputeOrientOrderAtom::compute_peratom()
double ComputeOrientOrderAtom::memory_usage()
{
double bytes = ncol*nmax * sizeof(double);
- bytes += (qmax*(2*qmax+1)+maxneigh*4) * sizeof(double);
- bytes += (nqlist+maxneigh) * sizeof(int);
+ bytes += (qmax*(2*qmax+1)+maxneigh*4) * sizeof(double);
+ bytes += (nqlist+maxneigh) * sizeof(int);
return bytes;
}
@@ -283,18 +308,18 @@ double ComputeOrientOrderAtom::memory_usage()
// Use no-op do while to create single statement
-#define SWAP(a,b) do { \
- tmp = a; a = b; b = tmp; \
+#define SWAP(a,b) do { \
+ tmp = a; a = b; b = tmp; \
} while(0)
-#define ISWAP(a,b) do { \
- itmp = a; a = b; b = itmp; \
+#define ISWAP(a,b) do { \
+ itmp = a; a = b; b = itmp; \
} while(0)
-#define SWAP3(a,b) do { \
- tmp = a[0]; a[0] = b[0]; b[0] = tmp; \
- tmp = a[1]; a[1] = b[1]; b[1] = tmp; \
- tmp = a[2]; a[2] = b[2]; b[2] = tmp; \
+#define SWAP3(a,b) do { \
+ tmp = a[0]; a[0] = b[0]; b[0] = tmp; \
+ tmp = a[1]; a[1] = b[1]; b[1] = tmp; \
+ tmp = a[2]; a[2] = b[2]; b[2] = tmp; \
} while(0)
/* ---------------------------------------------------------------------- */
@@ -313,7 +338,7 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
if (ir <= l+1) {
if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir]);
- ISWAP(iarr[l],iarr[ir]);
+ ISWAP(iarr[l],iarr[ir]);
SWAP3(arr3[l],arr3[ir]);
}
return;
@@ -325,17 +350,17 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir]);
ISWAP(iarr[l],iarr[ir]);
- SWAP3(arr3[l],arr3[ir]);
+ SWAP3(arr3[l],arr3[ir]);
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir]);
ISWAP(iarr[l+1],iarr[ir]);
- SWAP3(arr3[l+1],arr3[ir]);
+ SWAP3(arr3[l+1],arr3[ir]);
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1]);
ISWAP(iarr[l],iarr[l+1]);
- SWAP3(arr3[l],arr3[l+1]);
+ SWAP3(arr3[l],arr3[l+1]);
}
i = l+1;
j = ir;
@@ -350,7 +375,7 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
if (j < i) break;
SWAP(arr[i],arr[j]);
ISWAP(iarr[i],iarr[j]);
- SWAP3(arr3[i],arr3[j]);
+ SWAP3(arr3[i],arr3[j]);
}
arr[l+1] = arr[j];
arr[j] = a;
@@ -372,9 +397,9 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl
calculate the bond orientational order parameters
------------------------------------------------------------------------- */
-void ComputeOrientOrderAtom::calc_boop(double **rlist,
- int ncount, double qn[],
- int qlist[], int nqlist) {
+void ComputeOrientOrderAtom::calc_boop(double **rlist,
+ int ncount, double qn[],
+ int qlist[], int nqlist) {
for (int iw = 0; iw < nqlist; iw++) {
int n = qlist[iw];
@@ -412,37 +437,50 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist,
double expphim_r = expphi_r;
double expphim_i = expphi_i;
for(int m = 1; m <= +n; m++) {
- double prefactor = polar_prefactor(n, m, costheta);
- double c_r = prefactor * expphim_r;
- double c_i = prefactor * expphim_i;
- qnm_r[iw][m+n] += c_r;
- qnm_i[iw][m+n] += c_i;
- if(m & 1) {
- qnm_r[iw][-m+n] -= c_r;
- qnm_i[iw][-m+n] += c_i;
- } else {
- qnm_r[iw][-m+n] += c_r;
- qnm_i[iw][-m+n] -= c_i;
- }
- double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i;
- double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r;
- expphim_r = tmp_r;
- expphim_i = tmp_i;
+ double prefactor = polar_prefactor(n, m, costheta);
+ double c_r = prefactor * expphim_r;
+ double c_i = prefactor * expphim_i;
+ qnm_r[iw][m+n] += c_r;
+ qnm_i[iw][m+n] += c_i;
+ if(m & 1) {
+ qnm_r[iw][-m+n] -= c_r;
+ qnm_i[iw][-m+n] += c_i;
+ } else {
+ qnm_r[iw][-m+n] += c_r;
+ qnm_i[iw][-m+n] -= c_i;
+ }
+ double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i;
+ double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r;
+ expphim_r = tmp_r;
+ expphim_i = tmp_i;
}
}
}
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;
for(int m = 0; m < 2*n+1; m++) {
qm_sum += qnm_r[iw][m]*qnm_r[iw][m] + qnm_i[iw][m]*qnm_i[iw][m];
// printf("Ylm^2 = %d %d %g\n",n,m,
- // qnm_r[iw][m]*qnm_r[iw][m] + qnm_i[iw][m]*qnm_i[iw][m]);
+ // 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;
+ }
}
}
@@ -455,7 +493,7 @@ double ComputeOrientOrderAtom::dist(const double r[]) {
}
/* ----------------------------------------------------------------------
- polar prefactor for spherical harmonic Y_l^m, where
+ polar prefactor for spherical harmonic Y_l^m, where
Y_l^m (theta, phi) = prefactor(l, m, cos(theta)) * exp(i*m*phi)
------------------------------------------------------------------------- */
diff --git a/src/compute_orientorder_atom.h b/src/compute_orientorder_atom.h
index 9c5ec14f562efe1300c19719ce22456e6d001cc8..81b08dbddc17740538545d8ea7b6add7cd90ed7a 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,16 +43,13 @@ 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;
void select3(int, int, double *, int *, double **);
- void calc_boop(double **rlist, int numNeighbors,
+ void calc_boop(double **rlist, int numNeighbors,
double qn[], int nlist[], int nnlist);
double dist(const double r[]);