Merge pull request #332 from akohlmey/coord-atom-orientorder-atom-enhancements (973bef4d) · Commits · 郑智淋 / lammps

doc/src/compute_coord_atom.txt

+37 −8

Original line number	Diff line number	Diff line
		@@ -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 6.0 1 2
		compute 1 all coord/atom 6.0 24 58 * :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 24 58 *
		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

+22 −3

Original line number	Diff line number	Diff line
		@@ -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).

examples/prd/in.prd

+1 −1

Original line number	Diff line number	Diff line
		@@ -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

examples/tad/in.tad

+1 −1

Original line number	Diff line number	Diff line
		@@ -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

src/compute_coord_atom.cpp

+180 −52

Original line number	Diff line number	Diff line
		@@ -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,36 @@

		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;

		double cutoff = force->numeric(FLERR,arg[3]);
		if (strcmp(arg[3],"cutoff") == 0) {
		cstyle = CUTOFF;
		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])
		@@ -61,6 +69,32 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS lmp, int narg, char *arg) :
		}
		}

		} else if (strcmp(arg[3],"orientorder") == 0) {
		cstyle = ORIENT;
		if (narg != 6) error->all(FLERR,"Illegal compute coord/atom command");

		int 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;
		@@ -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(2l+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,10 @@ void ComputeCoordAtom::compute_peratom()
		int *type = atom->type;
		int *mask = atom->mask;

		if (cstyle == CUTOFF) {

		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 = delxdelx + delydely + 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,67 @@ void ComputeCoordAtom::compute_peratom()
		}
		}
		}

		} else if (cstyle == ORIENT) {

		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 = delxdelx + delydely + delz*delz;
		if (rsq < cutsq) {
		double dot_product = 0.0;
		for (int m=0; m < 2(2l+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(2l+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(2l+1); ++j) {
		normv[i][j] = buf[m++];
		}
		}

		}

		/* ----------------------------------------------------------------------

Admin message