integrate compute cnp/atom contributed by Paulo Branicio (USC)

This commit is contained in:
Axel Kohlmeyer 2017-05-31 00:38:44 -04:00
parent 286d4f2743
commit 85e917ae52
15 changed files with 557 additions and 4 deletions

BIN
doc/src/Eqs/cnp_cutoff.jpg Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 13 KiB

View File

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

BIN
doc/src/Eqs/cnp_cutoff2.jpg Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 2.5 KiB

View File

@ -0,0 +1,12 @@
\documentclass[12pt,article]{article}
\usepackage{indentfirst}
\usepackage{amsmath}
\begin{document}
$$
Rc + Rs > 2*{\rm cutoff}
$$
\end{document}

BIN
doc/src/Eqs/cnp_eq.jpg Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 23 KiB

9
doc/src/Eqs/cnp_eq.tex Normal file
View File

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

View File

@ -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,

View File

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

View File

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

View File

@ -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

View File

@ -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

4
src/.gitignore vendored
View File

@ -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

View File

@ -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

View File

@ -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 <string.h>
#include <stdlib.h>
#include <math.h>
#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;
}

View File

@ -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.
*/