ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

convert more styles to use new neighbor list request API, apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-03-05 21:09:19 -05:00
parent ffb367663d
commit 54e667e491
20 changed files with 1083 additions and 1122 deletions
Download Patch File Download Diff File Expand all files Collapse all files

230
src/compute_centro_atom.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -26,7 +25,6 @@
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "pair.h"
#include "update.h"
@ -38,15 +36,16 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeCentroAtom::ComputeCentroAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
distsq(nullptr), nearest(nullptr), centro(nullptr)
Compute(lmp, narg, arg), distsq(nullptr), nearest(nullptr), centro(nullptr)
{
if (narg < 4 || narg > 6)
error->all(FLERR,"Illegal compute centro/atom command");
if (narg < 4 || narg > 6) error->all(FLERR, "Illegal compute centro/atom command");
if (strcmp(arg[3],"fcc") == 0) nnn = 12;
else if (strcmp(arg[3],"bcc") == 0) nnn = 8;
else nnn = utils::inumeric(FLERR,arg[3],false,lmp);
if (strcmp(arg[3], "fcc") == 0)
nnn = 12;
else if (strcmp(arg[3], "bcc") == 0)
nnn = 8;
else
nnn = utils::inumeric(FLERR, arg[3], false, lmp);
// default values
@ -56,19 +55,22 @@ ComputeCentroAtom::ComputeCentroAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"axes") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute centro/atom command3");
axes_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
if (strcmp(arg[iarg], "axes") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute centro/atom command3");
axes_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else error->all(FLERR,"Illegal compute centro/atom command1");
} else
error->all(FLERR, "Illegal compute centro/atom command1");
}
if (nnn <= 0 || nnn % 2)
error->all(FLERR,"Illegal neighbor value for compute centro/atom command");
error->all(FLERR, "Illegal neighbor value for compute centro/atom command");
peratom_flag = 1;
if (!axes_flag) size_peratom_cols = 0;
else size_peratom_cols = 10;
if (!axes_flag)
size_peratom_cols = 0;
else
size_peratom_cols = 10;
nmax = 0;
maxneigh = 0;
@ -89,22 +91,16 @@ ComputeCentroAtom::~ComputeCentroAtom()
void ComputeCentroAtom::init()
{
if (force->pair == nullptr)
error->all(FLERR,"Compute centro/atom requires a pair style be defined");
error->all(FLERR, "Compute centro/atom requires a pair style be defined");
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"centro/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute centro/atom");
if (strcmp(modify->compute[i]->style, "centro/atom") == 0) count++;
if (count > 1 && comm->me == 0) error->warning(FLERR, "More than one compute centro/atom");
// 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;
neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
}
/* ---------------------------------------------------------------------- */
@ -118,9 +114,9 @@ void ComputeCentroAtom::init_list(int /*id*/, NeighList *ptr)
void ComputeCentroAtom::compute_peratom()
{
int i,j,k,ii,jj,kk,n,inum,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,value;
int *ilist,*jlist,*numneigh,**firstneigh;
int i, j, k, ii, jj, kk, n, inum, jnum;
double xtmp, ytmp, ztmp, delx, dely, delz, rsq, value;
int *ilist, *jlist, *numneigh, **firstneigh;
invoked_peratom = update->ntimestep;
@ -131,15 +127,14 @@ void ComputeCentroAtom::compute_peratom()
if (!axes_flag) {
memory->destroy(centro);
nmax = atom->nmax;
memory->create(centro,nmax,"centro/atom:centro");
memory->create(centro, nmax, "centro/atom:centro");
vector_atom = centro;
} else {
memory->destroy(centro);
memory->destroy(array_atom);
nmax = atom->nmax;
memory->create(centro,nmax,"centro/atom:centro");
memory->create(array_atom,nmax,size_peratom_cols,
"centro/atom:array_atom");
memory->create(centro, nmax, "centro/atom:centro");
memory->create(array_atom, nmax, size_peratom_cols, "centro/atom:array_atom");
}
}
@ -154,8 +149,8 @@ void ComputeCentroAtom::compute_peratom()
// npairs = number of unique pairs
int nhalf = nnn/2;
int npairs = nnn * (nnn-1) / 2;
int nhalf = nnn / 2;
int npairs = nnn * (nnn - 1) / 2;
double *pairs = new double[npairs];
// compute centro-symmetry parameter for each atom in group
@ -180,8 +175,8 @@ void ComputeCentroAtom::compute_peratom()
memory->destroy(distsq);
memory->destroy(nearest);
maxneigh = jnum;
memory->create(distsq,maxneigh,"centro/atom:distsq");
memory->create(nearest,maxneigh,"centro/atom:nearest");
memory->create(distsq, maxneigh, "centro/atom:distsq");
memory->create(nearest, maxneigh, "centro/atom:nearest");
}
// loop over list of all neighbors within force cutoff
@ -196,7 +191,7 @@ void ComputeCentroAtom::compute_peratom()
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cutsq) {
distsq[n] = rsq;
nearest[n++] = j;
@ -216,7 +211,7 @@ void ComputeCentroAtom::compute_peratom()
// store nnn nearest neighs in 1st nnn locations of distsq and nearest
select2(nnn,n,distsq,nearest);
select2(nnn, n, distsq, nearest);
// R = Ri + Rj for each of npairs i,j pairs among nnn neighbors
// pairs = squared length of each R
@ -224,13 +219,12 @@ void ComputeCentroAtom::compute_peratom()
n = 0;
for (j = 0; j < nnn; j++) {
jj = nearest[j];
for (k = j+1; k < nnn; k++) {
for (k = j + 1; k < nnn; k++) {
kk = nearest[k];
delx = x[jj][0] + x[kk][0] - 2.0*xtmp;
dely = x[jj][1] + x[kk][1] - 2.0*ytmp;
delz = x[jj][2] + x[kk][2] - 2.0*ztmp;
pairs[n++] = delx*delx + dely*dely + delz*delz;
delx = x[jj][0] + x[kk][0] - 2.0 * xtmp;
dely = x[jj][1] + x[kk][1] - 2.0 * ytmp;
delz = x[jj][2] + x[kk][2] - 2.0 * ztmp;
pairs[n++] = delx * delx + dely * dely + delz * delz;
}
}
@ -242,11 +236,11 @@ void ComputeCentroAtom::compute_peratom()
// R1, R2 are corresponding vectors Ri - Rj
// R3 is normal to R1, R2
double rsq1,rsq2;
double rsq1, rsq2;
double* r1 = &array_atom[i][1];
double* r2 = &array_atom[i][4];
double* r3 = &array_atom[i][7];
double *r1 = &array_atom[i][1];
double *r2 = &array_atom[i][4];
double *r3 = &array_atom[i][7];
if (n < nnn) {
centro[i] = 0.0;
@ -258,18 +252,18 @@ void ComputeCentroAtom::compute_peratom()
// store nnn nearest neighs in 1st nnn locations of distsq and nearest
select2(nnn,n,distsq,nearest);
select2(nnn, n, distsq, nearest);
n = 0;
rsq1 = rsq2 = cutsq;
for (j = 0; j < nnn; j++) {
jj = nearest[j];
for (k = j+1; k < nnn; k++) {
for (k = j + 1; k < nnn; k++) {
kk = nearest[k];
delx = x[jj][0] + x[kk][0] - 2.0*xtmp;
dely = x[jj][1] + x[kk][1] - 2.0*ytmp;
delz = x[jj][2] + x[kk][2] - 2.0*ztmp;
rsq = delx*delx + dely*dely + delz*delz;
delx = x[jj][0] + x[kk][0] - 2.0 * xtmp;
dely = x[jj][1] + x[kk][1] - 2.0 * ytmp;
delz = x[jj][2] + x[kk][2] - 2.0 * ztmp;
rsq = delx * delx + dely * dely + delz * delz;
pairs[n++] = rsq;
if (rsq < rsq2) {
@ -277,16 +271,16 @@ void ComputeCentroAtom::compute_peratom()
rsq2 = rsq1;
MathExtra::copy3(r1, r2);
rsq1 = rsq;
MathExtra::sub3(x[jj],x[kk],r1);
MathExtra::sub3(x[jj], x[kk], r1);
} else {
rsq2 = rsq;
MathExtra::sub3(x[jj],x[kk],r2);
MathExtra::sub3(x[jj], x[kk], r2);
}
}
}
}
MathExtra::cross3(r1,r2,r3);
MathExtra::cross3(r1, r2, r3);
MathExtra::norm3(r1);
MathExtra::norm3(r2);
MathExtra::norm3(r3);
@ -294,7 +288,7 @@ void ComputeCentroAtom::compute_peratom()
// store nhalf smallest pair distances in 1st nhalf locations of pairs
select(nhalf,npairs,pairs);
select(nhalf, npairs, pairs);
// centrosymmetry = sum of nhalf smallest squared values
@ -312,64 +306,62 @@ void ComputeCentroAtom::compute_peratom()
}
}
delete [] pairs;
delete[] pairs;
if (axes_flag)
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit)
array_atom[i][0] = centro[i];
if (mask[i] & groupbit) array_atom[i][0] = centro[i];
}
}
/* ----------------------------------------------------------------------
2 select routines from Numerical Recipes (slightly modified)
find k smallest values in array of length n
2nd routine sorts auxiliary array at same time
------------------------------------------------------------------------- */
#define SWAP(a,b) tmp = a; a = b; b = tmp;
#define ISWAP(a,b) itmp = a; a = b; b = itmp;
#define SWAP(a, b) \
tmp = a; \
a = b; \
b = tmp;
#define ISWAP(a, b) \
itmp = a; \
a = b; \
b = itmp;
void ComputeCentroAtom::select(int k, int n, double *arr)
{
int i,ir,j,l,mid;
double a,tmp;
int i, ir, j, l, mid;
double a, tmp;
arr--;
l = 1;
ir = n;
for (;;) {
if (ir <= l+1) {
if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir])
}
if (ir <= l + 1) {
if (ir == l + 1 && arr[ir] < arr[l]) { SWAP(arr[l], arr[ir]) }
return;
} else {
mid=(l+ir) >> 1;
SWAP(arr[mid],arr[l+1])
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir])
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir])
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1])
}
i = l+1;
mid = (l + ir) >> 1;
SWAP(arr[mid], arr[l + 1])
if (arr[l] > arr[ir]) { SWAP(arr[l], arr[ir]) }
if (arr[l + 1] > arr[ir]) { SWAP(arr[l + 1], arr[ir]) }
if (arr[l] > arr[l + 1]) { SWAP(arr[l], arr[l + 1]) }
i = l + 1;
j = ir;
a = arr[l+1];
a = arr[l + 1];
for (;;) {
do i++; while (arr[i] < a);
do j--; while (arr[j] > a);
do i++;
while (arr[i] < a);
do j--;
while (arr[j] > a);
if (j < i) break;
SWAP(arr[i],arr[j])
SWAP(arr[i], arr[j])
}
arr[l+1] = arr[j];
arr[l + 1] = arr[j];
arr[j] = a;
if (j >= k) ir = j-1;
if (j >= k) ir = j - 1;
if (j <= k) l = i;
}
}
@ -379,52 +371,54 @@ void ComputeCentroAtom::select(int k, int n, double *arr)
void ComputeCentroAtom::select2(int k, int n, double *arr, int *iarr)
{
int i,ir,j,l,mid,ia,itmp;
double a,tmp;
int i, ir, j, l, mid, ia, itmp;
double a, tmp;
arr--;
iarr--;
l = 1;
ir = n;
for (;;) {
if (ir <= l+1) {
if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir])
ISWAP(iarr[l],iarr[ir])
if (ir <= l + 1) {
if (ir == l + 1 && arr[ir] < arr[l]) {
SWAP(arr[l], arr[ir])
ISWAP(iarr[l], iarr[ir])
}
return;
} else {
mid=(l+ir) >> 1;
SWAP(arr[mid],arr[l+1])
ISWAP(iarr[mid],iarr[l+1])
mid = (l + ir) >> 1;
SWAP(arr[mid], arr[l + 1])
ISWAP(iarr[mid], iarr[l + 1])
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir])
ISWAP(iarr[l],iarr[ir])
SWAP(arr[l], arr[ir])
ISWAP(iarr[l], iarr[ir])
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir])
ISWAP(iarr[l+1],iarr[ir])
if (arr[l + 1] > arr[ir]) {
SWAP(arr[l + 1], arr[ir])
ISWAP(iarr[l + 1], iarr[ir])
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1])
ISWAP(iarr[l],iarr[l+1])
if (arr[l] > arr[l + 1]) {
SWAP(arr[l], arr[l + 1])
ISWAP(iarr[l], iarr[l + 1])
}
i = l+1;
i = l + 1;
j = ir;
a = arr[l+1];
ia = iarr[l+1];
a = arr[l + 1];
ia = iarr[l + 1];
for (;;) {
do i++; while (arr[i] < a);
do j--; while (arr[j] > a);
do i++;
while (arr[i] < a);
do j--;
while (arr[j] > a);
if (j < i) break;
SWAP(arr[i],arr[j])
ISWAP(iarr[i],iarr[j])
SWAP(arr[i], arr[j])
ISWAP(iarr[i], iarr[j])
}
arr[l+1] = arr[j];
arr[l + 1] = arr[j];
arr[j] = a;
iarr[l+1] = iarr[j];
iarr[l + 1] = iarr[j];
iarr[j] = ia;
if (j >= k) ir = j-1;
if (j >= k) ir = j - 1;
if (j <= k) l = i;
}
}
@ -436,7 +430,7 @@ void ComputeCentroAtom::select2(int k, int n, double *arr, int *iarr)
double ComputeCentroAtom::memory_usage()
{
double bytes = (double)nmax * sizeof(double);
if (axes_flag) bytes += (double)size_peratom_cols*nmax * sizeof(double);
double bytes = (double) nmax * sizeof(double);
if (axes_flag) bytes += (double) size_peratom_cols * nmax * sizeof(double);
return bytes;
}