jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Creatd new compute ADF for angular distribution function

This commit is contained in:
Aidan Thompson 2018-10-19 18:36:11 -06:00
parent 52f02f2bbb
commit 04a4a29fcf
3 changed files with 876 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

207
doc/src/compute_adf.txt Normal file
View File

@ -0,0 +1,207 @@
"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 adf command :h3
[Syntax:]
compute ID group-ID adf Nbin itype1 jtype1 ktype1 Rjinner1 Rjouter1 Rkinner1 Rkouter1 ... :pre
ID, group-ID are documented in "compute"_compute.html command :ulb,l
adf = style name of this compute command :l
Nbin = number of ADF bins :l
itypeN = central atom type for Nth ADF histogram (see asterisk form below) :l
jtypeN = J atom type for Nth ADF histogram (see asterisk form below) :l
ktypeN = K atom type for Nth ADF histogram (see asterisk form below) :l
RjinnerN = inner radius of J atom shell for Nth ADF histogram (distance units) :l
RjouterN = outer radius of J atom shell for Nth ADF histogram (distance units) :l
RkinnerN = inner radius of K atom shell for Nth ADF histogram (distance units) :l
RkouterN = outer radius of K atom shell for Nth ADF histogram (distance units) :l
zero or one keyword/value pairs may be appended :l
keyword = {ordinate} :l
{ordinate} value = {degree} or {radian} or {cosine}
Choose the ordinate parameter for the histogram :pre
:ule
[Examples:]
compute 1 fluid adf 32 1 1 1 0.0 1.2 0.0 1.2 &
1 1 2 0.0 1.2 0.0 1.5 &
1 2 2 0.0 1.5 0.0 1.5 &
2 1 1 0.0 1.2 0.0 1.2 &
2 1 2 0.0 1.5 2.0 3.5 &
2 2 2 2.0 3.5 2.0 3.5
compute 1 fluid adf 32 1*2 1*2 1*2 0.5 3.5
compute 1 fluid adf 32 :pre
[Description:]
Define a computation that calculates one or more angular distribution functions
(ADF) for a group of particles. Each ADF is calculated in histogram form
by measuring the angle formed by a central atom and two neighbor atoms and
binning these angles into {Nbin} bins.
Only neighbors for which {Rinner} < {R} < {Router} are counted, where
{Rinner} and {Router} are specified separately for the first and second
neighbor atom in each requested ADF.
NOTE: If you have a bonded system, then the settings of
"special_bonds"_special_bonds.html command can remove pairwise
interactions between atoms in the same bond, angle, or dihedral. This
is the default setting for the "special_bonds"_special_bonds.html
command, and means those pairwise interactions do not appear in the
neighbor list. Because this fix uses a neighbor list, it also means
those pairs will not be included in the ADF. This does not apply when
using long-range coulomb interactions ({coul/long}, {coul/msm},
{coul/wolf} or similar. One way to get around this would be to set
special_bond scaling factors to very tiny numbers that are not exactly
zero (e.g. 1.0e-50). Another workaround is to write a dump file, and
use the "rerun"_rerun.html command to compute the ADF for snapshots in
the dump file. The rerun script can use a
"special_bonds"_special_bonds.html command that includes all pairs in
the neighbor list.
NOTE: If you request any outer cutoff {Router} > force cutoff, or if no
pair style is defined, e.g. the "rerun"_rerun.html command is being used to
post-process a dump file of snapshots you must insure ghost atom information
out to the largest value of {Router} + {skin} is communicated, via the
"comm_modify cutoff"_comm_modify.html command, else the ADF computation
cannot be performed, and LAMMPS will give an error message. The {skin} value
is what is specified with the "neighbor"_neighbor.html command.
The {itypeN},{jtypeN},{ktypeN} settings can be specified in one of two
ways. An explicit numeric value can be used, as in the 1st example
above. Or a wild-card asterisk can be used to specify a range of atom
types as in the 2nd example above.
This takes the form "*" or "*n" or "n*" or "m*n". If N = the
number of atom types, then an asterisk with no numeric values means
all types from 1 to N. A leading asterisk means all types from 1 to n
(inclusive). A trailing asterisk means all types from n to N
(inclusive). A middle asterisk means all types from m to n
(inclusive).
If {itypeN}, {jtypeN}, and {ktypeN} are single values, as in the 1st example
above, this means that the ADF is computed where atoms of type {itypeN}
are the central atom, and neighbor atoms of type {jtypeN} and {ktypeN}
are forming the angle. If any of {itypeN}, {jtypeN}, or {ktypeN}
represent a range of values via
the wild-card asterisk, as in the 2nd example above, this means that the
ADF is computed where atoms of any of the range of types represented
by {itypeN} are the central atom, and the angle is formed by two neighbors,
one neighbor in the range of types represented by {jtypeN} and another neighbor
in the range of types represented by {ktypeN}.
If no {itypeN}, {jtypeN}, {ktypeN} settings are specified, then
LAMMPS will generate a single ADF for all atoms in the group.
The inner cutoff is set to zero and the outer cutoff is set
to the force cutoff. If no pair_style is specified, there is no
force cutoff and LAMMPS will give an error message. Note that
in most cases, generating an ADF for all atoms is not a good thing.
Such an ADF is both uninformative and
extremely expensive to compute. For example, with liquid water
with a 10 A force cutoff, there are 80,000 angles per atom.
In addition, most of the interesting angular structure occurs for
neighbors that are the closest to the central atom, involving
just a few dozen angles.
Angles for each ADF are generated by double-looping over the list of
neighbors of each central atom I,
just as they would be in the force calculation for
a threebody potential such as "Stillinger-Weber"_pair_sw.html.
The angle formed by central atom I and neighbor atoms J and K is included in an
ADF if the following criteria are met:
atoms I,J,K are all in the specified compute group
the distance between atoms I,J is between Rjinner and Rjouter
the distance between atoms I,K is between Rkinner and Rkouter
the type of the I atom matches itypeN (one or a range of types)
atoms I,J,K are distinct
the type of the J atom matches jtypeN (one or a range of types)
the type of the K atom matches ktypeN (one or a range of types) :ul
Each unique angle satisfying the above criteria is counted only once, regardless
of whether either or both of the neighbor atoms making up the
angle appear in both the J and K lists.
It is OK if a particular angle is included in more than
one individual histogram, due to the way the {itypeN}, {jtypeN}, {ktypeN}
arguments are specified.
The first ADF value for a bin is calculated from the histogram count by
dividing by the total number of triples satisfying the criteria,
so that the integral of the ADF w.r.t. angle is one i.e. the ADF
is a probability density function.
The second ADF value is reported as a cumulative sum of
all bins up to the current bins, averaged
over atoms of type {itypeN}. It represents the
number of angles per central atom with angle less
than or equal to the angle of the current bin,
analogous to the coordination
number radial distribution function.
The {ordinate} optional keyword determines
whether the bins are of uniform angular size from zero
to 180 ({degree}), zero to Pi ({radian}), or the
cosine of the angle uniform in the range \[-1,1\] ({cosine}).
Regardless of which is chosen, the first ADF will be normalized
so that the integral w.r.t. the ordinate is one.
The simplest way to output the results of the compute adf calculation
to a file is to use the "fix ave/time"_fix_ave_time.html command, for
example:
compute myADF all adf 32 2 2 2 0.5 3.5 0.5 3.5
fix 1 all ave/time 100 1 100 c_myADF\[*\] file tmp.adf mode vector :pre
[Output info:]
This compute calculates a global array with the number of rows =
{Nbins}, and the number of columns = 1 + 2*Ntriples, where Ntriples is the
number of I,J,K triples specified. The first column has the bin
coordinate (angle-related ordinate at midpoint of bin). Each subsequent column has
the two ADF values for a specific set of ({itypeN},{jtypeN},{ktypeN})
interactions, as described above. These values can be used
by any command that uses a global values from a compute as input. See
"Section 6.15"_Section_howto.html#howto_15 for an overview of
LAMMPS output options.
The array values calculated by this compute are all "intensive".
The first column of array values is the angle-related ordinate, either
the angle in degrees or radians, or the cosine of the angle. Each
subsequent pair of columns gives the first and second kinds of ADF
for a specific set of ({itypeN},{jtypeN},{ktypeN}). The values
in the first ADF column are normalized numbers >= 0.0,
whose integral w.r.t. the ordinate is one,
i.e. the first ADF is a normalized probability distribution.
The values in the second ADF column are also numbers >= 0.0.
They are the cumulative density distribution of angles per atom.
By definition, this ADF is monotonically increasing from zero to
a maximum value equal to the average total number of
angles per atom satisfying the ADF criteria.
[Restrictions:]
The ADF is not computed for neighbors outside the force cutoff,
since processors (in parallel) don't know about atom coordinates for
atoms further away than that distance. If you want an ADF for larger
distances, you can use the "rerun"_rerun.html command to post-process
a dump file and set the cutoff for the potential to be longer in the
rerun script. Note that in the rerun context, the force cutoff is
arbitrary, since you aren't running dynamics and thus are not changing
your model.
[Related commands:]
"compute rdf"_compute_rdf.html, "fix ave/time"_fix_ave_time.html, "compute_modify"_compute_modify.html
[Default:]
The keyword default is ordinate = degree.

580
src/compute_adf.cpp Normal file
View File

@ -0,0 +1,580 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Aidan P. Thompson (SNL)
------------------------------------------------------------------------- */
#include <mpi.h>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "compute_adf.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "neigh_list.h"
#include "group.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "comm.h"
using namespace LAMMPS_NS;
using namespace MathConst;
enum{DEGREE, RADIAN, COSINE};
/* ----------------------------------------------------------------------
compute angular distribution functions for I, J, K atoms
---------------------------------------------------------------------- */
ComputeADF::ComputeADF(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
ilo(NULL), ihi(NULL), jlo(NULL), jhi(NULL), klo(NULL), khi(NULL),
iatomcount(NULL), iatomcountall(NULL),
hist(NULL), histall(NULL),
iatomflag(NULL),
jatomflag(NULL), rcutinnerj(NULL), rcutouterj(NULL),
katomflag(NULL), rcutinnerk(NULL), rcutouterk(NULL),
maxjatom(NULL), numjatom(NULL), neighjatom(NULL),
maxkatom(NULL), numkatom(NULL), neighkatom(NULL),
maxjkatom(NULL), numjkatom(NULL), neighjkatom(NULL), bothjkatom(NULL)
{
int nargsperadf = 7;
if (narg < 4 ) error->all(FLERR,"Illegal compute adf command");
array_flag = 1;
extarray = 0;
// default values
ordinate_style = DEGREE;
cutflag = 0;
nbin = force->inumeric(FLERR,arg[3]);
if (nbin < 1) error->all(FLERR,"Illegal compute adf command");
// optional args
// nargtriple = # of triplewise args, starting at iarg = 4
int iarg;
for (iarg = 4; iarg < narg; iarg++)
if (strcmp(arg[iarg],"ordinate") == 0) break;
int nargtriple = iarg - 4;
// loop over keywords
while (iarg < narg) {
if (strcmp(arg[iarg],"ordinate") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute adf command");
if (strcmp(arg[iarg+1],"degree") == 0) ordinate_style = DEGREE;
else if (strcmp(arg[iarg+1],"radian") == 0) ordinate_style = RADIAN;
else if (strcmp(arg[iarg+1],"cosine") == 0) ordinate_style = COSINE;
else error->all(FLERR,"Illegal compute adf command");
iarg += 2;
} else error->all(FLERR,"Illegal compute adf command");
}
// triplewise args
if (!nargtriple) ntriples = 1;
else {
if (nargtriple % nargsperadf)
error->all(FLERR,"Illegal compute adf command");
ntriples = nargtriple/nargsperadf;
}
size_array_rows = nbin;
size_array_cols = 1 + 2*ntriples;
int ntypes = atom->ntypes;
memory->create(iatomflag,ntriples,ntypes+1,
"adf:iatomflag");
memory->create(jatomflag,ntriples,ntypes+1,
"adf:jatomflag");
memory->create(katomflag,ntriples,ntypes+1,
"adf:katomflag");
ilo = new int[ntriples];
ihi = new int[ntriples];
jlo = new int[ntriples];
jhi = new int[ntriples];
klo = new int[ntriples];
khi = new int[ntriples];
iatomcount = new int[ntriples];
iatomcountall = new int[ntriples];
rcutinnerj = new double[ntriples];
rcutouterj = new double[ntriples];
rcutinnerk = new double[ntriples];
rcutouterk = new double[ntriples];
if (!nargtriple) {
ilo[0] = 1; ihi[0] = ntypes;
jlo[0] = 1; jhi[0] = ntypes;
klo[0] = 1; khi[0] = ntypes;
} else {
cutflag = 1;
iarg = 4;
for (int itriple = 0; itriple < ntriples; itriple++) {
force->bounds(FLERR,arg[iarg],atom->ntypes,ilo[itriple],ihi[itriple]);
force->bounds(FLERR,arg[iarg+1],atom->ntypes,jlo[itriple],jhi[itriple]);
force->bounds(FLERR,arg[iarg+2],atom->ntypes,klo[itriple],khi[itriple]);
if (ilo[itriple] > ihi[itriple] ||
jlo[itriple] > jhi[itriple] ||
klo[itriple] > khi[itriple])
error->all(FLERR,"Illegal compute adf command");
rcutinnerj[itriple] = force->numeric(FLERR,arg[iarg+3]);
rcutouterj[itriple] = force->numeric(FLERR,arg[iarg+4]);
if (rcutinnerj[itriple] < 0.0 || rcutinnerj[itriple] >= rcutouterj[itriple])
error->all(FLERR,"Illegal compute adf command");
rcutinnerk[itriple] = force->numeric(FLERR,arg[iarg+5]);
rcutouterk[itriple] = force->numeric(FLERR,arg[iarg+6]);
if (rcutinnerk[itriple] < 0.0 || rcutinnerk[itriple] >= rcutouterk[itriple])
error->all(FLERR,"Illegal compute adf command");
iarg += nargsperadf;
}
}
// identify central atom types
int i,j,k;
for (int m = 0; m < ntriples; m++) {
for (i = 1; i <= ntypes; i++)
iatomflag[m][i] = 0;
for (i = ilo[m]; i <= ihi[m]; i++)
iatomflag[m][i] = 1;
}
// identify j atom types
for (int m = 0; m < ntriples; m++) {
for (j = 1; j <= ntypes; j++)
jatomflag[m][j] = 0;
for (j = jlo[m]; j <= jhi[m]; j++)
jatomflag[m][j] = 1;
}
// identify k atom types
for (int m = 0; m < ntriples; m++) {
for (k = 1; k <= ntypes; k++)
katomflag[m][k] = 0;
for (k = klo[m]; k <= khi[m]; k++)
katomflag[m][k] = 1;
}
memory->create(hist,ntriples,nbin,"adf:hist");
memory->create(histall,ntriples,nbin,"adf:histall");
memory->create(array,nbin,1+2*ntriples,"adf:array");
// initialize list of jatom neighbor lists
maxjatom = new int[ntriples];
numjatom = new int[ntriples];
neighjatom = new int*[ntriples];
for (int m = 0; m < ntriples; m++) {
maxjatom[m] = 10;
memory->create(neighjatom[m],maxjatom[m],"adf:neighjatom");
}
// initialize list of katom neighbor lists
maxkatom = new int[ntriples];
numkatom = new int[ntriples];
neighkatom = new int*[ntriples];
for (int m = 0; m < ntriples; m++) {
maxkatom[m] = 10;
memory->create(neighkatom[m],maxkatom[m],"adf:neighkatom");
}
// initialize list of short neighbor lists
maxjkatom = new int[ntriples];
numjkatom = new int[ntriples];
neighjkatom = new int*[ntriples];
bothjkatom = new int*[ntriples];
delrjkatom = new double**[ntriples];
for (int m = 0; m < ntriples; m++) {
maxjkatom[m] = 10;
memory->create(neighjkatom[m],maxjkatom[m],"adf:neighjkatom");
memory->create(bothjkatom[m],maxjkatom[m],"adf:bothjkatom");
memory->create(delrjkatom[m],maxjkatom[m],4,"adf:delrjkatom");
}
rad2deg = 180.0 / MY_PI;
}
/* ---------------------------------------------------------------------- */
ComputeADF::~ComputeADF()
{
memory->destroy(iatomflag);
delete [] ilo;
delete [] ihi;
delete [] jlo;
delete [] jhi;
delete [] klo;
delete [] khi;
delete [] iatomcount;
delete [] iatomcountall;
memory->destroy(hist);
memory->destroy(histall);
memory->destroy(array);
memory->destroy(jatomflag);
delete [] rcutinnerj;
delete [] rcutouterj;
delete [] maxjatom;
delete [] numjatom;
for (int m = 0; m < ntriples; m++)
memory->destroy(neighjatom[m]);
delete [] neighjatom;
memory->destroy(katomflag);
delete [] rcutinnerk;
delete [] rcutouterk;
delete [] maxkatom;
delete [] numkatom;
for (int m = 0; m < ntriples; m++)
memory->destroy(neighkatom[m]);
delete [] neighkatom;
delete [] maxjkatom;
delete [] numjkatom;
for (int m = 0; m < ntriples; m++)
memory->destroy(neighjkatom[m]);
delete [] neighjkatom;
for (int m = 0; m < ntriples; m++)
memory->destroy(bothjkatom[m]);
delete [] bothjkatom;
for (int m = 0; m < ntriples; m++)
memory->destroy(delrjkatom[m]);
delete [] delrjkatom;
}
/* ---------------------------------------------------------------------- */
void ComputeADF::init()
{
double mycutneigh;
double maxouter;
if (!cutflag) {
if (!force->pair)
error->all(FLERR,"Compute adf requires a pair style be defined "
"or outer cutoff specified");
rcutinnerj[0] = 0.0;
rcutinnerk[0] = 0.0;
rcutouterj[0] = force->pair->cutforce;
rcutouterk[0] = force->pair->cutforce;
maxouter = force->pair->cutforce;;
} else {
maxouter = 0.0;
for (int m = 0; m < ntriples; m++) {
maxouter = MAX(rcutouterj[m],maxouter);
maxouter = MAX(rcutouterk[m],maxouter);
}
}
// specify mycutneigh if force cutoff too small or non-existent
if (!(force->pair) || maxouter > force->pair->cutforce) {
double skin = neighbor->skin;
mycutneigh = maxouter + skin;
if (mycutneigh > comm->cutghostuser)
error->all(FLERR,"Compute adf outer cutoff exceeds ghost atom range - "
"use comm_modify cutoff command");
}
// assign ordinate values to 1st column of output array
int x0;
if (ordinate_style == DEGREE) {
deltax = MY_PI / nbin * rad2deg;
deltaxinv = nbin / MY_PI;
x0 = 0.0;
} else if (ordinate_style == RADIAN) {
deltax = MY_PI / nbin;
deltaxinv = nbin / MY_PI;
x0 = 0.0;
} else if (ordinate_style == COSINE) {
deltax = 2.0 / nbin;
deltaxinv = 1.0/deltax;
x0 = -1.0;
}
for (int i = 0; i < nbin; i++)
array[i][0] = x0 + (i+0.5) * deltax;
// need an occasional full neighbor list
// if mycutneigh specified, request a cutoff = maxouter + skin
// skin is included b/c Neighbor uses this value similar
// to its cutneighmax = force cutoff + skin
// also, this NeighList may be used by this compute for multiple steps
// (until next reneighbor), so it needs to contain atoms further
// than maxouter apart, just like a normal neighbor list does
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->occasional = 1;
if (mycutneigh > 0.0) {
neighbor->requests[irequest]->cut = 1;
neighbor->requests[irequest]->cutoff = mycutneigh;
}
}
/* ---------------------------------------------------------------------- */
void ComputeADF::init_list(int /*id*/, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeADF::compute_array()
{
int i,j,k,m,ii,jj,jatom,katom,jk,jjj,kkk;
int inum,jnum,itype,jtype,ibin,ihisto;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *ilist,*jlist,*klist,*numneigh,**firstneigh;
double factor_lj,factor_coul;
double delr1[3],delr2[3],rinv1,rinv2,rinv12,cs,theta;
invoked_array = update->ntimestep;
// 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;
// zero the histogram counts
for (i = 0; i < ntriples; i++)
for (j = 0; j < nbin; j++)
hist[i][j] = 0.0;
// zero the central atom counts
for (i = 0; i < ntriples; i++)
iatomcount[m] = 0;
// tally the ADFs
// all three atoms i, j, and k must be in fix group
// tally I,J,K triple only if I is central atom
// and J,K matches unordered neighbor types (JJ,KK)
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (!(mask[i] & groupbit)) continue;
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
// count atom i in each matching ADF
// zero the jatom, katom neighbor list counts
for (m = 0; m < ntriples; m++) {
if (iatomflag[m][itype]) {
iatomcount[m]++;
}
numjatom[m] = 0;
numkatom[m] = 0;
numjkatom[m] = 0;
}
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
// if both weighting factors are 0, skip this pair
// could be 0 and still be in neigh list for long-range Coulombics
// want consistency with non-charged triples which wouldn't be in list
if (factor_lj == 0.0 && factor_coul == 0.0) continue;
if (!(mask[j] & groupbit)) continue;
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;
for (m = 0; m < ntriples; m++) {
// check if itype, jtype, ktype match this ADF definition
// if yes, add j to jatom, katom, and jkatom lists
if (!iatomflag[m][itype]) continue;
int jflag = 0;
if (jatomflag[m][jtype] &&
rsq >= rcutinnerj[m]*rcutinnerj[m] &&
rsq <= rcutouterj[m]*rcutouterj[m]) {
jflag = 1;
jatom = numjatom[m]++;
neighjatom[m][jatom] = numjkatom[m];
if (numjatom[m] >= maxjatom[m]) {
maxjatom[m] += maxjatom[m]/2;
memory->grow(neighjatom[m],maxjatom[m],"adf:neighjatom");
}
}
int kflag = 0;
if (katomflag[m][jtype] &&
rsq >= rcutinnerk[m]*rcutinnerk[m] &&
rsq <= rcutouterk[m]*rcutouterk[m]) {
kflag = 1;
katom = numkatom[m]++;
neighkatom[m][katom] = numjkatom[m];
if (numkatom[m] >= maxkatom[m]) {
maxkatom[m] += maxkatom[m]/2;
memory->grow(neighkatom[m],maxkatom[m],"adf:neighkatom");
}
}
// if atom in either list, add to jklist
if (jflag || kflag) {
jk = numjkatom[m]++;
neighjkatom[m][jk] = j;
delrjkatom[m][jk][0] = delx;
delrjkatom[m][jk][1] = dely;
delrjkatom[m][jk][2] = delz;
delrjkatom[m][jk][3] = 1.0/sqrt(rsq);
if (numjkatom[m] >= maxjkatom[m]) {
maxjkatom[m] += maxjkatom[m]/2;
memory->grow(neighjkatom[m],maxjkatom[m],"adf:neighjkatom");
memory->grow(bothjkatom[m],maxjkatom[m],"adf:bothjkatom");
memory->grow(delrjkatom[m],maxjkatom[m],4,"adf:delrjkatom");
}
// indicate if atom in both lists
if (jflag && kflag)
bothjkatom[m][jk] = 1;
else
bothjkatom[m][jk] = 0;
}
}
}
// loop over ADFs
for (m = 0; m < ntriples; m++) {
// loop over (j,k) pairs
for (jatom = 0; jatom < numjatom[m]; jatom++) {
jjj = neighjatom[m][jatom];
j = neighjkatom[m][jjj];
delr1[0] = delrjkatom[m][jjj][0];
delr1[1] = delrjkatom[m][jjj][1];
delr1[2] = delrjkatom[m][jjj][2];
rinv1 = delrjkatom[m][jjj][3];
for (katom = 0; katom < numkatom[m]; katom++) {
kkk = neighkatom[m][katom];
k = neighjkatom[m][kkk];
// skip if j==k, or j > k and both are in both lists
if (k == j || (j > k && bothjkatom[m][jjj] && bothjkatom[m][kkk])) continue;
delr2[0] = delrjkatom[m][kkk][0];
delr2[1] = delrjkatom[m][kkk][1];
delr2[2] = delrjkatom[m][kkk][2];
rinv2 = delrjkatom[m][kkk][3];
rinv12 = rinv1*rinv2;
cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12;
if (ordinate_style == COSINE) {
ibin = static_cast<int> ((cs+1.0)*deltaxinv);
if (ibin >= nbin || ibin < 0) continue;
hist[m][ibin] += 1.0;
} else {
if (cs > 1.0) cs = 1.0;
if (cs < -1.0) cs = -1.0;
theta = acos(cs);
ibin = static_cast<int> (theta*deltaxinv);
if (ibin >= nbin || ibin < 0) continue;
hist[m][ibin] += 1.0;
}
}
}
}
}
// sum histograms across procs
MPI_Allreduce(hist[0],histall[0],ntriples*nbin,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(iatomcount,iatomcountall,ntriples,MPI_INT,MPI_SUM,world);
// convert counts to pdf(theta) and adf(theta)
// copy into output array
for (m = 0; m < ntriples; m++) {
double count = 0;
for (ibin = 0; ibin < nbin; ibin++)
count += histall[m][ibin];
double normfac1, pdftheta, normfac2, adftheta;
if (count > 0.0) normfac1 = deltaxinv/count;
else normfac1 = 0.0;
if (iatomcountall[m] > 0.0) normfac2 = 1.0/iatomcountall[m];
else normfac2 = 0.0;
adftheta = 0.0;
for (ibin = 0; ibin < nbin; ibin++) {
pdftheta = histall[m][ibin] * normfac1;
adftheta += histall[m][ibin] * normfac2;
array[ibin][1+2*m] = pdftheta;
array[ibin][2+2*m] = adftheta;
}
}
}

89
src/compute_adf.h Normal file
View File

@ -0,0 +1,89 @@
/* -*- 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(adf,ComputeADF)
#else
#ifndef LMP_COMPUTE_ADF_H
#define LMP_COMPUTE_ADF_H
#include <cstdio>
#include "compute.h"
namespace LAMMPS_NS {
class ComputeADF : public Compute {
public:
ComputeADF(class LAMMPS *, int, char **);
~ComputeADF();
void init();
void init_list(int, class NeighList *);
void compute_array();
private:
int nbin; // # of adf bins
int ntriples; // # of adf triples
double deltax,deltaxinv; // bin width and inverse-width
int *ilo,*ihi,*jlo,*jhi,*klo,*khi;
double **hist; // histogram bins
double **histall; // summed histogram bins across all procs
double *rcutinnerj, *rcutinnerk; // list of inner cutoffs
double *rcutouterj, *rcutouterk; // list of outer cutoffs
class NeighList *list; // full neighbor list
int *iatomcount; // local number of central atoms
int *iatomcountall; // total number of central atoms
int **iatomflag; // 1 if type is central atom in ADF
int *maxjatom, *maxkatom; // allocated size jatom, katom neighlist
int *numjatom, *numkatom; // actual size of jatom, katom neighlist
int **neighjatom, **neighkatom;// list of short neighbor lists
int **jatomflag, **katomflag; // 1 if type is neighbor atom in ADF
int *maxjkatom; // allocated size short neighlist
int *numjkatom; // actual size of short neighlist
int **neighjkatom; // list of short neighbor lists
int **bothjkatom; // 1 if atom is in both jatom and katom lists
double ***delrjkatom; // list of 4-vectors: delx, dely, delx, and 1/r
double rad2deg; // conversion factor from radians to degrees
int ordinate_style; // DEGREE, RADIAN, or COSINE
int cutflag; // 1 if at least one outer cutoff specified
};
}
#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 adf requires a pair style be defined or outer cutoff specified
Self-explanatory.
E: Compute adf cutoff exceeds ghost atom range - use comm_modify cutoff command
Self-explanatary.
*/