diff --git a/doc/src/Eqs/cnp_cutoff.jpg b/doc/src/Eqs/cnp_cutoff.jpg new file mode 100644 index 0000000000..fae5c6b636 Binary files /dev/null and b/doc/src/Eqs/cnp_cutoff.jpg differ diff --git a/doc/src/Eqs/cnp_cutoff.tex b/doc/src/Eqs/cnp_cutoff.tex new file mode 100644 index 0000000000..f74166e608 --- /dev/null +++ b/doc/src/Eqs/cnp_cutoff.tex @@ -0,0 +1,14 @@ +\documentclass[12pt,article]{article} + +\usepackage{indentfirst} +\usepackage{amsmath} + +\begin{document} + +\begin{eqnarray*} + r_{c}^{fcc} & = & \frac{1}{2} \left(\frac{\sqrt{2}}{2} + 1\right) \mathrm{a} \simeq 0.8536 \:\mathrm{a} \\ + r_{c}^{bcc} & = & \frac{1}{2}(\sqrt{2} + 1) \mathrm{a} \simeq 1.207 \:\mathrm{a} \\ + r_{c}^{hcp} & = & \frac{1}{2}\left(1+\sqrt{\frac{4+2x^{2}}{3}}\right) \mathrm{a} +\end{eqnarray*} + +\end{document} diff --git a/doc/src/Eqs/cnp_cutoff2.jpg b/doc/src/Eqs/cnp_cutoff2.jpg new file mode 100644 index 0000000000..744b61e9b4 Binary files /dev/null and b/doc/src/Eqs/cnp_cutoff2.jpg differ diff --git a/doc/src/Eqs/cnp_cutoff2.tex b/doc/src/Eqs/cnp_cutoff2.tex new file mode 100644 index 0000000000..fcec31fd24 --- /dev/null +++ b/doc/src/Eqs/cnp_cutoff2.tex @@ -0,0 +1,12 @@ +\documentclass[12pt,article]{article} + +\usepackage{indentfirst} +\usepackage{amsmath} + +\begin{document} + +$$ + Rc + Rs > 2*{\rm cutoff} +$$ + +\end{document} diff --git a/doc/src/Eqs/cnp_eq.jpg b/doc/src/Eqs/cnp_eq.jpg new file mode 100644 index 0000000000..d421314442 Binary files /dev/null and b/doc/src/Eqs/cnp_eq.jpg differ diff --git a/doc/src/Eqs/cnp_eq.tex b/doc/src/Eqs/cnp_eq.tex new file mode 100644 index 0000000000..e5f157e6ba --- /dev/null +++ b/doc/src/Eqs/cnp_eq.tex @@ -0,0 +1,9 @@ +\documentclass[12pt]{article} + +\begin{document} + +$$ + Q_{i} = \frac{1}{n_i}\sum_{j = 1}^{n_i} | \sum_{k = 1}^{n_{ij}} \vec{R}_{ik} + \vec{R}_{jk} |^2 +$$ + +\end{document} diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt index dc7ddebe58..4895a0729c 100644 --- a/doc/src/Section_commands.txt +++ b/doc/src/Section_commands.txt @@ -831,6 +831,7 @@ package"_Section_start.html#start_3. "ackland/atom"_compute_ackland_atom.html, "basal/atom"_compute_basal_atom.html, +"cnp/atom"_compute_cnp_atom.html, "dpd"_compute_dpd.html, "dpd/atom"_compute_dpd_atom.html, "fep"_compute_fep.html, diff --git a/doc/src/compute_cna_atom.txt b/doc/src/compute_cna_atom.txt index 74240b515d..23289b0132 100644 --- a/doc/src/compute_cna_atom.txt +++ b/doc/src/compute_cna_atom.txt @@ -26,7 +26,7 @@ Define a computation that calculates the CNA (Common Neighbor Analysis) pattern for each atom in the group. In solid-state systems the CNA pattern is a useful measure of the local crystal structure around an atom. The CNA methodology is described in "(Faken)"_#Faken -and "(Tsuzuki)"_#Tsuzuki. +and "(Tsuzuki)"_#Tsuzuki1. Currently, there are five kinds of CNA patterns LAMMPS recognizes: @@ -93,5 +93,5 @@ above. :link(Faken) [(Faken)] Faken, Jonsson, Comput Mater Sci, 2, 279 (1994). -:link(Tsuzuki) +:link(Tsuzuki1) [(Tsuzuki)] Tsuzuki, Branicio, Rino, Comput Phys Comm, 177, 518 (2007). diff --git a/doc/src/compute_cnp_atom.txt b/doc/src/compute_cnp_atom.txt new file mode 100644 index 0000000000..2e9950b4c1 --- /dev/null +++ b/doc/src/compute_cnp_atom.txt @@ -0,0 +1,111 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +compute cnp/atom command :h3 + +[Syntax:] + +compute ID group-ID cnp/atom cutoff :pre + +ID, group-ID are documented in "compute"_compute.html command +cnp/atom = style name of this compute command +cutoff = cutoff distance for nearest neighbors (distance units) :ul + +[Examples:] + +compute 1 all cnp/atom 3.08 :pre + +[Description:] + +Define a computation that calculates the Common Neighborhood +Parameter (CNP) for each atom in the group. In solid-state systems +the CNP is a useful measure of the local crystal structure +around an atom and can be used to characterize whether the +atom is part of a perfect lattice, a local defect (e.g. a dislocation +or stacking fault), or at a surface. + +The value of the CNP parameter will be 0.0 for atoms not +in the specified compute group. Note that normally a CNP calculation should only be +performed on mono-component systems. + +This parameter is computed using the following formula from +"(Tsuzuki)"_#Tsuzuki2 + +:c,image(Eqs/cnp_eq.jpg) + +where the index {j} goes over the {n}i nearest neighbors of atom +{i}, and the index {k} goes over the {n}ij common nearest neighbors +between atom {i} and atom {j}. Rik and Rjk are the vectors connecting atom +{k} to atoms {i} and {j}. The quantity in the double sum is computed +for each atom. + +The CNP calculation is sensitive to the specified cutoff value. +You should ensure that the appropriate nearest neighbors of an atom are +found within the cutoff distance for the presumed crystal structure. +E.g. 12 nearest neighbor for perfect FCC and HCP crystals, 14 nearest +neighbors for perfect BCC crystals. These formulas can be used to +obtain a good cutoff distance: + +:c,image(Eqs/cnp_cutoff.jpg) + +where a is the lattice constant for the crystal structure concerned +and in the HCP case, x = (c/a) / 1.633, where 1.633 is the ideal c/a +for HCP crystals. + +Also note that since the CNP calculation in LAMMPS uses the neighbors +of an owned atom to find the nearest neighbors of a ghost atom, the +following relation should also be satisfied: + +:c,image(Eqs/cnp_cutoff2.jpg) + +where Rc is the cutoff distance of the potential, Rs is the skin +distance as specified by the "neighbor"_neighbor.html command, and +cutoff is the argument used with the compute cnp/atom command. LAMMPS +will issue a warning if this is not the case. + +The neighbor list needed to compute this quantity is constructed each +time the calculation is performed (e.g. each time a snapshot of atoms +is dumped). Thus it can be inefficient to compute/dump this quantity +too frequently or to have multiple compute/dump commands, each with a +{cnp/atom} style. + +[Output info:] + +This compute calculates a per-atom vector, which 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 options. + +The per-atom vector values will be real positive numbers. Some typical CNP +values: + +FCC lattice = 0.0 +BCC lattice = 0.0 +HCP lattice = 4.4 :pre + +FCC (111) surface ~ 13.0 +FCC (100) surface ~ 26.5 +FCC dislocation core ~ 11 :pre + +[Restrictions:] + +This compute is part of the USER-MISC package. It is only enabled if +LAMMPS was built with that package. See the "Making +LAMMPS"_Section_start.html#start_3 section for more info. + +[Related commands:] + +"compute cna/atom"_compute_cna_atom.html +"compute centro/atom"_compute_centro_atom.html + +[Default:] none + +:line + +:link(Tsuzuki2) +[(Tsuzuki)] Tsuzuki, Branicio, Rino, Comput Phys Comm, 177, 518 (2007). diff --git a/doc/src/computes.txt b/doc/src/computes.txt index 1d01798791..5a6ca66c46 100644 --- a/doc/src/computes.txt +++ b/doc/src/computes.txt @@ -17,6 +17,7 @@ Computes :h1 compute_chunk_atom compute_cluster_atom compute_cna_atom + compute_cnp_atom compute_com compute_com_chunk compute_contact_atom diff --git a/doc/src/lammps.book b/doc/src/lammps.book index 1769f29825..79ff5b4a75 100644 --- a/doc/src/lammps.book +++ b/doc/src/lammps.book @@ -301,6 +301,7 @@ compute_centro_atom.html compute_chunk_atom.html compute_cluster_atom.html compute_cna_atom.html +compute_cnp_atom.html compute_com.html compute_com_chunk.html compute_contact_atom.html diff --git a/src/.gitignore b/src/.gitignore index 0cddfa6951..d6c535ebe9 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -177,8 +177,8 @@ /compute_basal_atom.h /compute_body_local.cpp /compute_body_local.h -/compute_cna_atom2.cpp -/compute_cna_atom2.h +/compute_cnp_atom.cpp +/compute_cnp_atom.h /compute_damage_atom.cpp /compute_damage_atom.h /compute_dilatation_atom.cpp diff --git a/src/USER-MISC/README b/src/USER-MISC/README index cacee41e0c..f501c81c17 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -28,6 +28,7 @@ bond_style harmonic/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 bond_style harmonic/shift/cut, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007 compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013 +compute cnp/atom, Paulo Branicio (USC), branicio at usc.edu, 31 May 2017 compute temp/rotate, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 compute PRESSURE/GREM, David Stelter, dstelter@bu.edu, 22 Nov 16 dihedral_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 diff --git a/src/USER-MISC/compute_cnp_atom.cpp b/src/USER-MISC/compute_cnp_atom.cpp new file mode 100644 index 0000000000..6dce41291e --- /dev/null +++ b/src/USER-MISC/compute_cnp_atom.cpp @@ -0,0 +1,311 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + + Common Neighbor Parameter as proposed in: + Tsuzuki, Branicio, Rino, Comput Phys Comm, 177, 518 (2007) + Cite: http://dx.doi.org/10.1063/1.2197987 + +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Paulo Branicio (University of Southern California) + branicio@usc.edu +------------------------------------------------------------------------- */ + +#include +#include +#include + +#include "compute_cnp_atom.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "force.h" +#include "pair.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +//define maximum values +#define MAXNEAR 24 +#define MAXCOMMON 12 + +enum{NCOMMON}; + +/* ---------------------------------------------------------------------- */ + +ComputeCNPAtom::ComputeCNPAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), + nearest(NULL), nnearest(NULL), cnpv(NULL) +{ + if (narg != 4) error->all(FLERR,"Illegal compute cnp/atom command"); + + peratom_flag = 1; + size_peratom_cols = 0; + + double cutoff = force->numeric(FLERR,arg[3]); + if (cutoff < 0.0) error->all(FLERR,"Illegal compute cnp/atom command"); + cutsq = cutoff*cutoff; + + nmax = 0; +} + +/* ---------------------------------------------------------------------- */ + +ComputeCNPAtom::~ComputeCNPAtom() +{ + memory->destroy(nearest); + memory->destroy(nnearest); + memory->destroy(cnpv); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeCNPAtom::init() +{ + if (force->pair == NULL) + error->all(FLERR,"Compute cnp/atom requires a pair style be defined"); + + if (sqrt(cutsq) > force->pair->cutforce) + error->all(FLERR,"Compute cnp/atom cutoff is longer than pairwise cutoff"); + + if (2.0*sqrt(cutsq) > force->pair->cutforce + neighbor->skin && + comm->me == 0) + error->warning(FLERR,"Compute cnp/atom cutoff may be too large to find " + "ghost atom neighbors"); + + int count = 0; + for (int i = 0; i < modify->ncompute; i++) + if (strcmp(modify->compute[i]->style,"cnp/atom") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning(FLERR,"More than one compute cnp/atom defined"); + + // need an occasional full neighbor list + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + neighbor->requests[irequest]->occasional = 1; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeCNPAtom::init_list(int id, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeCNPAtom::compute_peratom() +{ + int i,j,k,ii,jj,kk,m,n,inum,jnum,inear,jnear; + int firstflag,ncommon; + int *ilist,*jlist,*numneigh,**firstneigh; + int cnp[MAXNEAR][4],onenearest[MAXNEAR]; + int common[MAXCOMMON]; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double xjtmp,yjtmp,zjtmp,rjkx,rjky,rjkz; + + invoked_peratom = update->ntimestep; + + // grow arrays if necessary + + if (atom->nmax > nmax) { + memory->destroy(nearest); + memory->destroy(nnearest); + memory->destroy(cnpv); + nmax = atom->nmax; + memory->create(nearest,nmax,MAXNEAR,"cnp:nearest"); + memory->create(nnearest,nmax,"cnp:nnearest"); + memory->create(cnpv,nmax,"cnp:cnp_cnpv"); + vector_atom = cnpv; + } + + // invoke full neighbor list (will copy or build if necessary) + + neighbor->build_one(list); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // find the neigbors of each atom within cutoff using full neighbor list + // nearest[] = atom indices of nearest neighbors, up to MAXNEAR + // do this for all atoms, not just compute group + // since CNP calculation requires neighbors of neighbors + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int nerror = 0; + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + 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) { + if (n < MAXNEAR) nearest[i][n++] = j; + else { + nerror++; + break; + } + } + } + nnearest[i] = n; + } + + // warning message + + int nerrorall; + MPI_Allreduce(&nerror,&nerrorall,1,MPI_INT,MPI_SUM,world); + if (nerrorall && comm->me == 0) { + char str[128]; + sprintf(str,"Too many neighbors in CNP for %d atoms",nerrorall); + error->warning(FLERR,str,0); + } + + // compute CNP value for each atom in group + // only performed if # of nearest neighbors = 12 or 14 (fcc,hcp) + + nerror = 0; + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + // reset cnpv + cnpv[i] = 0.0; + + // loop over nearest neighbors of I to build cnp data structure + // cnp[k][NCOMMON] = # of common neighbors of I with each of its neighbors + for (m = 0; m < nnearest[i]; m++) { + j = nearest[i][m]; + xjtmp = x[j][0]; + yjtmp = x[j][1]; + zjtmp = x[j][2]; + + // common = list of neighbors common to atom I and atom J + // if J is an owned atom, use its near neighbor list to find them + // if J is a ghost atom, use full neighbor list of I to find them + // in latter case, must exclude J from I's neighbor list + + // find common neighbors of i and j using near neighbor list + if (j < nlocal) { + firstflag = 1; + ncommon = 0; + for (inear = 0; inear < nnearest[i]; inear++) + for (jnear = 0; jnear < nnearest[j]; jnear++) + if (nearest[i][inear] == nearest[j][jnear]) { + if (ncommon < MAXCOMMON) common[ncommon++] = nearest[i][inear]; + else if (firstflag) { + nerror++; + firstflag = 0; + } + } + + // find common neighbors of i and j using full neighbor list + } else { + jlist = firstneigh[i]; + jnum = numneigh[i]; + + n = 0; + for (kk = 0; kk < jnum; kk++) { + k = jlist[kk]; + k &= NEIGHMASK; + if (k == j) continue; + + delx = xjtmp - x[k][0]; + dely = yjtmp - x[k][1]; + delz = zjtmp - x[k][2]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq < cutsq) { + if (n < MAXNEAR) onenearest[n++] = k; + else break; + } + } + + firstflag = 1; + ncommon = 0; + for (inear = 0; inear < nnearest[i]; inear++) + for (jnear = 0; (jnear < n) && (n < MAXNEAR); jnear++) + if (nearest[i][inear] == onenearest[jnear]) { + if (ncommon < MAXCOMMON) common[ncommon++] = nearest[i][inear]; + else if (firstflag) { + nerror++; + firstflag = 0; + } + } + } + + // Calculate and update sum |Rik+Rjk|ˆ2 + rjkx = 0.0; + rjky = 0.0; + rjkz = 0.0; + for (kk = 0; kk < ncommon; kk++) { + k = common[kk]; + rjkx += 2.0*x[k][0] - xjtmp - xtmp; + rjky += 2.0*x[k][1] - yjtmp - ytmp; + rjkz += 2.0*x[k][2] - zjtmp - ztmp; + } + // update cnpv with summed (valuejk) + cnpv[i] += rjkx*rjkx + rjky*rjky + rjkz*rjkz; + + // end of loop over j atoms + } + + // normalize cnp by the number of nearest neighbors + cnpv[i] = cnpv[i] / nnearest[i]; + + // end of loop over i atoms + } + + // warning message + MPI_Allreduce(&nerror,&nerrorall,1,MPI_INT,MPI_SUM,world); + if (nerrorall && comm->me == 0) { + char str[128]; + sprintf(str,"Too many common neighbors in CNP %d times",nerrorall); + error->warning(FLERR,str); + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeCNPAtom::memory_usage() +{ + double bytes = nmax * sizeof(int); + bytes += nmax * MAXNEAR * sizeof(int); + bytes += nmax * sizeof(double); + return bytes; +} diff --git a/src/USER-MISC/compute_cnp_atom.h b/src/USER-MISC/compute_cnp_atom.h new file mode 100644 index 0000000000..4fdb3954f2 --- /dev/null +++ b/src/USER-MISC/compute_cnp_atom.h @@ -0,0 +1,92 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(cnp/atom,ComputeCNPAtom) + +#else + +#ifndef LMP_COMPUTE_CNP_ATOM_H +#define LMP_COMPUTE_CNP_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeCNPAtom : public Compute { + public: + ComputeCNPAtom(class LAMMPS *, int, char **); + ~ComputeCNPAtom(); + void init(); + void init_list(int, class NeighList *); + void compute_peratom(); + double memory_usage(); + + private: +//revise + int nmax; + double cutsq; + class NeighList *list; + int **nearest; + int *nnearest; + double *cnpv; +// int nmax; +// double cutsq; +// class NeighList *list; +// int **nearest; +// int *nnearest; +// double *pattern; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Compute cnp/atom requires a pair style be defined + +Self-explanatory. + +E: Compute cnp/atom cutoff is longer than pairwise cutoff + +Self-explanatory. + +W: Compute cnp/atom cutoff may be too large to find ghost atom neighbors + +The neighbor cutoff used may not encompass enough ghost atoms +to perform this operation correctly. + +W: More than one compute cnp/atom defined + +It is not efficient to use compute cnp/atom more than once. + +W: Too many neighbors in CNP for %d atoms + +More than the maximum # of neighbors was found multiple times. This +was unexpected. + +W: Too many common neighbors in CNP %d times + +More than the maximum # of neighbors was found multiple times. This +was unexpected. + +*/