diff --git a/src/KOKKOS/compute_orientorder_atom_kokkos.cpp b/src/KOKKOS/compute_orientorder_atom_kokkos.cpp index f7ec879e3c..285f212c1a 100644 --- a/src/KOKKOS/compute_orientorder_atom_kokkos.cpp +++ b/src/KOKKOS/compute_orientorder_atom_kokkos.cpp @@ -85,15 +85,10 @@ void ComputeOrientOrderAtomKokkos::init() // need an occasional full neighbor list - // irequest = neigh request made by parent class - - int irequest = neighbor->nrequest - 1; - - neighbor->requests[irequest]-> - kokkos_host = std::is_same::value && - !std::is_same::value; - neighbor->requests[irequest]-> - kokkos_device = std::is_same::value; + auto request = neighbor->find_request(this); + request->set_kokkos_host(std::is_same::value && + !std::is_same::value); + request->set_kokkos_device(std::is_same::value); } /* ---------------------------------------------------------------------- */ diff --git a/src/compute_aggregate_atom.cpp b/src/compute_aggregate_atom.cpp index a378baf3ac..6811ae2d87 100644 --- a/src/compute_aggregate_atom.cpp +++ b/src/compute_aggregate_atom.cpp @@ -26,7 +26,6 @@ #include "memory.h" #include "modify.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include "pair.h" #include "update.h" diff --git a/src/compute_centro_atom.cpp b/src/compute_centro_atom.cpp index 90f535613b..535720c575 100644 --- a/src/compute_centro_atom.cpp +++ b/src/compute_centro_atom.cpp @@ -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; } diff --git a/src/compute_cluster_atom.cpp b/src/compute_cluster_atom.cpp index c429ba5cac..07a0c700ee 100644 --- a/src/compute_cluster_atom.cpp +++ b/src/compute_cluster_atom.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -22,7 +21,6 @@ #include "memory.h" #include "modify.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include "pair.h" #include "update.h" @@ -32,18 +30,17 @@ using namespace LAMMPS_NS; -enum{CLUSTER,MASK,COORDS}; +enum { CLUSTER, MASK, COORDS }; /* ---------------------------------------------------------------------- */ ComputeClusterAtom::ComputeClusterAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - clusterID(nullptr) + Compute(lmp, narg, arg), clusterID(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute cluster/atom command"); + if (narg != 4) error->all(FLERR, "Illegal compute cluster/atom command"); - double cutoff = utils::numeric(FLERR,arg[3],false,lmp); - cutsq = cutoff*cutoff; + double cutoff = utils::numeric(FLERR, arg[3], false, lmp); + cutsq = cutoff * cutoff; peratom_flag = 1; size_peratom_cols = 0; @@ -64,28 +61,21 @@ ComputeClusterAtom::~ComputeClusterAtom() void ComputeClusterAtom::init() { if (atom->tag_enable == 0) - error->all(FLERR,"Cannot use compute cluster/atom unless atoms have IDs"); + error->all(FLERR, "Cannot use compute cluster/atom unless atoms have IDs"); if (force->pair == nullptr) - error->all(FLERR,"Compute cluster/atom requires a pair style to be defined"); + error->all(FLERR, "Compute cluster/atom requires a pair style to be defined"); if (sqrt(cutsq) > force->pair->cutforce) - error->all(FLERR, - "Compute cluster/atom cutoff is longer than pairwise cutoff"); + error->all(FLERR, "Compute cluster/atom cutoff is longer than pairwise cutoff"); // need an occasional full neighbor list // full required so that pair of atoms on 2 procs both set their clusterID - 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); int count = 0; for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style,"cluster/atom") == 0) count++; - if (count > 1 && comm->me == 0) - error->warning(FLERR,"More than one compute cluster/atom"); + if (strcmp(modify->compute[i]->style, "cluster/atom") == 0) count++; + if (count > 1 && comm->me == 0) error->warning(FLERR, "More than one compute cluster/atom"); } /* ---------------------------------------------------------------------- */ @@ -99,9 +89,9 @@ void ComputeClusterAtom::init_list(int /*id*/, NeighList *ptr) void ComputeClusterAtom::compute_peratom() { - int i,j,ii,jj,inum,jnum; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum; + double xtmp, ytmp, ztmp, delx, dely, delz, rsq; + int *ilist, *jlist, *numneigh, **firstneigh; invoked_peratom = update->ntimestep; @@ -110,15 +100,17 @@ void ComputeClusterAtom::compute_peratom() if (atom->nmax > nmax) { memory->destroy(clusterID); nmax = atom->nmax; - memory->create(clusterID,nmax,"cluster/atom:clusterID"); + memory->create(clusterID, nmax, "cluster/atom:clusterID"); vector_atom = clusterID; } // invoke full neighbor list (will copy or build if necessary) // on the first step of a run, set preflag to one in neighbor->build_one(...) - if (update->firststep == update->ntimestep) neighbor->build_one(list,1); - else neighbor->build_one(list); + if (update->firststep == update->ntimestep) + neighbor->build_one(list, 1); + else + neighbor->build_one(list); inum = list->inum; ilist = list->ilist; @@ -148,8 +140,10 @@ void ComputeClusterAtom::compute_peratom() for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - if (mask[i] & groupbit) clusterID[i] = tag[i]; - else clusterID[i] = 0; + if (mask[i] & groupbit) + clusterID[i] = tag[i]; + else + clusterID[i] = 0; } // loop until no more changes on any proc: @@ -162,7 +156,7 @@ void ComputeClusterAtom::compute_peratom() commflag = CLUSTER; double **x = atom->x; - int change,done,anychange; + int change, done, anychange; while (true) { comm->forward_comm(this); @@ -189,9 +183,9 @@ void ComputeClusterAtom::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) { - clusterID[i] = clusterID[j] = MIN(clusterID[i],clusterID[j]); + clusterID[i] = clusterID[j] = MIN(clusterID[i], clusterID[j]); done = 0; } } @@ -202,17 +196,17 @@ void ComputeClusterAtom::compute_peratom() // stop if all procs are done - MPI_Allreduce(&change,&anychange,1,MPI_INT,MPI_MAX,world); + MPI_Allreduce(&change, &anychange, 1, MPI_INT, MPI_MAX, world); if (!anychange) break; } } /* ---------------------------------------------------------------------- */ -int ComputeClusterAtom::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int ComputeClusterAtom::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { - int i,j,m; + int i, j, m; m = 0; if (commflag == CLUSTER) { @@ -243,7 +237,7 @@ int ComputeClusterAtom::pack_forward_comm(int n, int *list, double *buf, void ComputeClusterAtom::unpack_forward_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; m = 0; last = first + n; @@ -268,6 +262,6 @@ void ComputeClusterAtom::unpack_forward_comm(int n, int first, double *buf) double ComputeClusterAtom::memory_usage() { - double bytes = (double)nmax * sizeof(double); + double bytes = (double) nmax * sizeof(double); return bytes; } diff --git a/src/compute_cna_atom.cpp b/src/compute_cna_atom.cpp index e002e54983..4b9c28362c 100644 --- a/src/compute_cna_atom.cpp +++ b/src/compute_cna_atom.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -25,36 +24,34 @@ #include "memory.h" #include "modify.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include "pair.h" #include "update.h" -#include #include +#include using namespace LAMMPS_NS; -#define MAXNEAR 16 -#define MAXCOMMON 8 +static constexpr int MAXNEAR = 16; +static constexpr int MAXCOMMON = 8; -enum{UNKNOWN,FCC,HCP,BCC,ICOS,OTHER}; -enum{NCOMMON,NBOND,MAXBOND,MINBOND}; +enum { UNKNOWN, FCC, HCP, BCC, ICOS, OTHER }; +enum { NCOMMON, NBOND, MAXBOND, MINBOND }; /* ---------------------------------------------------------------------- */ ComputeCNAAtom::ComputeCNAAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - list(nullptr), nearest(nullptr), nnearest(nullptr), pattern(nullptr) + Compute(lmp, narg, arg), list(nullptr), nearest(nullptr), nnearest(nullptr), pattern(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute cna/atom command"); + if (narg != 4) error->all(FLERR, "Illegal compute cna/atom command"); peratom_flag = 1; size_peratom_cols = 0; - double cutoff = utils::numeric(FLERR,arg[3],false,lmp); - if (cutoff < 0.0) error->all(FLERR,"Illegal compute cna/atom command"); - cutsq = cutoff*cutoff; + double cutoff = utils::numeric(FLERR, arg[3], false, lmp); + if (cutoff < 0.0) error->all(FLERR, "Illegal compute cna/atom command"); + cutsq = cutoff * cutoff; nmax = 0; } @@ -73,31 +70,25 @@ ComputeCNAAtom::~ComputeCNAAtom() void ComputeCNAAtom::init() { if (force->pair == nullptr) - error->all(FLERR,"Compute cna/atom requires a pair style be defined"); + error->all(FLERR, "Compute cna/atom requires a pair style be defined"); if (sqrt(cutsq) > force->pair->cutforce) - error->all(FLERR,"Compute cna/atom cutoff is longer than pairwise cutoff"); + error->all(FLERR, "Compute cna/atom cutoff is longer than pairwise cutoff"); // cannot use neighbor->cutneighmax b/c neighbor has not yet been init - if (2.0*sqrt(cutsq) > force->pair->cutforce + neighbor->skin && - comm->me == 0) - error->warning(FLERR,"Compute cna/atom cutoff may be too large to find " + if (2.0 * sqrt(cutsq) > force->pair->cutforce + neighbor->skin && comm->me == 0) + error->warning(FLERR, + "Compute cna/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,"cna/atom") == 0) count++; - if (count > 1 && comm->me == 0) - error->warning(FLERR,"More than one compute cna/atom defined"); + if (strcmp(modify->compute[i]->style, "cna/atom") == 0) count++; + if (count > 1 && comm->me == 0) error->warning(FLERR, "More than one compute cna/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; + neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL); } /* ---------------------------------------------------------------------- */ @@ -111,13 +102,13 @@ void ComputeCNAAtom::init_list(int /*id*/, NeighList *ptr) void ComputeCNAAtom::compute_peratom() { - int i,j,k,ii,jj,kk,m,n,inum,jnum,inear,jnear; - int firstflag,ncommon,nbonds,maxbonds,minbonds; - int nfcc,nhcp,nbcc4,nbcc6,nico,cj,ck,cl,cm; - int *ilist,*jlist,*numneigh,**firstneigh; - int cna[MAXNEAR][4],onenearest[MAXNEAR]; - int common[MAXCOMMON],bonds[MAXCOMMON]; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int i, j, k, ii, jj, kk, m, n, inum, jnum, inear, jnear; + int firstflag, ncommon, nbonds, maxbonds, minbonds; + int nfcc, nhcp, nbcc4, nbcc6, nico, cj, ck, cl, cm; + int *ilist, *jlist, *numneigh, **firstneigh; + int cna[MAXNEAR][4], onenearest[MAXNEAR]; + int common[MAXCOMMON], bonds[MAXCOMMON]; + double xtmp, ytmp, ztmp, delx, dely, delz, rsq; invoked_peratom = update->ntimestep; @@ -129,9 +120,9 @@ void ComputeCNAAtom::compute_peratom() memory->destroy(pattern); nmax = atom->nmax; - memory->create(nearest,nmax,MAXNEAR,"cna:nearest"); - memory->create(nnearest,nmax,"cna:nnearest"); - memory->create(pattern,nmax,"cna:cna_pattern"); + memory->create(nearest, nmax, MAXNEAR, "cna:nearest"); + memory->create(nnearest, nmax, "cna:nnearest"); + memory->create(pattern, nmax, "cna:cna_pattern"); vector_atom = pattern; } @@ -170,9 +161,10 @@ void ComputeCNAAtom::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) { - if (n < MAXNEAR) nearest[i][n++] = j; + if (n < MAXNEAR) + nearest[i][n++] = j; else { nerror++; break; @@ -185,9 +177,9 @@ void ComputeCNAAtom::compute_peratom() // warning message int nerrorall; - MPI_Allreduce(&nerror,&nerrorall,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&nerror, &nerrorall, 1, MPI_INT, MPI_SUM, world); if (nerrorall && comm->me == 0) - error->warning(FLERR,"Too many neighbors in CNA for {} atoms",nerrorall); + error->warning(FLERR, "Too many neighbors in CNA for {} atoms", nerrorall); // compute CNA for each atom in group // only performed if # of nearest neighbors = 12 or 14 (fcc,hcp) @@ -226,7 +218,8 @@ void ComputeCNAAtom::compute_peratom() 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]; + if (ncommon < MAXCOMMON) + common[ncommon++] = nearest[i][inear]; else if (firstflag) { nerror++; firstflag = 0; @@ -249,10 +242,12 @@ void ComputeCNAAtom::compute_peratom() delx = xtmp - x[k][0]; dely = ytmp - x[k][1]; delz = ztmp - x[k][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; if (rsq < cutsq) { - if (n < MAXNEAR) onenearest[n++] = k; - else break; + if (n < MAXNEAR) + onenearest[n++] = k; + else + break; } } @@ -261,7 +256,8 @@ void ComputeCNAAtom::compute_peratom() 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]; + if (ncommon < MAXCOMMON) + common[ncommon++] = nearest[i][inear]; else if (firstflag) { nerror++; firstflag = 0; @@ -278,17 +274,17 @@ void ComputeCNAAtom::compute_peratom() for (n = 0; n < ncommon; n++) bonds[n] = 0; nbonds = 0; - for (jj = 0; jj < ncommon-1; jj++) { + for (jj = 0; jj < ncommon - 1; jj++) { j = common[jj]; xtmp = x[j][0]; ytmp = x[j][1]; ztmp = x[j][2]; - for (kk = jj+1; kk < ncommon; kk++) { + for (kk = jj + 1; kk < ncommon; kk++) { k = common[kk]; delx = xtmp - x[k][0]; dely = ytmp - x[k][1]; delz = ztmp - x[k][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; if (rsq < cutsq) { nbonds++; bonds[jj]++; @@ -302,8 +298,8 @@ void ComputeCNAAtom::compute_peratom() maxbonds = 0; minbonds = MAXCOMMON; for (n = 0; n < ncommon; n++) { - maxbonds = MAX(bonds[n],maxbonds); - minbonds = MIN(bonds[n],minbonds); + maxbonds = MAX(bonds[n], maxbonds); + minbonds = MIN(bonds[n], minbonds); } cna[m][MAXBOND] = maxbonds; cna[m][MINBOND] = minbonds; @@ -320,13 +316,19 @@ void ComputeCNAAtom::compute_peratom() ck = cna[inear][NBOND]; cl = cna[inear][MAXBOND]; cm = cna[inear][MINBOND]; - if (cj == 4 && ck == 2 && cl == 1 && cm == 1) nfcc++; - else if (cj == 4 && ck == 2 && cl == 2 && cm == 0) nhcp++; - else if (cj == 5 && ck == 5 && cl == 2 && cm == 2) nico++; + if (cj == 4 && ck == 2 && cl == 1 && cm == 1) + nfcc++; + else if (cj == 4 && ck == 2 && cl == 2 && cm == 0) + nhcp++; + else if (cj == 5 && ck == 5 && cl == 2 && cm == 2) + nico++; } - if (nfcc == 12) pattern[i] = FCC; - else if (nfcc == 6 && nhcp == 6) pattern[i] = HCP; - else if (nico == 12) pattern[i] = ICOS; + if (nfcc == 12) + pattern[i] = FCC; + else if (nfcc == 6 && nhcp == 6) + pattern[i] = HCP; + else if (nico == 12) + pattern[i] = ICOS; } else if (nnearest[i] == 14) { for (inear = 0; inear < 14; inear++) { @@ -334,8 +336,10 @@ void ComputeCNAAtom::compute_peratom() ck = cna[inear][NBOND]; cl = cna[inear][MAXBOND]; cm = cna[inear][MINBOND]; - if (cj == 4 && ck == 4 && cl == 2 && cm == 2) nbcc4++; - else if (cj == 6 && ck == 6 && cl == 2 && cm == 2) nbcc6++; + if (cj == 4 && ck == 4 && cl == 2 && cm == 2) + nbcc4++; + else if (cj == 6 && ck == 6 && cl == 2 && cm == 2) + nbcc6++; } if (nbcc4 == 6 && nbcc6 == 8) pattern[i] = BCC; } @@ -343,9 +347,9 @@ void ComputeCNAAtom::compute_peratom() // warning message - MPI_Allreduce(&nerror,&nerrorall,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&nerror, &nerrorall, 1, MPI_INT, MPI_SUM, world); if (nerrorall && comm->me == 0) - error->warning(FLERR,"Too many common neighbors in CNA: {}x", nerrorall); + error->warning(FLERR, "Too many common neighbors in CNA: {}x", nerrorall); } /* ---------------------------------------------------------------------- @@ -354,8 +358,8 @@ void ComputeCNAAtom::compute_peratom() double ComputeCNAAtom::memory_usage() { - double bytes = (double)nmax * sizeof(int); - bytes += (double)nmax * MAXNEAR * sizeof(int); - bytes += (double)nmax * sizeof(double); + double bytes = (double) nmax * sizeof(int); + bytes += (double) nmax * MAXNEAR * sizeof(int); + bytes += (double) nmax * sizeof(double); return bytes; } diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp index 45784205d3..4487e82985 100644 --- a/src/compute_coord_atom.cpp +++ b/src/compute_coord_atom.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -23,7 +22,6 @@ #include "memory.h" #include "modify.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include "pair.h" #include "update.h" @@ -33,36 +31,33 @@ using namespace LAMMPS_NS; - /* ---------------------------------------------------------------------- */ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - typelo(nullptr), typehi(nullptr), cvec(nullptr), carray(nullptr), - group2(nullptr), id_orientorder(nullptr), normv(nullptr) + Compute(lmp, narg, arg), typelo(nullptr), typehi(nullptr), cvec(nullptr), carray(nullptr), + group2(nullptr), id_orientorder(nullptr), normv(nullptr) { - if (narg < 5) error->all(FLERR,"Illegal compute coord/atom command"); + if (narg < 5) error->all(FLERR, "Illegal compute coord/atom command"); jgroup = group->find("all"); jgroupbit = group->bitmask[jgroup]; cstyle = NONE; - if (strcmp(arg[3],"cutoff") == 0) { + if (strcmp(arg[3], "cutoff") == 0) { cstyle = CUTOFF; - double cutoff = utils::numeric(FLERR,arg[4],false,lmp); - cutsq = cutoff*cutoff; + double cutoff = utils::numeric(FLERR, arg[4], false, lmp); + cutsq = cutoff * cutoff; int iarg = 5; - if ((narg > 6) && (strcmp(arg[5],"group") == 0)) { + if ((narg > 6) && (strcmp(arg[5], "group") == 0)) { group2 = utils::strdup(arg[6]); iarg += 2; jgroup = group->find(group2); - if (jgroup == -1) - error->all(FLERR,"Compute coord/atom group2 ID does not exist"); + if (jgroup == -1) error->all(FLERR, "Compute coord/atom group2 ID does not exist"); jgroupbit = group->bitmask[jgroup]; } - ncol = narg-iarg + 1; + ncol = narg - iarg + 1; int ntypes = atom->ntypes; typelo = new int[ncol]; typehi = new int[ncol]; @@ -74,29 +69,27 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) : } else { ncol = 0; while (iarg < narg) { - utils::bounds(FLERR,arg[iarg],1,ntypes,typelo[ncol],typehi[ncol],error); - if (typelo[ncol] > typehi[ncol]) - error->all(FLERR,"Illegal compute coord/atom command"); + utils::bounds(FLERR, arg[iarg], 1, ntypes, typelo[ncol], typehi[ncol], error); + if (typelo[ncol] > typehi[ncol]) error->all(FLERR, "Illegal compute coord/atom command"); ncol++; iarg++; } } - } else if (strcmp(arg[3],"orientorder") == 0) { + } else if (strcmp(arg[3], "orientorder") == 0) { cstyle = ORIENT; - if (narg != 6) error->all(FLERR,"Illegal compute coord/atom command"); + if (narg != 6) error->all(FLERR, "Illegal compute coord/atom command"); id_orientorder = utils::strdup(arg[4]); int iorientorder = modify->find_compute(id_orientorder); - if (iorientorder < 0) - error->all(FLERR,"Could not find compute coord/atom compute ID"); - if (!utils::strmatch(modify->compute[iorientorder]->style,"^orientorder/atom")) - error->all(FLERR,"Compute coord/atom compute ID is not orientorder/atom"); + if (iorientorder < 0) error->all(FLERR, "Could not find compute coord/atom compute ID"); + if (!utils::strmatch(modify->compute[iorientorder]->style, "^orientorder/atom")) + error->all(FLERR, "Compute coord/atom compute ID is not orientorder/atom"); - threshold = utils::numeric(FLERR,arg[5],false,lmp); + threshold = utils::numeric(FLERR, arg[5], false, lmp); if (threshold <= -1.0 || threshold >= 1.0) - error->all(FLERR,"Compute coord/atom threshold not between -1 and 1"); + error->all(FLERR, "Compute coord/atom threshold not between -1 and 1"); ncol = 1; typelo = new int[ncol]; @@ -104,11 +97,14 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) : typelo[0] = 1; typehi[0] = atom->ntypes; - } else error->all(FLERR,"Invalid cstyle in compute coord/atom"); + } 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; + if (ncol == 1) + size_peratom_cols = 0; + else + size_peratom_cols = ncol; nmax = 0; } @@ -119,12 +115,12 @@ ComputeCoordAtom::~ComputeCoordAtom() { if (copymode) return; - delete [] group2; - delete [] typelo; - delete [] typehi; + delete[] group2; + delete[] typelo; + delete[] typehi; memory->destroy(cvec); memory->destroy(carray); - delete [] id_orientorder; + delete[] id_orientorder; } /* ---------------------------------------------------------------------- */ @@ -133,30 +129,24 @@ void ComputeCoordAtom::init() { if (cstyle == ORIENT) { int iorientorder = modify->find_compute(id_orientorder); - c_orientorder = (ComputeOrientOrderAtom*)(modify->compute[iorientorder]); + 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*(2*l+1); + comm_forward = 2 * (2 * l + 1); if (!(c_orientorder->qlcompflag)) - error->all(FLERR,"Compute coord/atom requires components " - "option in compute orientorder/atom"); + error->all(FLERR, + "Compute coord/atom requires components option in compute orientorder/atom"); } if (force->pair == nullptr) - error->all(FLERR,"Compute coord/atom requires a pair style be defined"); + error->all(FLERR, "Compute coord/atom requires a pair style be defined"); if (sqrt(cutsq) > force->pair->cutforce) - error->all(FLERR, - "Compute coord/atom cutoff is longer than pairwise cutoff"); + error->all(FLERR, "Compute coord/atom cutoff is longer than pairwise cutoff"); // 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); } /* ---------------------------------------------------------------------- */ @@ -170,9 +160,9 @@ void ComputeCoordAtom::init_list(int /*id*/, NeighList *ptr) void ComputeCoordAtom::compute_peratom() { - int i,j,m,ii,jj,inum,jnum,jtype,n; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, m, ii, jj, inum, jnum, jtype, n; + double xtmp, ytmp, ztmp, delx, dely, delz, rsq; + int *ilist, *jlist, *numneigh, **firstneigh; double *count; invoked_peratom = update->ntimestep; @@ -183,12 +173,12 @@ void ComputeCoordAtom::compute_peratom() if (ncol == 1) { memory->destroy(cvec); nmax = atom->nmax; - memory->create(cvec,nmax,"coord/atom:cvec"); + memory->create(cvec, nmax, "coord/atom:cvec"); vector_atom = cvec; } else { memory->destroy(carray); nmax = atom->nmax; - memory->create(carray,nmax,ncol,"coord/atom:carray"); + memory->create(carray, nmax, ncol, "coord/atom:carray"); array_atom = carray; } } @@ -242,14 +232,14 @@ void ComputeCoordAtom::compute_peratom() 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 && jtype >= typelo[0] && jtype <= typehi[0]) - n++; + rsq = delx * delx + dely * dely + delz * delz; + if (rsq < cutsq && jtype >= typelo[0] && jtype <= typehi[0]) n++; } } cvec[i] = n; - } else cvec[i] = 0.0; + } else + cvec[i] = 0.0; } } else { @@ -265,7 +255,6 @@ void ComputeCoordAtom::compute_peratom() jlist = firstneigh[i]; jnum = numneigh[i]; - for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; @@ -274,11 +263,10 @@ void ComputeCoordAtom::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) { for (m = 0; m < ncol; m++) - if (jtype >= typelo[m] && jtype <= typehi[m]) - count[m] += 1.0; + if (jtype >= typelo[m] && jtype <= typehi[m]) count[m] += 1.0; } } } @@ -303,31 +291,30 @@ void ComputeCoordAtom::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) { double dot_product = 0.0; - for (m=0; m < 2*(2*l+1); m++) { - dot_product += normv[i][nqlist+m]*normv[j][nqlist+m]; + for (m = 0; m < 2 * (2 * l + 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; + } else + cvec[i] = 0.0; } } } /* ---------------------------------------------------------------------- */ -int ComputeCoordAtom::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int ComputeCoordAtom::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { - int i,m=0,j; + int i, m = 0, j; for (i = 0; i < n; ++i) { - for (j = nqlist; j < nqlist + 2*(2*l+1); ++j) { - buf[m++] = normv[list[i]][j]; - } + for (j = nqlist; j < nqlist + 2 * (2 * l + 1); ++j) { buf[m++] = normv[list[i]][j]; } } return m; @@ -337,12 +324,10 @@ int ComputeCoordAtom::pack_forward_comm(int n, int *list, double *buf, void ComputeCoordAtom::unpack_forward_comm(int n, int first, double *buf) { - int i,last,m=0,j; + int i, last, m = 0, j; last = first + n; for (i = first; i < last; ++i) { - for (j = nqlist; j < nqlist + 2*(2*l+1); ++j) { - normv[i][j] = buf[m++]; - } + for (j = nqlist; j < nqlist + 2 * (2 * l + 1); ++j) { normv[i][j] = buf[m++]; } } } @@ -352,6 +337,6 @@ void ComputeCoordAtom::unpack_forward_comm(int n, int first, double *buf) double ComputeCoordAtom::memory_usage() { - double bytes = (double)ncol*nmax * sizeof(double); + double bytes = (double) ncol * nmax * sizeof(double); return bytes; } diff --git a/src/compute_group_group.cpp b/src/compute_group_group.cpp index db26a75880..2fcd3b9c22 100644 --- a/src/compute_group_group.cpp +++ b/src/compute_group_group.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -28,28 +27,26 @@ #include "kspace.h" #include "math_const.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include "pair.h" #include "update.h" -#include #include +#include using namespace LAMMPS_NS; using namespace MathConst; #define SMALL 0.00001 -enum{OFF,INTER,INTRA}; +enum { OFF, INTER, INTRA }; /* ---------------------------------------------------------------------- */ ComputeGroupGroup::ComputeGroupGroup(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - group2(nullptr) + Compute(lmp, narg, arg), group2(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute group/group command"); + if (narg < 4) error->all(FLERR, "Illegal compute group/group command"); scalar_flag = vector_flag = 1; size_vector = 3; @@ -58,8 +55,7 @@ ComputeGroupGroup::ComputeGroupGroup(LAMMPS *lmp, int narg, char **arg) : group2 = utils::strdup(arg[3]); jgroup = group->find(group2); - if (jgroup == -1) - error->all(FLERR,"Compute group/group group ID does not exist"); + if (jgroup == -1) error->all(FLERR, "Compute group/group group ID does not exist"); jgroupbit = group->bitmask[jgroup]; pairflag = 1; @@ -69,28 +65,33 @@ ComputeGroupGroup::ComputeGroupGroup(LAMMPS *lmp, int narg, char **arg) : int iarg = 4; while (iarg < narg) { - if (strcmp(arg[iarg],"pair") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute group/group command"); - pairflag = utils::logical(FLERR,arg[iarg+1],false,lmp); + if (strcmp(arg[iarg], "pair") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute group/group command"); + pairflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"kspace") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute group/group command"); - kspaceflag = utils::logical(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "kspace") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute group/group command"); + kspaceflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"boundary") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute group/group command"); - boundaryflag = utils::logical(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "boundary") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute group/group command"); + boundaryflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"molecule") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute group/group command"); - if (strcmp(arg[iarg+1],"off") == 0) molflag = OFF; - else if (strcmp(arg[iarg+1],"inter") == 0) molflag = INTER; - else if (strcmp(arg[iarg+1],"intra") == 0) molflag = INTRA; - else error->all(FLERR,"Illegal compute group/group command"); + } else if (strcmp(arg[iarg], "molecule") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute group/group command"); + if (strcmp(arg[iarg + 1], "off") == 0) + molflag = OFF; + else if (strcmp(arg[iarg + 1], "inter") == 0) + molflag = INTER; + else if (strcmp(arg[iarg + 1], "intra") == 0) + molflag = INTRA; + else + error->all(FLERR, "Illegal compute group/group command"); if (molflag != OFF && atom->molecule_flag == 0) - error->all(FLERR,"Compute group/group molecule requires molecule IDs"); + error->all(FLERR, "Compute group/group molecule requires molecule IDs"); iarg += 2; - } else error->all(FLERR,"Illegal compute group/group command"); + } else + error->all(FLERR, "Illegal compute group/group command"); } vector = new double[size_vector]; @@ -100,8 +101,8 @@ ComputeGroupGroup::ComputeGroupGroup(LAMMPS *lmp, int narg, char **arg) : ComputeGroupGroup::~ComputeGroupGroup() { - delete [] group2; - delete [] vector; + delete[] group2; + delete[] vector; } /* ---------------------------------------------------------------------- */ @@ -112,50 +113,47 @@ void ComputeGroupGroup::init() // if hybrid, let hybrid determine if sub-style sets single_enable = 0 if (pairflag && force->pair == nullptr) - error->all(FLERR,"No pair style defined for compute group/group"); - if (force->pair_match("^hybrid",0) == nullptr - && force->pair->single_enable == 0) - error->all(FLERR,"Pair style does not support compute group/group"); + error->all(FLERR, "No pair style defined for compute group/group"); + if (force->pair_match("^hybrid", 0) == nullptr && force->pair->single_enable == 0) + error->all(FLERR, "Pair style does not support compute group/group"); // error if Kspace style does not compute group/group interactions if (kspaceflag && force->kspace == nullptr) - error->all(FLERR,"No Kspace style defined for compute group/group"); + error->all(FLERR, "No Kspace style defined for compute group/group"); if (kspaceflag && force->kspace->group_group_enable == 0) - error->all(FLERR,"Kspace style does not support compute group/group"); + error->all(FLERR, "Kspace style does not support compute group/group"); if (pairflag) { pair = force->pair; cutsq = force->pair->cutsq; - } else pair = nullptr; + } else + pair = nullptr; - if (kspaceflag) kspace = force->kspace; - else kspace = nullptr; + if (kspaceflag) + kspace = force->kspace; + else + kspace = nullptr; // compute Kspace correction terms if (kspaceflag) { kspace_correction(); if ((fabs(e_correction) > SMALL) && (comm->me == 0)) - error->warning(FLERR,"Both groups in compute group/group have a net charge; " + error->warning(FLERR, + "Both groups in compute group/group have a net charge; " "the Kspace boundary correction to energy will be non-zero"); } // recheck that group 2 has not been deleted jgroup = group->find(group2); - if (jgroup == -1) - error->all(FLERR,"Compute group/group group ID does not exist"); + if (jgroup == -1) error->all(FLERR, "Compute group/group group ID does not exist"); jgroupbit = group->bitmask[jgroup]; // need an occasional half neighbor list - if (pairflag) { - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->pair = 0; - neighbor->requests[irequest]->compute = 1; - neighbor->requests[irequest]->occasional = 1; - } + if (pairflag) neighbor->add_request(this, NeighConst::REQ_OCCASIONAL); } /* ---------------------------------------------------------------------- */ @@ -197,10 +195,10 @@ void ComputeGroupGroup::compute_vector() void ComputeGroupGroup::pair_contribution() { - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz; - double rsq,eng,fpair,factor_coul,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz; + double rsq, eng, fpair, factor_coul, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; tagint *molecule = atom->molecule; @@ -269,11 +267,11 @@ void ComputeGroupGroup::pair_contribution() 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; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); + eng = pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); // energy only computed once so tally full amount // force tally is jgroup acting on igroup @@ -281,25 +279,25 @@ void ComputeGroupGroup::pair_contribution() if (newton_pair || j < nlocal) { one[0] += eng; if (ij_flag) { - one[1] += delx*fpair; - one[2] += dely*fpair; - one[3] += delz*fpair; + one[1] += delx * fpair; + one[2] += dely * fpair; + one[3] += delz * fpair; } if (ji_flag) { - one[1] -= delx*fpair; - one[2] -= dely*fpair; - one[3] -= delz*fpair; + one[1] -= delx * fpair; + one[2] -= dely * fpair; + one[3] -= delz * fpair; } - // energy computed twice so tally half amount - // only tally force if I own igroup atom + // energy computed twice so tally half amount + // only tally force if I own igroup atom } else { - one[0] += 0.5*eng; + one[0] += 0.5 * eng; if (ij_flag) { - one[1] += delx*fpair; - one[2] += dely*fpair; - one[3] += delz*fpair; + one[1] += delx * fpair; + one[2] += dely * fpair; + one[3] += delz * fpair; } } } @@ -307,9 +305,11 @@ void ComputeGroupGroup::pair_contribution() } double all[4]; - MPI_Allreduce(one,all,4,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(one, all, 4, MPI_DOUBLE, MPI_SUM, world); scalar += all[0]; - vector[0] += all[1]; vector[1] += all[2]; vector[2] += all[3]; + vector[0] += all[1]; + vector[1] += all[2]; + vector[2] += all[3]; } /* ---------------------------------------------------------------------- */ @@ -318,8 +318,8 @@ void ComputeGroupGroup::kspace_contribution() { double *vector_kspace = force->kspace->f2group; - force->kspace->compute_group_group(groupbit,jgroupbit,0); - scalar += 2.0*force->kspace->e2group; + force->kspace->compute_group_group(groupbit, jgroupbit, 0); + scalar += 2.0 * force->kspace->e2group; vector[0] += vector_kspace[0]; vector[1] += vector_kspace[1]; vector[2] += vector_kspace[2]; @@ -328,7 +328,7 @@ void ComputeGroupGroup::kspace_contribution() // real-space style of compute group-group // add extra Kspace term to energy - force->kspace->compute_group_group(groupbit,jgroupbit,1); + force->kspace->compute_group_group(groupbit, jgroupbit, 1); scalar -= force->kspace->e2group; // self energy correction term @@ -345,8 +345,8 @@ void ComputeGroupGroup::kspace_contribution() // adjustment of z dimension for 2d slab Ewald // 3d Ewald just uses zprd since slab_volfactor = 1.0 - double volume = xprd*yprd*zprd*force->kspace->slab_volfactor; - scalar -= e_correction/volume; + double volume = xprd * yprd * zprd * force->kspace->slab_volfactor; + scalar -= e_correction / volume; } } @@ -357,7 +357,7 @@ void ComputeGroupGroup::kspace_correction() // total charge of groups A & B, needed for correction term - double qsqsum_group,qsum_A,qsum_B; + double qsqsum_group, qsum_A, qsum_B; qsqsum_group = qsum_A = qsum_B = 0.0; double *q = atom->q; @@ -366,20 +366,19 @@ void ComputeGroupGroup::kspace_correction() int groupbit_B = jgroupbit; for (int i = 0; i < atom->nlocal; i++) { - if ((mask[i] & groupbit_A) && (mask[i] & groupbit_B)) - qsqsum_group += q[i]*q[i]; + if ((mask[i] & groupbit_A) && (mask[i] & groupbit_B)) qsqsum_group += q[i] * q[i]; if (mask[i] & groupbit_A) qsum_A += q[i]; if (mask[i] & groupbit_B) qsum_B += q[i]; } double tmp; - MPI_Allreduce(&qsqsum_group,&tmp,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&qsqsum_group, &tmp, 1, MPI_DOUBLE, MPI_SUM, world); qsqsum_group = tmp; - MPI_Allreduce(&qsum_A,&tmp,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&qsum_A, &tmp, 1, MPI_DOUBLE, MPI_SUM, world); qsum_A = tmp; - MPI_Allreduce(&qsum_B,&tmp,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&qsum_B, &tmp, 1, MPI_DOUBLE, MPI_SUM, world); qsum_B = tmp; double g_ewald = force->kspace->g_ewald; @@ -389,29 +388,28 @@ void ComputeGroupGroup::kspace_correction() // self-energy correction - e_self = qscale * g_ewald*qsqsum_group/MY_PIS; - e_correction = 2.0*qsum_A*qsum_B; + e_self = qscale * g_ewald * qsqsum_group / MY_PIS; + e_correction = 2.0 * qsum_A * qsum_B; // subtract extra AA terms qsum_A = qsum_B = 0.0; for (int i = 0; i < atom->nlocal; i++) { - if (!((mask[i] & groupbit_A) && (mask[i] & groupbit_B))) - continue; + if (!((mask[i] & groupbit_A) && (mask[i] & groupbit_B))) continue; if (mask[i] & groupbit_A) qsum_A += q[i]; if (mask[i] & groupbit_B) qsum_B += q[i]; } - MPI_Allreduce(&qsum_A,&tmp,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&qsum_A, &tmp, 1, MPI_DOUBLE, MPI_SUM, world); qsum_A = tmp; - MPI_Allreduce(&qsum_B,&tmp,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&qsum_B, &tmp, 1, MPI_DOUBLE, MPI_SUM, world); qsum_B = tmp; // k=0 energy correction term (still need to divide by volume above) - e_correction -= qsum_A*qsum_B; - e_correction *= qscale * MY_PI2 / (g_ewald*g_ewald); + e_correction -= qsum_A * qsum_B; + e_correction *= qscale * MY_PI2 / (g_ewald * g_ewald); } diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp index fd020e1865..163f68e5ac 100644 --- a/src/compute_orientorder_atom.cpp +++ b/src/compute_orientorder_atom.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -28,23 +27,21 @@ #include "memory.h" #include "modify.h" #include "neigh_list.h" -#include "neigh_request.h" #include "neighbor.h" #include "pair.h" #include "update.h" -#include #include +#include using namespace LAMMPS_NS; using namespace MathConst; using namespace MathSpecial; - #ifdef DBL_EPSILON - #define MY_EPSILON (10.0*DBL_EPSILON) +#define MY_EPSILON (10.0 * DBL_EPSILON) #else - #define MY_EPSILON (10.0*2.220446049250313e-16) +#define MY_EPSILON (10.0 * 2.220446049250313e-16) #endif #define QEPSILON 1.0e-6 @@ -52,11 +49,10 @@ using namespace MathSpecial; /* ---------------------------------------------------------------------- */ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - qlist(nullptr), distsq(nullptr), nearest(nullptr), rlist(nullptr), - qnarray(nullptr), qnm_r(nullptr), qnm_i(nullptr), cglist(nullptr) + Compute(lmp, narg, arg), qlist(nullptr), distsq(nullptr), nearest(nullptr), rlist(nullptr), + qnarray(nullptr), qnm_r(nullptr), qnm_i(nullptr), cglist(nullptr) { - if (narg < 3 ) error->all(FLERR,"Illegal compute orientorder/atom command"); + if (narg < 3) error->all(FLERR, "Illegal compute orientorder/atom command"); // set default values for optional args @@ -70,7 +66,7 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg // specify which orders to request nqlist = 5; - memory->create(qlist,nqlist,"orientorder/atom:qlist"); + memory->create(qlist, nqlist, "orientorder/atom:qlist"); qlist[0] = 4; qlist[1] = 6; qlist[2] = 8; @@ -82,69 +78,69 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg int iarg = 3; while (iarg < narg) { - if (strcmp(arg[iarg],"nnn") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); - if (strcmp(arg[iarg+1],"NULL") == 0) { + if (strcmp(arg[iarg], "nnn") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute orientorder/atom command"); + if (strcmp(arg[iarg + 1], "NULL") == 0) { nnn = 0; } else { - nnn = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (nnn <= 0) error->all(FLERR,"Illegal compute orientorder/atom command"); + nnn = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (nnn <= 0) error->all(FLERR, "Illegal compute orientorder/atom command"); } iarg += 2; - } else if (strcmp(arg[iarg],"degrees") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); - nqlist = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (nqlist <= 0) error->all(FLERR,"Illegal compute orientorder/atom command"); + } else if (strcmp(arg[iarg], "degrees") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute orientorder/atom command"); + nqlist = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (nqlist <= 0) error->all(FLERR, "Illegal compute orientorder/atom command"); memory->destroy(qlist); - memory->create(qlist,nqlist,"orientorder/atom:qlist"); + memory->create(qlist, nqlist, "orientorder/atom:qlist"); iarg += 2; - if (iarg+nqlist > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); + if (iarg + nqlist > narg) error->all(FLERR, "Illegal compute orientorder/atom command"); qmax = 0; for (int il = 0; il < nqlist; il++) { - qlist[il] = utils::numeric(FLERR,arg[iarg+il],false,lmp); - if (qlist[il] < 0) error->all(FLERR,"Illegal compute orientorder/atom command"); + qlist[il] = utils::numeric(FLERR, arg[iarg + il], false, lmp); + if (qlist[il] < 0) error->all(FLERR, "Illegal compute orientorder/atom command"); if (qlist[il] > qmax) qmax = qlist[il]; } iarg += nqlist; - } else if (strcmp(arg[iarg],"wl") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute orientorder/atom command"); - wlflag = utils::logical(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "wl") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute orientorder/atom command"); + wlflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"wl/hat") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); - wlhatflag = utils::logical(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "wl/hat") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute orientorder/atom command"); + wlhatflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"components") == 0) { + } else if (strcmp(arg[iarg], "components") == 0) { qlcompflag = 1; - if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); - qlcomp = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute orientorder/atom command"); + qlcomp = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iqlcomp = -1; for (int il = 0; il < nqlist; il++) if (qlcomp == qlist[il]) { iqlcomp = il; break; } - if (iqlcomp == -1) error->all(FLERR,"Illegal compute orientorder/atom command"); + if (iqlcomp == -1) error->all(FLERR, "Illegal compute orientorder/atom command"); iarg += 2; - } else if (strcmp(arg[iarg],"cutoff") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); - double cutoff = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (cutoff <= 0.0) error->all(FLERR,"Illegal compute orientorder/atom command"); - cutsq = cutoff*cutoff; + } else if (strcmp(arg[iarg], "cutoff") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute orientorder/atom command"); + double cutoff = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (cutoff <= 0.0) error->all(FLERR, "Illegal compute orientorder/atom command"); + cutsq = cutoff * cutoff; iarg += 2; - } else if (strcmp(arg[iarg],"chunksize") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); - chunksize = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (chunksize <= 0) error->all(FLERR,"Illegal compute orientorder/atom command"); + } else if (strcmp(arg[iarg], "chunksize") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute orientorder/atom command"); + chunksize = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (chunksize <= 0) error->all(FLERR, "Illegal compute orientorder/atom command"); iarg += 2; - } else error->all(FLERR,"Illegal compute orientorder/atom command"); + } else + error->all(FLERR, "Illegal compute orientorder/atom command"); } ncol = nqlist; if (wlflag) ncol += nqlist; if (wlhatflag) ncol += nqlist; - if (qlcompflag) ncol += 2*(2*qlcomp+1); + if (qlcompflag) ncol += 2 * (2 * qlcomp + 1); peratom_flag = 1; size_peratom_cols = ncol; @@ -174,32 +170,23 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom() void ComputeOrientOrderAtom::init() { if (force->pair == nullptr) - error->all(FLERR,"Compute orientorder/atom requires a " - "pair style be defined"); - if (cutsq == 0.0) cutsq = force->pair->cutforce * force->pair->cutforce; + error->all(FLERR, "Compute orientorder/atom requires a pair style be defined"); + if (cutsq == 0.0) + cutsq = force->pair->cutforce * force->pair->cutforce; else if (sqrt(cutsq) > force->pair->cutforce) - error->all(FLERR,"Compute orientorder/atom cutoff is " - "longer than pairwise cutoff"); + error->all(FLERR, "Compute orientorder/atom cutoff is longer than pairwise cutoff"); memory->destroy(qnm_r); memory->destroy(qnm_i); - memory->create(qnm_r,nqlist,2*qmax+1,"orientorder/atom:qnm_r"); - memory->create(qnm_i,nqlist,2*qmax+1,"orientorder/atom:qnm_i"); + memory->create(qnm_r, nqlist, 2 * qmax + 1, "orientorder/atom:qnm_r"); + memory->create(qnm_i, nqlist, 2 * qmax + 1, "orientorder/atom:qnm_i"); // 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); - int count = 0; - for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style,"orientorder/atom") == 0) count++; - if (count > 1 && comm->me == 0) - error->warning(FLERR,"More than one compute orientorder/atom"); + if ((modify->get_compute_by_style("orientorder/atom").size() > 1) && (comm->me == 0)) + error->warning(FLERR, "More than one instance of compute orientorder/atom"); if (wlflag || wlhatflag) init_clebsch_gordan(); } @@ -215,9 +202,9 @@ void ComputeOrientOrderAtom::init_list(int /*id*/, NeighList *ptr) void ComputeOrientOrderAtom::compute_peratom() { - int i,j,ii,jj,inum,jnum; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum; + double xtmp, ytmp, ztmp, delx, dely, delz, rsq; + int *ilist, *jlist, *numneigh, **firstneigh; invoked_peratom = update->ntimestep; @@ -226,7 +213,7 @@ void ComputeOrientOrderAtom::compute_peratom() if (atom->nmax > nmax) { memory->destroy(qnarray); nmax = atom->nmax; - memory->create(qnarray,nmax,ncol,"orientorder/atom:qnarray"); + memory->create(qnarray, nmax, ncol, "orientorder/atom:qnarray"); array_atom = qnarray; } @@ -244,11 +231,11 @@ void ComputeOrientOrderAtom::compute_peratom() double **x = atom->x; int *mask = atom->mask; - memset(&qnarray[0][0],0,sizeof(double)*nmax*ncol); + memset(&qnarray[0][0], 0, sizeof(double) * nmax * ncol); for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - double* qn = qnarray[i]; + double *qn = qnarray[i]; if (mask[i] & groupbit) { xtmp = x[i][0]; ytmp = x[i][1]; @@ -263,9 +250,9 @@ void ComputeOrientOrderAtom::compute_peratom() memory->destroy(rlist); memory->destroy(nearest); maxneigh = jnum; - memory->create(distsq,maxneigh,"orientorder/atom:distsq"); - memory->create(rlist,maxneigh,3,"orientorder/atom:rlist"); - memory->create(nearest,maxneigh,"orientorder/atom:nearest"); + memory->create(distsq, maxneigh, "orientorder/atom:distsq"); + memory->create(rlist, maxneigh, 3, "orientorder/atom:rlist"); + memory->create(nearest, maxneigh, "orientorder/atom:nearest"); } // loop over list of all neighbors within force cutoff @@ -281,7 +268,7 @@ void ComputeOrientOrderAtom::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[ncount] = rsq; rlist[ncount][0] = delx; @@ -294,15 +281,14 @@ void ComputeOrientOrderAtom::compute_peratom() // if not nnn neighbors, order parameter = 0; if ((ncount == 0) || (ncount < nnn)) { - for (jj = 0; jj < ncol; jj++) - qn[jj] = 0.0; + for (jj = 0; jj < ncol; jj++) qn[jj] = 0.0; continue; } // if nnn > 0, use only nearest nnn neighbors if (nnn > 0) { - select3(nnn,ncount,distsq,nearest,rlist); + select3(nnn, ncount, distsq, nearest, rlist); ncount = nnn; } @@ -317,9 +303,9 @@ void ComputeOrientOrderAtom::compute_peratom() double ComputeOrientOrderAtom::memory_usage() { - double bytes = (double)ncol*nmax * sizeof(double); - bytes += (double)(qmax*(2*qmax+1)+maxneigh*4) * sizeof(double); - bytes += (double)(nqlist+maxneigh) * sizeof(int); + double bytes = (double) ncol * nmax * sizeof(double); + bytes += (double) (qmax * (2 * qmax + 1) + maxneigh * 4) * sizeof(double); + bytes += (double) (nqlist + maxneigh) * sizeof(int); return bytes; } @@ -331,26 +317,39 @@ double ComputeOrientOrderAtom::memory_usage() // Use no-op do while to create single statement -#define SWAP(a,b) do { \ - tmp = a; a = b; b = tmp; \ +#define SWAP(a, b) \ + do { \ + tmp = a; \ + a = b; \ + b = tmp; \ } while (0) -#define ISWAP(a,b) do { \ - itmp = a; a = b; b = itmp; \ +#define ISWAP(a, b) \ + do { \ + itmp = a; \ + a = b; \ + b = itmp; \ } while (0) -#define SWAP3(a,b) do { \ - tmp = a[0]; a[0] = b[0]; b[0] = tmp; \ - tmp = a[1]; a[1] = b[1]; b[1] = tmp; \ - tmp = a[2]; a[2] = b[2]; b[2] = tmp; \ +#define SWAP3(a, b) \ + do { \ + tmp = a[0]; \ + a[0] = b[0]; \ + b[0] = tmp; \ + tmp = a[1]; \ + a[1] = b[1]; \ + b[1] = tmp; \ + tmp = a[2]; \ + a[2] = b[2]; \ + b[2] = tmp; \ } while (0) /* ---------------------------------------------------------------------- */ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, double **arr3) { - int i,ir,j,l,mid,ia,itmp; - double a,tmp,a3[3]; + int i, ir, j, l, mid, ia, itmp; + double a, tmp, a3[3]; arr--; iarr--; @@ -358,59 +357,61 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl 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]); - SWAP3(arr3[l],arr3[ir]); + if (ir <= l + 1) { + if (ir == l + 1 && arr[ir] < arr[l]) { + SWAP(arr[l], arr[ir]); + ISWAP(iarr[l], iarr[ir]); + SWAP3(arr3[l], arr3[ir]); } return; } else { - mid=(l+ir) >> 1; - SWAP(arr[mid],arr[l+1]); - ISWAP(iarr[mid],iarr[l+1]); - SWAP3(arr3[mid],arr3[l+1]); + mid = (l + ir) >> 1; + SWAP(arr[mid], arr[l + 1]); + ISWAP(iarr[mid], iarr[l + 1]); + SWAP3(arr3[mid], arr3[l + 1]); if (arr[l] > arr[ir]) { - SWAP(arr[l],arr[ir]); - ISWAP(iarr[l],iarr[ir]); - SWAP3(arr3[l],arr3[ir]); + SWAP(arr[l], arr[ir]); + ISWAP(iarr[l], iarr[ir]); + SWAP3(arr3[l], arr3[ir]); } - if (arr[l+1] > arr[ir]) { - SWAP(arr[l+1],arr[ir]); - ISWAP(iarr[l+1],iarr[ir]); - SWAP3(arr3[l+1],arr3[ir]); + if (arr[l + 1] > arr[ir]) { + SWAP(arr[l + 1], arr[ir]); + ISWAP(iarr[l + 1], iarr[ir]); + SWAP3(arr3[l + 1], arr3[ir]); } - if (arr[l] > arr[l+1]) { - SWAP(arr[l],arr[l+1]); - ISWAP(iarr[l],iarr[l+1]); - SWAP3(arr3[l],arr3[l+1]); + if (arr[l] > arr[l + 1]) { + SWAP(arr[l], arr[l + 1]); + ISWAP(iarr[l], iarr[l + 1]); + SWAP3(arr3[l], arr3[l + 1]); } - i = l+1; + i = l + 1; j = ir; - a = arr[l+1]; - ia = iarr[l+1]; - a3[0] = arr3[l+1][0]; - a3[1] = arr3[l+1][1]; - a3[2] = arr3[l+1][2]; + a = arr[l + 1]; + ia = iarr[l + 1]; + a3[0] = arr3[l + 1][0]; + a3[1] = arr3[l + 1][1]; + a3[2] = arr3[l + 1][2]; 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]); - SWAP3(arr3[i],arr3[j]); + SWAP(arr[i], arr[j]); + ISWAP(iarr[i], iarr[j]); + SWAP3(arr3[i], arr3[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; - arr3[l+1][0] = arr3[j][0]; - arr3[l+1][1] = arr3[j][1]; - arr3[l+1][2] = arr3[j][2]; + arr3[l + 1][0] = arr3[j][0]; + arr3[l + 1][1] = arr3[j][1]; + arr3[l + 1][2] = arr3[j][2]; arr3[j][0] = a3[0]; arr3[j][1] = a3[1]; arr3[j][2] = a3[2]; - if (j >= k) ir = j-1; + if (j >= k) ir = j - 1; if (j <= k) l = i; } } @@ -420,34 +421,32 @@ void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, doubl calculate the bond orientational order parameters ------------------------------------------------------------------------- */ -void ComputeOrientOrderAtom::calc_boop(double **rlist, - int ncount, double qn[], - int qlist[], int nqlist) { +void ComputeOrientOrderAtom::calc_boop(double **rlist, int ncount, double qn[], int qlist[], + int nqlist) +{ for (int il = 0; il < nqlist; il++) { int l = qlist[il]; - for (int m = 0; m < 2*l+1; m++) { + for (int m = 0; m < 2 * l + 1; m++) { qnm_r[il][m] = 0.0; qnm_i[il][m] = 0.0; } } for (int ineigh = 0; ineigh < ncount; ineigh++) { - const double * const r = rlist[ineigh]; - double rmag = dist(r); - if (rmag <= MY_EPSILON) { - return; - } + const double *const r = rlist[ineigh]; + double rmag = sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]); + if (rmag <= MY_EPSILON) { return; } double costheta = r[2] / rmag; double expphi_r = r[0]; double expphi_i = r[1]; - double rxymag = sqrt(expphi_r*expphi_r+expphi_i*expphi_i); + double rxymag = sqrt(expphi_r * expphi_r + expphi_i * expphi_i); if (rxymag <= MY_EPSILON) { expphi_r = 1.0; expphi_i = 0.0; } else { - double rxymaginv = 1.0/rxymag; + double rxymaginv = 1.0 / rxymag; expphi_r *= rxymaginv; expphi_i *= rxymaginv; } @@ -467,21 +466,20 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, double prefactor = polar_prefactor(l, m, costheta); double ylm_r = prefactor * expphim_r; double ylm_i = prefactor * expphim_i; - qnm_r[il][m+l] += ylm_r; - qnm_i[il][m+l] += ylm_i; + qnm_r[il][m + l] += ylm_r; + qnm_i[il][m + l] += ylm_i; if (m & 1) { - qnm_r[il][-m+l] -= ylm_r; - qnm_i[il][-m+l] += ylm_i; + qnm_r[il][-m + l] -= ylm_r; + qnm_i[il][-m + l] += ylm_i; } else { - qnm_r[il][-m+l] += ylm_r; - qnm_i[il][-m+l] -= ylm_i; + qnm_r[il][-m + l] += ylm_r; + qnm_i[il][-m + l] -= ylm_i; } - double tmp_r = expphim_r*expphi_r - expphim_i*expphi_i; - double tmp_i = expphim_r*expphi_i + expphim_i*expphi_r; + double tmp_r = expphim_r * expphi_r - expphim_i * expphi_i; + double tmp_i = expphim_r * expphi_i + expphim_i * expphi_r; expphim_r = tmp_r; expphim_i = tmp_i; } - } } @@ -490,7 +488,7 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, double facn = 1.0 / ncount; for (int il = 0; il < nqlist; il++) { int l = qlist[il]; - for (int m = 0; m < 2*l+1; m++) { + for (int m = 0; m < 2 * l + 1; m++) { qnm_r[il][m] *= facn; qnm_i[il][m] *= facn; } @@ -502,10 +500,10 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, int jj = 0; for (int il = 0; il < nqlist; il++) { int l = qlist[il]; - double qnormfac = sqrt(MY_4PI/(2*l+1)); + double qnormfac = sqrt(MY_4PI / (2 * l + 1)); double qm_sum = 0.0; - for (int m = 0; m < 2*l+1; m++) - qm_sum += qnm_r[il][m]*qnm_r[il][m] + qnm_i[il][m]*qnm_i[il][m]; + for (int m = 0; m < 2 * l + 1; m++) + qm_sum += qnm_r[il][m] * qnm_r[il][m] + qnm_i[il][m] * qnm_i[il][m]; qn[jj++] = qnormfac * sqrt(qm_sum); } @@ -516,16 +514,16 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, for (int il = 0; il < nqlist; il++) { int l = qlist[il]; double wlsum = 0.0; - for (int m1 = 0; m1 < 2*l+1; m1++) { - for (int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) { + for (int m1 = 0; m1 < 2 * l + 1; m1++) { + for (int m2 = MAX(0, l - m1); m2 < MIN(2 * l + 1, 3 * l - m1 + 1); m2++) { int m = m1 + m2 - l; - double qm1qm2_r = qnm_r[il][m1]*qnm_r[il][m2] - qnm_i[il][m1]*qnm_i[il][m2]; - double qm1qm2_i = qnm_r[il][m1]*qnm_i[il][m2] + qnm_i[il][m1]*qnm_r[il][m2]; - wlsum += (qm1qm2_r*qnm_r[il][m] + qm1qm2_i*qnm_i[il][m])*cglist[idxcg_count]; + double qm1qm2_r = qnm_r[il][m1] * qnm_r[il][m2] - qnm_i[il][m1] * qnm_i[il][m2]; + double qm1qm2_i = qnm_r[il][m1] * qnm_i[il][m2] + qnm_i[il][m1] * qnm_r[il][m2]; + wlsum += (qm1qm2_r * qnm_r[il][m] + qm1qm2_i * qnm_i[il][m]) * cglist[idxcg_count]; idxcg_count++; } } - qn[jj++] = wlsum/sqrt(2*l+1); + qn[jj++] = wlsum / sqrt(2 * l + 1); } } @@ -536,21 +534,21 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, for (int il = 0; il < nqlist; il++) { int l = qlist[il]; double wlsum = 0.0; - for (int m1 = 0; m1 < 2*l+1; m1++) { - for (int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) { + for (int m1 = 0; m1 < 2 * l + 1; m1++) { + for (int m2 = MAX(0, l - m1); m2 < MIN(2 * l + 1, 3 * l - m1 + 1); m2++) { int m = m1 + m2 - l; - double qm1qm2_r = qnm_r[il][m1]*qnm_r[il][m2] - qnm_i[il][m1]*qnm_i[il][m2]; - double qm1qm2_i = qnm_r[il][m1]*qnm_i[il][m2] + qnm_i[il][m1]*qnm_r[il][m2]; - wlsum += (qm1qm2_r*qnm_r[il][m] + qm1qm2_i*qnm_i[il][m])*cglist[idxcg_count]; + double qm1qm2_r = qnm_r[il][m1] * qnm_r[il][m2] - qnm_i[il][m1] * qnm_i[il][m2]; + double qm1qm2_i = qnm_r[il][m1] * qnm_i[il][m2] + qnm_i[il][m1] * qnm_r[il][m2]; + wlsum += (qm1qm2_r * qnm_r[il][m] + qm1qm2_i * qnm_i[il][m]) * cglist[idxcg_count]; idxcg_count++; } } if (qn[il] < QEPSILON) qn[jj++] = 0.0; else { - double qnormfac = sqrt(MY_4PI/(2*l+1)); - double qnfac = qnormfac/qn[il]; - qn[jj++] = wlsum/sqrt(2*l+1)*(qnfac*qnfac*qnfac); + double qnormfac = sqrt(MY_4PI / (2 * l + 1)); + double qnfac = qnormfac / qn[il]; + qn[jj++] = wlsum / sqrt(2 * l + 1) * (qnfac * qnfac * qnfac); } } } @@ -561,29 +559,19 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, int il = iqlcomp; int l = qlcomp; if (qn[il] < QEPSILON) - for (int m = 0; m < 2*l+1; m++) { + for (int m = 0; m < 2 * l + 1; m++) { qn[jj++] = 0.0; qn[jj++] = 0.0; } else { - double qnormfac = sqrt(MY_4PI/(2*l+1)); - double qnfac = qnormfac/qn[il]; - for (int m = 0; m < 2*l+1; m++) { + double qnormfac = sqrt(MY_4PI / (2 * l + 1)); + double qnfac = qnormfac / qn[il]; + for (int m = 0; m < 2 * l + 1; m++) { qn[jj++] = qnm_r[il][m] * qnfac; qn[jj++] = qnm_i[il][m] * qnfac; } } } - -} - -/* ---------------------------------------------------------------------- - calculate scalar distance -------------------------------------------------------------------------- */ - -double ComputeOrientOrderAtom::dist(const double r[]) -{ - return sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]); } /* ---------------------------------------------------------------------- @@ -596,11 +584,10 @@ double ComputeOrientOrderAtom::polar_prefactor(int l, int m, double costheta) const int mabs = abs(m); double prefactor = 1.0; - for (int i=l-mabs+1; i < l+mabs+1; ++i) - prefactor *= static_cast(i); + for (int i = l - mabs + 1; i < l + mabs + 1; ++i) prefactor *= static_cast(i); - prefactor = sqrt(static_cast(2*l+1)/(MY_4PI*prefactor)) - * associated_legendre(l,mabs,costheta); + prefactor = sqrt(static_cast(2 * l + 1) / (MY_4PI * prefactor)) * + associated_legendre(l, mabs, costheta); if ((m < 0) && (m % 2)) prefactor = -prefactor; @@ -619,16 +606,15 @@ double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x) double p(1.0), pm1(0.0), pm2(0.0); if (m != 0) { - const double msqx = -sqrt(1.0-x*x); - for (int i=1; i < m+1; ++i) - p *= static_cast(2*i-1) * msqx; + const double msqx = -sqrt(1.0 - x * x); + for (int i = 1; i < m + 1; ++i) p *= static_cast(2 * i - 1) * msqx; } - for (int i=m+1; i < l+1; ++i) { + for (int i = m + 1; i < l + 1; ++i) { pm2 = pm1; pm1 = p; - p = (static_cast(2*i-1)*x*pm1 - - static_cast(i+m-1)*pm2) / static_cast(i-m); + p = (static_cast(2 * i - 1) * x * pm1 - static_cast(i + m - 1) * pm2) / + static_cast(i - m); } return p; @@ -642,16 +628,15 @@ double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x) void ComputeOrientOrderAtom::init_clebsch_gordan() { - double sum,dcg,sfaccg, sfac1, sfac2; + double sum, dcg, sfaccg, sfac1, sfac2; int m, aa2, bb2, cc2; int ifac, idxcg_count; idxcg_count = 0; for (int il = 0; il < nqlist; il++) { int l = qlist[il]; - for (int m1 = 0; m1 < 2*l+1; m1++) - for (int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) - idxcg_count++; + for (int m1 = 0; m1 < 2 * l + 1; m1++) + for (int m2 = MAX(0, l - m1); m2 < MIN(2 * l + 1, 3 * l - m1 + 1); m2++) idxcg_count++; } idxcg_max = idxcg_count; memory->destroy(cglist); @@ -660,41 +645,44 @@ void ComputeOrientOrderAtom::init_clebsch_gordan() idxcg_count = 0; for (int il = 0; il < nqlist; il++) { int l = qlist[il]; - for (int m1 = 0; m1 < 2*l+1; m1++) { - aa2 = m1 - l; - for (int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) { - bb2 = m2 - l; - m = aa2 + bb2 + l; + for (int m1 = 0; m1 < 2 * l + 1; m1++) { + aa2 = m1 - l; + for (int m2 = MAX(0, l - m1); m2 < MIN(2 * l + 1, 3 * l - m1 + 1); m2++) { + bb2 = m2 - l; + m = aa2 + bb2 + l; - sum = 0.0; - for (int z = MAX(0, MAX(-aa2, bb2)); - z <= MIN(l, MIN(l - aa2, l + bb2)); z++) { - ifac = z % 2 ? -1 : 1; - sum += ifac / - (factorial(z) * - factorial(l - z) * - factorial(l - aa2 - z) * - factorial(l + bb2 - z) * - factorial(aa2 + z) * - factorial(-bb2 + z)); - } + // clang-format off - cc2 = m - l; - sfaccg = sqrt(factorial(l + aa2) * - factorial(l - aa2) * - factorial(l + bb2) * - factorial(l - bb2) * - factorial(l + cc2) * - factorial(l - cc2) * - (2*l + 1)); - - sfac1 = factorial(3*l + 1); - sfac2 = factorial(l); - dcg = sqrt(sfac2*sfac2*sfac2 / sfac1); - - cglist[idxcg_count] = sum * dcg * sfaccg; - idxcg_count++; + sum = 0.0; + for (int z = MAX(0, MAX(-aa2, bb2)); + z <= MIN(l, MIN(l - aa2, l + bb2)); z++) { + ifac = z % 2 ? -1 : 1; + sum += ifac / + (factorial(z) * + factorial(l - z) * + factorial(l - aa2 - z) * + factorial(l + bb2 - z) * + factorial(aa2 + z) * + factorial(-bb2 + z)); } + + cc2 = m - l; + sfaccg = sqrt(factorial(l + aa2) * + factorial(l - aa2) * + factorial(l + bb2) * + factorial(l - bb2) * + factorial(l + cc2) * + factorial(l - cc2) * + (2*l + 1)); + // clang-format on + + sfac1 = factorial(3 * l + 1); + sfac2 = factorial(l); + dcg = sqrt(sfac2 * sfac2 * sfac2 / sfac1); + + cglist[idxcg_count] = sum * dcg * sfaccg; + idxcg_count++; } + } } } diff --git a/src/compute_orientorder_atom.h b/src/compute_orientorder_atom.h index b88b01fbc9..c762c0e74b 100644 --- a/src/compute_orientorder_atom.h +++ b/src/compute_orientorder_atom.h @@ -50,7 +50,6 @@ class ComputeOrientOrderAtom : public Compute { void select3(int, int, double *, int *, double **); void calc_boop(double **rlist, int numNeighbors, double qn[], int nlist[], int nnlist); - double dist(const double r[]); double polar_prefactor(int, int, double); double associated_legendre(int, int, double); diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index 3dc37707c6..708b31c370 100644 --- a/src/compute_pair_local.cpp +++ b/src/compute_pair_local.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -31,16 +30,15 @@ using namespace LAMMPS_NS; #define DELTA 10000 -enum{DIST,ENG,FORCE,FX,FY,FZ,PN,DX,DY,DZ}; -enum{TYPE,RADIUS}; +enum { DIST, ENG, FORCE, FX, FY, FZ, PN, DX, DY, DZ }; +enum { TYPE, RADIUS }; /* ---------------------------------------------------------------------- */ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - pstyle(nullptr), pindex(nullptr), vlocal(nullptr), alocal(nullptr) + Compute(lmp, narg, arg), pstyle(nullptr), pindex(nullptr), vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute pair/local command"); + if (narg < 4) error->all(FLERR, "Illegal compute pair/local command"); local_flag = 1; nvalues = narg - 3; @@ -50,23 +48,32 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : nvalues = 0; int iarg = 3; while (iarg < narg) { - if (strcmp(arg[iarg],"dist") == 0) pstyle[nvalues++] = DIST; - else if (strcmp(arg[iarg],"eng") == 0) pstyle[nvalues++] = ENG; - else if (strcmp(arg[iarg],"force") == 0) pstyle[nvalues++] = FORCE; - else if (strcmp(arg[iarg],"fx") == 0) pstyle[nvalues++] = FX; - else if (strcmp(arg[iarg],"fy") == 0) pstyle[nvalues++] = FY; - else if (strcmp(arg[iarg],"fz") == 0) pstyle[nvalues++] = FZ; - else if (strcmp(arg[iarg],"dx") == 0) pstyle[nvalues++] = DX; - else if (strcmp(arg[iarg],"dy") == 0) pstyle[nvalues++] = DY; - else if (strcmp(arg[iarg],"dz") == 0) pstyle[nvalues++] = DZ; + if (strcmp(arg[iarg], "dist") == 0) + pstyle[nvalues++] = DIST; + else if (strcmp(arg[iarg], "eng") == 0) + pstyle[nvalues++] = ENG; + else if (strcmp(arg[iarg], "force") == 0) + pstyle[nvalues++] = FORCE; + else if (strcmp(arg[iarg], "fx") == 0) + pstyle[nvalues++] = FX; + else if (strcmp(arg[iarg], "fy") == 0) + pstyle[nvalues++] = FY; + else if (strcmp(arg[iarg], "fz") == 0) + pstyle[nvalues++] = FZ; + else if (strcmp(arg[iarg], "dx") == 0) + pstyle[nvalues++] = DX; + else if (strcmp(arg[iarg], "dy") == 0) + pstyle[nvalues++] = DY; + else if (strcmp(arg[iarg], "dz") == 0) + pstyle[nvalues++] = DZ; else if (arg[iarg][0] == 'p') { int n = atoi(&arg[iarg][1]); - if (n <= 0) error->all(FLERR, - "Invalid keyword in compute pair/local command"); + if (n <= 0) error->all(FLERR, "Invalid keyword in compute pair/local command"); pstyle[nvalues] = PN; - pindex[nvalues++] = n-1; + pindex[nvalues++] = n - 1; - } else break; + } else + break; iarg++; } @@ -76,20 +83,23 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : cutstyle = TYPE; while (iarg < narg) { - if (strcmp(arg[iarg],"cutoff") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute pair/local command"); - if (strcmp(arg[iarg+1],"type") == 0) cutstyle = TYPE; - else if (strcmp(arg[iarg+1],"radius") == 0) cutstyle = RADIUS; - else error->all(FLERR,"Illegal compute pair/local command"); + if (strcmp(arg[iarg], "cutoff") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute pair/local command"); + if (strcmp(arg[iarg + 1], "type") == 0) + cutstyle = TYPE; + else if (strcmp(arg[iarg + 1], "radius") == 0) + cutstyle = RADIUS; + else + error->all(FLERR, "Illegal compute pair/local command"); iarg += 2; - } else error->all(FLERR,"Illegal compute pair/local command"); + } else + error->all(FLERR, "Illegal compute pair/local command"); } // error check if (cutstyle == RADIUS && !atom->radius_flag) - error->all(FLERR,"Compute pair/local requires atom attribute radius"); + error->all(FLERR, "Compute pair/local requires atom attribute radius"); // set singleflag if need to call pair->single() @@ -97,8 +107,10 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : for (int i = 0; i < nvalues; i++) if (pstyle[i] != DIST && pstyle[i] != DX && pstyle[i] != DY && pstyle[i] != DZ) singleflag = 1; - if (nvalues == 1) size_local_cols = 0; - else size_local_cols = nvalues; + if (nvalues == 1) + size_local_cols = 0; + else + size_local_cols = nvalues; nmax = 0; vlocal = nullptr; @@ -111,8 +123,8 @@ ComputePairLocal::~ComputePairLocal() { memory->destroy(vlocal); memory->destroy(alocal); - delete [] pstyle; - delete [] pindex; + delete[] pstyle; + delete[] pindex; } /* ---------------------------------------------------------------------- */ @@ -120,25 +132,22 @@ ComputePairLocal::~ComputePairLocal() void ComputePairLocal::init() { if (singleflag && force->pair == nullptr) - error->all(FLERR,"No pair style is defined for compute pair/local"); + error->all(FLERR, "No pair style is defined for compute pair/local"); if (singleflag && force->pair->single_enable == 0) - error->all(FLERR,"Pair style does not support compute pair/local"); + error->all(FLERR, "Pair style does not support compute pair/local"); for (int i = 0; i < nvalues; i++) if (pstyle[i] == PN && pindex[i] >= force->pair->single_extra) - error->all(FLERR,"Pair style does not have extra field" - " requested by compute pair/local"); + error->all(FLERR, "Pair style does not have extra field requested by compute pair/local"); // need an occasional half neighbor list // set size to same value as request made by force->pair // this should enable it to always be a copy list (e.g. for granular pstyle) - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->pair = 0; - neighbor->requests[irequest]->compute = 1; - neighbor->requests[irequest]->occasional = 1; - NeighRequest *pairrequest = neighbor->find_request((void *) force->pair); - if (pairrequest) neighbor->requests[irequest]->size = pairrequest->size; + int neighflags = NeighConst::REQ_OCCASIONAL; + auto pairrequest = neighbor->find_request(force->pair); + if (pairrequest && pairrequest->get_size()) neighflags |= NeighConst::REQ_SIZE; + neighbor->add_request(this, neighflags); } /* ---------------------------------------------------------------------- */ @@ -171,11 +180,11 @@ void ComputePairLocal::compute_local() int ComputePairLocal::compute_pairs(int flag) { - int i,j,m,n,ii,jj,inum,jnum,itype,jtype; - tagint itag,jtag; - double xtmp,ytmp,ztmp,delx,dely,delz; - double rsq,radsum,eng,fpair,factor_coul,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, m, n, ii, jj, inum, jnum, itype, jtype; + tagint itag, jtag; + double xtmp, ytmp, ztmp, delx, dely, delz; + double rsq, radsum, eng, fpair, factor_coul, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; double *ptr; double **x = atom->x; @@ -234,9 +243,9 @@ int ComputePairLocal::compute_pairs(int flag) if (newton_pair == 0 && j >= nlocal) { if (itag > jtag) { - if ((itag+jtag) % 2 == 0) continue; + if ((itag + jtag) % 2 == 0) continue; } else if (itag < jtag) { - if ((itag+jtag) % 2 == 1) continue; + if ((itag + jtag) % 2 == 1) continue; } else { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { @@ -249,59 +258,62 @@ int ComputePairLocal::compute_pairs(int flag) 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; jtype = type[j]; if (cutstyle == TYPE) { if (rsq >= cutsq[itype][jtype]) continue; } else { radsum = radius[i] + radius[j]; - if (rsq >= radsum*radsum) continue; + if (rsq >= radsum * radsum) continue; } if (flag) { if (singleflag) - eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); - else eng = fpair = 0.0; + eng = pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); + else + eng = fpair = 0.0; - if (nvalues == 1) ptr = &vlocal[m]; - else ptr = alocal[m]; + if (nvalues == 1) + ptr = &vlocal[m]; + else + ptr = alocal[m]; // to make sure dx, dy and dz are always from the lower to the higher id double directionCorrection = itag > jtag ? -1.0 : 1.0; for (n = 0; n < nvalues; n++) { switch (pstyle[n]) { - case DIST: - ptr[n] = sqrt(rsq); - break; - case DX: - ptr[n] = delx*directionCorrection; - break; - case DY: - ptr[n] = dely*directionCorrection; - break; - case DZ: - ptr[n] = delz*directionCorrection; - break; - case ENG: - ptr[n] = eng; - break; - case FORCE: - ptr[n] = sqrt(rsq)*fpair; - break; - case FX: - ptr[n] = delx*fpair; - break; - case FY: - ptr[n] = dely*fpair; - break; - case FZ: - ptr[n] = delz*fpair; - break; - case PN: - ptr[n] = pair->svector[pindex[n]]; - break; + case DIST: + ptr[n] = sqrt(rsq); + break; + case DX: + ptr[n] = delx * directionCorrection; + break; + case DY: + ptr[n] = dely * directionCorrection; + break; + case DZ: + ptr[n] = delz * directionCorrection; + break; + case ENG: + ptr[n] = eng; + break; + case FORCE: + ptr[n] = sqrt(rsq) * fpair; + break; + case FX: + ptr[n] = delx * fpair; + break; + case FY: + ptr[n] = dely * fpair; + break; + case FZ: + ptr[n] = delz * fpair; + break; + case PN: + ptr[n] = pair->svector[pindex[n]]; + break; } } } @@ -323,11 +335,11 @@ void ComputePairLocal::reallocate(int n) if (nvalues == 1) { memory->destroy(vlocal); - memory->create(vlocal,nmax,"pair/local:vector_local"); + memory->create(vlocal, nmax, "pair/local:vector_local"); vector_local = vlocal; } else { memory->destroy(alocal); - memory->create(alocal,nmax,nvalues,"pair/local:array_local"); + memory->create(alocal, nmax, nvalues, "pair/local:array_local"); array_local = alocal; } } @@ -338,6 +350,6 @@ void ComputePairLocal::reallocate(int n) double ComputePairLocal::memory_usage() { - double bytes = (double)nmax*nvalues * sizeof(double); + double bytes = (double) nmax * nvalues * sizeof(double); return bytes; } diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index fce18fecd1..b00535d5ed 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -264,12 +264,10 @@ void ComputePropertyLocal::init() // this should enable it to always be a copy list (e.g. for granular pstyle) if (kindflag == NEIGH || kindflag == PAIR) { - int irequest = neighbor->request(this, instance_me); - neighbor->requests[irequest]->pair = 0; - neighbor->requests[irequest]->compute = 1; - neighbor->requests[irequest]->occasional = 1; - NeighRequest *pairrequest = neighbor->find_request((void *) force->pair); - if (pairrequest) neighbor->requests[irequest]->size = pairrequest->size; + int neighflags = NeighConst::REQ_OCCASIONAL; + auto pairrequest = neighbor->find_request(force->pair); + if (pairrequest && pairrequest->get_size()) neighflags |= NeighConst::REQ_SIZE; + neighbor->add_request(this, neighflags); } // do initial memory allocation so that memory_usage() is correct diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp index 8be704c07c..947d1383a6 100644 --- a/src/neigh_request.cpp +++ b/src/neigh_request.cpp @@ -296,3 +296,10 @@ void NeighRequest::set_kokkos_host(int flag) { kokkos_host = flag; } + +void NeighRequest::set_skip(int *_iskip, int **_ijskip) +{ + skip = 1; + iskip = _iskip; + ijskip = _ijskip; +} diff --git a/src/neigh_request.h b/src/neigh_request.h index 1bf3caa9da..59453c88c8 100644 --- a/src/neigh_request.h +++ b/src/neigh_request.h @@ -126,8 +126,10 @@ class NeighRequest : protected Pointers { void set_kokkos_device(int); void set_kokkos_host(int); void set_cutoff(double); + void set_skip(int *, int **); int get_size() const { return size; } + void *get_requestor() const { return requestor; } }; } // namespace LAMMPS_NS diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 45dc5cfecd..0935b7a9e7 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -1742,6 +1742,18 @@ NeighRequest *Neighbor::find_request(void *classptr) return nullptr; } +/* ---------------------------------------------------------------------- + return vector with neighbor list requests from pair styles +------------------------------------------------------------------------- */ + +const std::vector Neighbor::get_pair_requests() const +{ + std::vector matches; + for (int i=0; i < nrequest; ++i) + if (requests[i]->pair) matches.push_back(requests[i]); + return matches; +} + /* ---------------------------------------------------------------------- find and return list requested by classptr if not found or classptr = nullptr, return nullptr diff --git a/src/neighbor.h b/src/neighbor.h index 3f3a75aa1e..b98e48dd66 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -152,6 +152,7 @@ class Neighbor : protected Pointers { int exclude_setting(); // return exclude value to accelerator pkg NeighList *find_list(void *); // find a neighbor request NeighRequest *find_request(void *); // find a neighbor request + const std::vector get_pair_requests() const; int any_full(); // Check if any old requests had full neighbor lists void build_collection(int); // build peratom collection array starting at the given index diff --git a/src/pair_buck_coul_cut.cpp b/src/pair_buck_coul_cut.cpp index 95cec11eff..065086ed2d 100644 --- a/src/pair_buck_coul_cut.cpp +++ b/src/pair_buck_coul_cut.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -68,13 +67,13 @@ PairBuckCoulCut::~PairBuckCoulCut() void PairBuckCoulCut::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; - double rsq,r2inv,r6inv,r,rexp,forcecoul,forcebuck,factor_coul,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair; + double rsq, r2inv, r6inv, r, rexp, forcecoul, forcebuck, factor_coul, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; evdwl = ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -112,47 +111,49 @@ void PairBuckCoulCut::compute(int eflag, int vflag) 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; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; r = sqrt(rsq); if (rsq < cut_coulsq[itype][jtype]) - forcecoul = qqrd2e * qtmp*q[j]/r; - else forcecoul = 0.0; + forcecoul = qqrd2e * qtmp * q[j] / r; + else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - r6inv = r2inv*r2inv*r2inv; - rexp = exp(-r*rhoinv[itype][jtype]); - forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv; - } else forcebuck = 0.0; + r6inv = r2inv * r2inv * r2inv; + rexp = exp(-r * rhoinv[itype][jtype]); + forcebuck = buck1[itype][jtype] * r * rexp - buck2[itype][jtype] * r6inv; + } else + forcebuck = 0.0; - fpair = (factor_coul*forcecoul + factor_lj*forcebuck) * r2inv; + fpair = (factor_coul * forcecoul + factor_lj * forcebuck) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { if (rsq < cut_coulsq[itype][jtype]) - ecoul = factor_coul * qqrd2e * qtmp*q[j]/r; - else ecoul = 0.0; + ecoul = factor_coul * qqrd2e * qtmp * q[j] / r; + else + ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv - - offset[itype][jtype]; + evdwl = a[itype][jtype] * rexp - c[itype][jtype] * r6inv - offset[itype][jtype]; evdwl *= factor_lj; - } else evdwl = 0.0; + } else + evdwl = 0.0; } - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz); } } } @@ -167,27 +168,26 @@ void PairBuckCoulCut::compute(int eflag, int vflag) void PairBuckCoulCut::allocate() { allocated = 1; - int n = atom->ntypes; + int np1 = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); + memory->create(setflag, np1, np1, "pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + for (int i = 1; i < np1; i++) + for (int j = i; j < np1; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq, np1, np1, "pair:cutsq"); - memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); - memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); - memory->create(cut_coul,n+1,n+1,"pair:cut_coul"); - memory->create(cut_coulsq,n+1,n+1,"pair:cut_coulsq"); - memory->create(a,n+1,n+1,"pair:a"); - memory->create(rho,n+1,n+1,"pair:rho"); - memory->create(c,n+1,n+1,"pair:c"); - memory->create(rhoinv,n+1,n+1,"pair:rhoinv"); - memory->create(buck1,n+1,n+1,"pair:buck1"); - memory->create(buck2,n+1,n+1,"pair:buck2"); - memory->create(offset,n+1,n+1,"pair:offset"); + memory->create(cut_lj, np1, np1, "pair:cut_lj"); + memory->create(cut_ljsq, np1, np1, "pair:cut_ljsq"); + memory->create(cut_coul, np1, np1, "pair:cut_coul"); + memory->create(cut_coulsq, np1, np1, "pair:cut_coulsq"); + memory->create(a, np1, np1, "pair:a"); + memory->create(rho, np1, np1, "pair:rho"); + memory->create(c, np1, np1, "pair:c"); + memory->create(rhoinv, np1, np1, "pair:rhoinv"); + memory->create(buck1, np1, np1, "pair:buck1"); + memory->create(buck2, np1, np1, "pair:buck2"); + memory->create(offset, np1, np1, "pair:offset"); } /* ---------------------------------------------------------------------- @@ -196,16 +196,18 @@ void PairBuckCoulCut::allocate() void PairBuckCoulCut::settings(int narg, char **arg) { - if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command"); + if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command"); - cut_lj_global = utils::numeric(FLERR,arg[0],false,lmp); - if (narg == 1) cut_coul_global = cut_lj_global; - else cut_coul_global = utils::numeric(FLERR,arg[1],false,lmp); + cut_lj_global = utils::numeric(FLERR, arg[0], false, lmp); + if (narg == 1) + cut_coul_global = cut_lj_global; + else + cut_coul_global = utils::numeric(FLERR, arg[1], false, lmp); // reset cutoffs that have been explicitly set if (allocated) { - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) if (setflag[i][j]) { @@ -221,27 +223,26 @@ void PairBuckCoulCut::settings(int narg, char **arg) void PairBuckCoulCut::coeff(int narg, char **arg) { - if (narg < 5 || narg > 7) - error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 5 || narg > 7) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); - double a_one = utils::numeric(FLERR,arg[2],false,lmp); - double rho_one = utils::numeric(FLERR,arg[3],false,lmp); - if (rho_one <= 0) error->all(FLERR,"Incorrect args for pair coefficients"); - double c_one = utils::numeric(FLERR,arg[4],false,lmp); + double a_one = utils::numeric(FLERR, arg[2], false, lmp); + double rho_one = utils::numeric(FLERR, arg[3], false, lmp); + if (rho_one <= 0) error->all(FLERR, "Incorrect args for pair coefficients"); + double c_one = utils::numeric(FLERR, arg[4], false, lmp); double cut_lj_one = cut_lj_global; double cut_coul_one = cut_coul_global; - if (narg >= 6) cut_coul_one = cut_lj_one = utils::numeric(FLERR,arg[5],false,lmp); - if (narg == 7) cut_coul_one = utils::numeric(FLERR,arg[6],false,lmp); + if (narg >= 6) cut_coul_one = cut_lj_one = utils::numeric(FLERR, arg[5], false, lmp); + if (narg == 7) cut_coul_one = utils::numeric(FLERR, arg[6], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { a[i][j] = a_one; rho[i][j] = rho_one; c[i][j] = c_one; @@ -252,7 +253,7 @@ void PairBuckCoulCut::coeff(int narg, char **arg) } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -261,10 +262,9 @@ void PairBuckCoulCut::coeff(int narg, char **arg) void PairBuckCoulCut::init_style() { - if (!atom->q_flag) - error->all(FLERR,"Pair style buck/coul/cut requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair style buck/coul/cut requires atom attribute q"); - neighbor->request(this,instance_me); + neighbor->add_request(this, NeighConst::REQ_DEFAULT); } /* ---------------------------------------------------------------------- @@ -273,20 +273,21 @@ void PairBuckCoulCut::init_style() double PairBuckCoulCut::init_one(int i, int j) { - if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set"); - double cut = MAX(cut_lj[i][j],cut_coul[i][j]); + double cut = MAX(cut_lj[i][j], cut_coul[i][j]); cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j]; - rhoinv[i][j] = 1.0/rho[i][j]; - buck1[i][j] = a[i][j]/rho[i][j]; - buck2[i][j] = 6.0*c[i][j]; + rhoinv[i][j] = 1.0 / rho[i][j]; + buck1[i][j] = a[i][j] / rho[i][j]; + buck2[i][j] = 6.0 * c[i][j]; if (offset_flag && (cut_lj[i][j] > 0.0)) { - double rexp = exp(-cut_lj[i][j]/rho[i][j]); - offset[i][j] = a[i][j]*rexp - c[i][j]/pow(cut_lj[i][j],6.0); - } else offset[i][j] = 0.0; + double rexp = exp(-cut_lj[i][j] / rho[i][j]); + offset[i][j] = a[i][j] * rexp - c[i][j] / pow(cut_lj[i][j], 6.0); + } else + offset[i][j] = 0.0; cut_ljsq[j][i] = cut_ljsq[i][j]; cut_coulsq[j][i] = cut_coulsq[i][j]; @@ -304,26 +305,26 @@ double PairBuckCoulCut::init_one(int i, int j) int *type = atom->type; int nlocal = atom->nlocal; - double count[2],all[2]; + double count[2], all[2]; count[0] = count[1] = 0.0; for (int k = 0; k < nlocal; k++) { if (type[k] == i) count[0] += 1.0; if (type[k] == j) count[1] += 1.0; } - MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world); double rho1 = rho[i][j]; - double rho2 = rho1*rho1; - double rho3 = rho2*rho1; + double rho2 = rho1 * rho1; + double rho3 = rho2 * rho1; double rc = cut_lj[i][j]; - double rc2 = rc*rc; - double rc3 = rc2*rc; - etail_ij = 2.0*MY_PI*all[0]*all[1]* - (a[i][j]*exp(-rc/rho1)*rho1*(rc2 + 2.0*rho1*rc + 2.0*rho2) - - c[i][j]/(3.0*rc3)); - ptail_ij = (-1/3.0)*2.0*MY_PI*all[0]*all[1]* - (-a[i][j]*exp(-rc/rho1)* - (rc3 + 3.0*rho1*rc2 + 6.0*rho2*rc + 6.0*rho3) + 2.0*c[i][j]/rc3); + double rc2 = rc * rc; + double rc3 = rc2 * rc; + etail_ij = 2.0 * MY_PI * all[0] * all[1] * + (a[i][j] * exp(-rc / rho1) * rho1 * (rc2 + 2.0 * rho1 * rc + 2.0 * rho2) - + c[i][j] / (3.0 * rc3)); + ptail_ij = (-1 / 3.0) * 2.0 * MY_PI * all[0] * all[1] * + (-a[i][j] * exp(-rc / rho1) * (rc3 + 3.0 * rho1 * rc2 + 6.0 * rho2 * rc + 6.0 * rho3) + + 2.0 * c[i][j] / rc3); } return cut; @@ -337,16 +338,16 @@ void PairBuckCoulCut::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); + fwrite(&setflag[i][j], sizeof(int), 1, fp); if (setflag[i][j]) { - fwrite(&a[i][j],sizeof(double),1,fp); - fwrite(&rho[i][j],sizeof(double),1,fp); - fwrite(&c[i][j],sizeof(double),1,fp); - fwrite(&cut_lj[i][j],sizeof(double),1,fp); - fwrite(&cut_coul[i][j],sizeof(double),1,fp); + fwrite(&a[i][j], sizeof(double), 1, fp); + fwrite(&rho[i][j], sizeof(double), 1, fp); + fwrite(&c[i][j], sizeof(double), 1, fp); + fwrite(&cut_lj[i][j], sizeof(double), 1, fp); + fwrite(&cut_coul[i][j], sizeof(double), 1, fp); } } } @@ -361,25 +362,25 @@ void PairBuckCoulCut::read_restart(FILE *fp) allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); if (setflag[i][j]) { if (me == 0) { - utils::sfread(FLERR,&a[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&rho[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&c[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_lj[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul[i][j],sizeof(double),1,fp,nullptr,error); + utils::sfread(FLERR, &a[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &rho[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &c[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_lj[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul[i][j], sizeof(double), 1, fp, nullptr, error); } - MPI_Bcast(&a[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&rho[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&c[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&a[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&rho[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&c[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_lj[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -390,11 +391,11 @@ void PairBuckCoulCut::read_restart(FILE *fp) void PairBuckCoulCut::write_restart_settings(FILE *fp) { - fwrite(&cut_lj_global,sizeof(double),1,fp); - fwrite(&cut_coul_global,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&cut_lj_global, sizeof(double), 1, fp); + fwrite(&cut_coul_global, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); + fwrite(&tail_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -404,17 +405,17 @@ void PairBuckCoulCut::write_restart_settings(FILE *fp) void PairBuckCoulCut::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&cut_lj_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_lj_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); - MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&cut_lj_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- @@ -424,7 +425,7 @@ void PairBuckCoulCut::read_restart_settings(FILE *fp) void PairBuckCoulCut::write_data(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d %g %g %g\n",i,a[i][i],rho[i][i],c[i][i]); + fprintf(fp, "%d %g %g %g\n", i, a[i][i], rho[i][i], c[i][i]); } /* ---------------------------------------------------------------------- @@ -435,40 +436,39 @@ void PairBuckCoulCut::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g %g %g %g %g\n",i,j, - a[i][j],rho[i][j],c[i][j],cut_lj[i][j],cut_coul[i][j]); + fprintf(fp, "%d %d %g %g %g %g %g\n", i, j, a[i][j], rho[i][j], c[i][j], cut_lj[i][j], + cut_coul[i][j]); } /* ---------------------------------------------------------------------- */ -double PairBuckCoulCut::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) +double PairBuckCoulCut::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, + double factor_lj, double &fforce) { - double r2inv,r6inv,r,rexp,forcecoul,forcebuck,phicoul,phibuck; + double r2inv, r6inv, r, rexp, forcecoul, forcebuck, phicoul, phibuck; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; if (rsq < cut_coulsq[itype][jtype]) - forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv); - else forcecoul = 0.0; + forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv); + else + forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { - r6inv = r2inv*r2inv*r2inv; + r6inv = r2inv * r2inv * r2inv; r = sqrt(rsq); - rexp = exp(-r*rhoinv[itype][jtype]); - forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv; - } else forcebuck = 0.0; - fforce = (factor_coul*forcecoul + factor_lj*forcebuck) * r2inv; + rexp = exp(-r * rhoinv[itype][jtype]); + forcebuck = buck1[itype][jtype] * r * rexp - buck2[itype][jtype] * r6inv; + } else + forcebuck = 0.0; + fforce = (factor_coul * forcecoul + factor_lj * forcebuck) * r2inv; double eng = 0.0; if (rsq < cut_coulsq[itype][jtype]) { - phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv); - eng += factor_coul*phicoul; + phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv); + eng += factor_coul * phicoul; } if (rsq < cut_ljsq[itype][jtype]) { - phibuck = a[itype][jtype]*rexp - c[itype][jtype]*r6inv - - offset[itype][jtype]; - eng += factor_lj*phibuck; + phibuck = a[itype][jtype] * rexp - c[itype][jtype] * r6inv - offset[itype][jtype]; + eng += factor_lj * phibuck; } return eng; } @@ -478,7 +478,7 @@ double PairBuckCoulCut::single(int i, int j, int itype, int jtype, void *PairBuckCoulCut::extract(const char *str, int &dim) { dim = 2; - if (strcmp(str,"a") == 0) return (void *) a; - if (strcmp(str,"c") == 0) return (void *) c; + if (strcmp(str, "a") == 0) return (void *) a; + if (strcmp(str, "c") == 0) return (void *) c; return nullptr; } diff --git a/src/pair_coul_cut.cpp b/src/pair_coul_cut.cpp index ea29358174..628de8185e 100644 --- a/src/pair_coul_cut.cpp +++ b/src/pair_coul_cut.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -29,7 +28,8 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairCoulCut::PairCoulCut(LAMMPS *lmp) : Pair(lmp) { +PairCoulCut::PairCoulCut(LAMMPS *lmp) : Pair(lmp) +{ writedata = 1; } @@ -52,13 +52,13 @@ PairCoulCut::~PairCoulCut() void PairCoulCut::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; - double rsq,r2inv,rinv,forcecoul,factor_coul; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, ecoul, fpair; + double rsq, r2inv, rinv, forcecoul, factor_coul; + int *ilist, *jlist, *numneigh, **firstneigh; ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -94,29 +94,27 @@ void PairCoulCut::compute(int eflag, int vflag) 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; jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; rinv = sqrt(r2inv); - forcecoul = qqrd2e * scale[itype][jtype] * qtmp*q[j]*rinv; - fpair = factor_coul*forcecoul * r2inv; + forcecoul = qqrd2e * scale[itype][jtype] * qtmp * q[j] * rinv; + fpair = factor_coul * forcecoul * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } - if (eflag) - ecoul = factor_coul * qqrd2e * scale[itype][jtype] * qtmp*q[j]*rinv; + if (eflag) ecoul = factor_coul * qqrd2e * scale[itype][jtype] * qtmp * q[j] * rinv; - if (evflag) ev_tally(i,j,nlocal,newton_pair, - 0.0,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, ecoul, fpair, delx, dely, delz); } } } @@ -131,19 +129,19 @@ void PairCoulCut::compute(int eflag, int vflag) void PairCoulCut::allocate() { allocated = 1; - int n = atom->ntypes; + int np1 = atom->ntypes + 1; - memory->create(setflag,n+1,n+1,"pair:setflag"); - memory->create(scale,n+1,n+1,"pair:scale"); - for (int i = 1; i <= n; i++) { - for (int j = i; j <= n; j++) { + memory->create(setflag, np1, np1, "pair:setflag"); + memory->create(scale, np1, np1, "pair:scale"); + for (int i = 1; i < np1; i++) { + for (int j = i; j < np1; j++) { setflag[i][j] = 0; scale[i][j] = 1.0; } } - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(cutsq, np1, np1, "pair:cutsq"); + memory->create(cut, np1, np1, "pair:cut"); } /* ---------------------------------------------------------------------- @@ -152,14 +150,14 @@ void PairCoulCut::allocate() void PairCoulCut::settings(int narg, char **arg) { - if (narg != 1) error->all(FLERR,"Illegal pair_style command"); + if (narg != 1) error->all(FLERR, "Illegal pair_style command"); - cut_global = utils::numeric(FLERR,arg[0],false,lmp); + cut_global = utils::numeric(FLERR, arg[0], false, lmp); // reset cutoffs that have been explicitly set if (allocated) { - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) if (setflag[i][j]) cut[i][j] = cut_global; @@ -172,20 +170,19 @@ void PairCoulCut::settings(int narg, char **arg) void PairCoulCut::coeff(int narg, char **arg) { - if (narg < 2 || narg > 3) - error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 2 || narg > 3) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); double cut_one = cut_global; - if (narg == 3) cut_one = utils::numeric(FLERR,arg[2],false,lmp); + if (narg == 3) cut_one = utils::numeric(FLERR, arg[2], false, lmp); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { cut[i][j] = cut_one; scale[i][j] = 1.0; setflag[i][j] = 1; @@ -193,20 +190,18 @@ void PairCoulCut::coeff(int narg, char **arg) } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } - /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairCoulCut::init_style() { - if (!atom->q_flag) - error->all(FLERR,"Pair style coul/cut requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair style coul/cut requires atom attribute q"); - neighbor->request(this,instance_me); + neighbor->add_request(this, NeighConst::REQ_DEFAULT); } /* ---------------------------------------------------------------------- @@ -216,7 +211,7 @@ void PairCoulCut::init_style() double PairCoulCut::init_one(int i, int j) { if (setflag[i][j] == 0) { - cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); scale[i][j] = 1.0; } @@ -233,14 +228,12 @@ void PairCoulCut::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) { for (j = i; j <= atom->ntypes; j++) { - fwrite(&scale[i][j],sizeof(double),1,fp); - fwrite(&setflag[i][j],sizeof(int),1,fp); - if (setflag[i][j]) { - fwrite(&cut[i][j],sizeof(double),1,fp); - } + fwrite(&scale[i][j], sizeof(double), 1, fp); + fwrite(&setflag[i][j], sizeof(int), 1, fp); + if (setflag[i][j]) { fwrite(&cut[i][j], sizeof(double), 1, fp); } } } } @@ -254,20 +247,19 @@ void PairCoulCut::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) { for (j = i; j <= atom->ntypes; j++) { if (me == 0) { - utils::sfread(FLERR,&scale[i][j],sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &scale[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&scale[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + MPI_Bcast(&scale[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); if (setflag[i][j]) { - if (me == 0) - utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + if (me == 0) utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); } } } @@ -279,9 +271,9 @@ void PairCoulCut::read_restart(FILE *fp) void PairCoulCut::write_restart_settings(FILE *fp) { - fwrite(&cut_global,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); + fwrite(&cut_global, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -291,13 +283,13 @@ void PairCoulCut::write_restart_settings(FILE *fp) void PairCoulCut::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- @@ -306,8 +298,7 @@ void PairCoulCut::read_restart_settings(FILE *fp) void PairCoulCut::write_data(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp,"%d\n",i); + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d\n", i); } /* ---------------------------------------------------------------------- @@ -317,25 +308,23 @@ void PairCoulCut::write_data(FILE *fp) void PairCoulCut::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) - for (int j = i; j <= atom->ntypes; j++) - fprintf(fp,"%d %d %g\n",i,j,cut[i][j]); + for (int j = i; j <= atom->ntypes; j++) fprintf(fp, "%d %d %g\n", i, j, cut[i][j]); } /* ---------------------------------------------------------------------- */ -double PairCoulCut::single(int i, int j, int /*itype*/, int /*jtype*/, - double rsq, double factor_coul, double /*factor_lj*/, - double &fforce) +double PairCoulCut::single(int i, int j, int /*itype*/, int /*jtype*/, double rsq, + double factor_coul, double /*factor_lj*/, double &fforce) { - double r2inv,rinv,forcecoul,phicoul; + double r2inv, rinv, forcecoul, phicoul; - r2inv = 1.0/rsq; + r2inv = 1.0 / rsq; rinv = sqrt(r2inv); - forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*rinv; - fforce = factor_coul*forcecoul * r2inv; + forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * rinv; + fforce = factor_coul * forcecoul * r2inv; - phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*rinv; - return factor_coul*phicoul; + phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * rinv; + return factor_coul * phicoul; } /* ---------------------------------------------------------------------- */ @@ -343,7 +332,7 @@ double PairCoulCut::single(int i, int j, int /*itype*/, int /*jtype*/, void *PairCoulCut::extract(const char *str, int &dim) { dim = 2; - if (strcmp(str,"cut_coul") == 0) return (void *) cut; - if (strcmp(str,"scale") == 0) return (void *) scale; + if (strcmp(str, "cut_coul") == 0) return (void *) cut; + if (strcmp(str, "scale") == 0) return (void *) scale; return nullptr; } diff --git a/src/pair_coul_dsf.cpp b/src/pair_coul_dsf.cpp index 98c9f7f5d6..c84a257107 100644 --- a/src/pair_coul_dsf.cpp +++ b/src/pair_coul_dsf.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -19,27 +18,28 @@ #include "pair_coul_dsf.h" -#include -#include #include "atom.h" #include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "memory.h" -#include "math_const.h" #include "error.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; -#define EWALD_F 1.12837917 -#define EWALD_P 0.3275911 -#define A1 0.254829592 -#define A2 -0.284496736 -#define A3 1.421413741 -#define A4 -1.453152027 -#define A5 1.061405429 +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 /* ---------------------------------------------------------------------- */ @@ -61,14 +61,14 @@ PairCoulDSF::~PairCoulDSF() void PairCoulDSF::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; - double r,rsq,forcecoul,factor_coul; - double prefactor,erfcc,erfcd,t; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, ecoul, fpair; + double r, rsq, forcecoul, factor_coul; + double prefactor, erfcc, erfcd, t; + int *ilist, *jlist, *numneigh, **firstneigh; ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -95,8 +95,8 @@ void PairCoulDSF::compute(int eflag, int vflag) jnum = numneigh[i]; if (eflag) { - double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e; - ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0); + double e_self = -(e_shift / 2.0 + alpha / MY_PIS) * qtmp * qtmp * qqrd2e; + ev_tally(i, i, nlocal, 0, 0.0, e_self, 0.0, 0.0, 0.0, 0.0); } for (jj = 0; jj < jnum; jj++) { @@ -107,36 +107,35 @@ void PairCoulDSF::compute(int eflag, int vflag) 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 < cut_coulsq) { r = sqrt(rsq); - prefactor = qqrd2e*qtmp*q[j]/r; - erfcd = exp(-alpha*alpha*rsq); - t = 1.0 / (1.0 + EWALD_P*alpha*r); - erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; + prefactor = qqrd2e * qtmp * q[j] / r; + erfcd = exp(-alpha * alpha * rsq); + t = 1.0 / (1.0 + EWALD_P * alpha * r); + erfcc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * erfcd; - forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd + - r*f_shift) * r; - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + forcecoul = prefactor * (erfcc / r + 2.0 * alpha / MY_PIS * erfcd + r * f_shift) * r; + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; fpair = forcecoul / rsq; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { - ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); - if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; - } else ecoul = 0.0; + ecoul = prefactor * (erfcc - r * e_shift - rsq * f_shift); + if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; + } else + ecoul = 0.0; - if (evflag) ev_tally(i,j,nlocal,newton_pair, - 0.0,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, ecoul, fpair, delx, dely, delz); } } } @@ -153,12 +152,11 @@ void PairCoulDSF::allocate() allocated = 1; int n = atom->ntypes; - memory->create(setflag,n+1,n+1,"pair:setflag"); + memory->create(setflag, n + 1, n + 1, "pair:setflag"); for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + for (int j = i; j <= n; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); } /* ---------------------------------------------------------------------- @@ -167,10 +165,10 @@ void PairCoulDSF::allocate() void PairCoulDSF::settings(int narg, char **arg) { - if (narg != 2) error->all(FLERR,"Illegal pair_style command"); + if (narg != 2) error->all(FLERR, "Illegal pair_style command"); - alpha = utils::numeric(FLERR,arg[0],false,lmp); - cut_coul = utils::numeric(FLERR,arg[1],false,lmp); + alpha = utils::numeric(FLERR, arg[0], false, lmp); + cut_coul = utils::numeric(FLERR, arg[1], false, lmp); } /* ---------------------------------------------------------------------- @@ -179,22 +177,22 @@ void PairCoulDSF::settings(int narg, char **arg) void PairCoulDSF::coeff(int narg, char **arg) { - if (narg != 2) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg != 2) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { setflag[i][j] = 1; count++; } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -203,16 +201,15 @@ void PairCoulDSF::coeff(int narg, char **arg) void PairCoulDSF::init_style() { - if (!atom->q_flag) - error->all(FLERR,"Pair style coul/dsf requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair style coul/dsf requires atom attribute q"); - neighbor->request(this,instance_me); + neighbor->add_request(this, NeighConst::REQ_DEFAULT); cut_coulsq = cut_coul * cut_coul; - double erfcc = erfc(alpha*cut_coul); - double erfcd = exp(-alpha*alpha*cut_coul*cut_coul); - f_shift = -(erfcc/cut_coulsq + 2.0/MY_PIS*alpha*erfcd/cut_coul); - e_shift = erfcc/cut_coul - f_shift*cut_coul; + double erfcc = erfc(alpha * cut_coul); + double erfcd = exp(-alpha * alpha * cut_coul * cut_coul); + f_shift = -(erfcc / cut_coulsq + 2.0 / MY_PIS * alpha * erfcd / cut_coul); + e_shift = erfcc / cut_coul - f_shift * cut_coul; } /* ---------------------------------------------------------------------- @@ -232,11 +229,9 @@ void PairCoulDSF::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); - } + for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j], sizeof(int), 1, fp); } } /* ---------------------------------------------------------------------- @@ -248,12 +243,12 @@ void PairCoulDSF::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); } } @@ -263,10 +258,10 @@ void PairCoulDSF::read_restart(FILE *fp) void PairCoulDSF::write_restart_settings(FILE *fp) { - fwrite(&alpha,sizeof(double),1,fp); - fwrite(&cut_coul,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); + fwrite(&alpha, sizeof(double), 1, fp); + fwrite(&cut_coul, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -276,40 +271,38 @@ void PairCoulDSF::write_restart_settings(FILE *fp) void PairCoulDSF::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&alpha,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &alpha, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&alpha,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&alpha, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- */ double PairCoulDSF::single(int i, int j, int /*itype*/, int /*jtype*/, double rsq, - double factor_coul, double /*factor_lj*/, - double &fforce) + double factor_coul, double /*factor_lj*/, double &fforce) { - double r,erfcc,erfcd,prefactor,t; - double forcecoul,phicoul; + double r, erfcc, erfcd, prefactor, t; + double forcecoul, phicoul; forcecoul = phicoul = 0.0; if (rsq < cut_coulsq) { r = sqrt(rsq); - prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; - erfcd = exp(-alpha*alpha*rsq); - t = 1.0 / (1.0 + EWALD_P*alpha*r); - erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd; + prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r; + erfcd = exp(-alpha * alpha * rsq); + t = 1.0 / (1.0 + EWALD_P * alpha * r); + erfcc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * erfcd; - forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS*erfcd + - r*f_shift) * r; - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + forcecoul = prefactor * (erfcc / r + 2.0 * alpha / MY_PIS * erfcd + r * f_shift) * r; + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; - phicoul = prefactor * (erfcc - r*e_shift - rsq*f_shift); - if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + phicoul = prefactor * (erfcc - r * e_shift - rsq * f_shift); + if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor; } fforce = forcecoul / rsq; @@ -321,7 +314,7 @@ double PairCoulDSF::single(int i, int j, int /*itype*/, int /*jtype*/, double rs void *PairCoulDSF::extract(const char *str, int &dim) { - if (strcmp(str,"cut_coul") == 0) { + if (strcmp(str, "cut_coul") == 0) { dim = 0; return (void *) &cut_coul; } diff --git a/src/pair_coul_wolf.cpp b/src/pair_coul_wolf.cpp index 3a260abf28..eb4947383f 100644 --- a/src/pair_coul_wolf.cpp +++ b/src/pair_coul_wolf.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -18,16 +17,16 @@ #include "pair_coul_wolf.h" -#include #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" -#include "neighbor.h" -#include "neigh_list.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "neigh_list.h" +#include "neighbor.h" +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -36,7 +35,7 @@ using namespace MathConst; PairCoulWolf::PairCoulWolf(LAMMPS *lmp) : Pair(lmp) { - single_enable = 0; // NOTE: single() method below is not yet correct + single_enable = 0; // NOTE: single() method below is not yet correct } /* ---------------------------------------------------------------------- */ @@ -55,16 +54,16 @@ PairCoulWolf::~PairCoulWolf() void PairCoulWolf::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; - double rsq,forcecoul,factor_coul; + int i, j, ii, jj, inum, jnum; + double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, ecoul, fpair; + double rsq, forcecoul, factor_coul; double prefactor; double r; - int *ilist,*jlist,*numneigh,**firstneigh; - double erfcc,erfcd,v_sh,dvdrr,e_self,e_shift,f_shift,qisq; + int *ilist, *jlist, *numneigh, **firstneigh; + double erfcc, erfcd, v_sh, dvdrr, e_self, e_shift, f_shift, qisq; ecoul = 0.0; - ev_init(eflag,vflag); + ev_init(eflag, vflag); double **x = atom->x; double **f = atom->f; @@ -77,9 +76,8 @@ void PairCoulWolf::compute(int eflag, int vflag) // self and shifted coulombic energy e_self = v_sh = 0.0; - e_shift = erfc(alf*cut_coul)/cut_coul; - f_shift = -(e_shift+ 2.0*alf/MY_PIS * exp(-alf*alf*cut_coul*cut_coul)) / - cut_coul; + e_shift = erfc(alf * cut_coul) / cut_coul; + f_shift = -(e_shift + 2.0 * alf / MY_PIS * exp(-alf * alf * cut_coul * cut_coul)) / cut_coul; inum = list->inum; ilist = list->ilist; @@ -97,9 +95,9 @@ void PairCoulWolf::compute(int eflag, int vflag) jlist = firstneigh[i]; jnum = numneigh[i]; - qisq = qtmp*qtmp; - e_self = -(e_shift/2.0 + alf/MY_PIS) * qisq*qqrd2e; - if (evflag) ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0); + qisq = qtmp * qtmp; + e_self = -(e_shift / 2.0 + alf / MY_PIS) * qisq * qqrd2e; + if (evflag) ev_tally(i, i, nlocal, 0, 0.0, e_self, 0.0, 0.0, 0.0, 0.0); for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; @@ -109,35 +107,35 @@ void PairCoulWolf::compute(int eflag, int vflag) 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 < cut_coulsq) { r = sqrt(rsq); - prefactor = qqrd2e*qtmp*q[j]/r; - erfcc = erfc(alf*r); - erfcd = exp(-alf*alf*r*r); - v_sh = (erfcc - e_shift*r) * prefactor; - dvdrr = (erfcc/rsq + 2.0*alf/MY_PIS * erfcd/r) + f_shift; - forcecoul = dvdrr*rsq*prefactor; - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + prefactor = qqrd2e * qtmp * q[j] / r; + erfcc = erfc(alf * r); + erfcd = exp(-alf * alf * r * r); + v_sh = (erfcc - e_shift * r) * prefactor; + dvdrr = (erfcc / rsq + 2.0 * alf / MY_PIS * erfcd / r) + f_shift; + forcecoul = dvdrr * rsq * prefactor; + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; fpair = forcecoul / rsq; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; } if (eflag) { ecoul = v_sh; - if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; - } else ecoul = 0.0; + if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor; + } else + ecoul = 0.0; - if (evflag) ev_tally(i,j,nlocal,newton_pair, - 0.0,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, ecoul, fpair, delx, dely, delz); } } } @@ -154,12 +152,11 @@ void PairCoulWolf::allocate() allocated = 1; int n = atom->ntypes; - memory->create(setflag,n+1,n+1,"pair:setflag"); + memory->create(setflag, n + 1, n + 1, "pair:setflag"); for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + for (int j = i; j <= n; j++) setflag[i][j] = 0; - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); } /* ---------------------------------------------------------------------- @@ -170,10 +167,10 @@ void PairCoulWolf::allocate() void PairCoulWolf::settings(int narg, char **arg) { - if (narg != 2) error->all(FLERR,"Illegal pair_style command"); + if (narg != 2) error->all(FLERR, "Illegal pair_style command"); - alf = utils::numeric(FLERR,arg[0],false,lmp); - cut_coul = utils::numeric(FLERR,arg[1],false,lmp); + alf = utils::numeric(FLERR, arg[0], false, lmp); + cut_coul = utils::numeric(FLERR, arg[1], false, lmp); } /* ---------------------------------------------------------------------- @@ -182,22 +179,22 @@ void PairCoulWolf::settings(int narg, char **arg) void PairCoulWolf::coeff(int narg, char **arg) { - if (narg != 2) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg != 2) error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo,ihi,jlo,jhi; - utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); - utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); + int ilo, ihi, jlo, jhi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); int count = 0; for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { + for (int j = MAX(jlo, i); j <= jhi; j++) { setflag[i][j] = 1; count++; } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- @@ -206,12 +203,11 @@ void PairCoulWolf::coeff(int narg, char **arg) void PairCoulWolf::init_style() { - if (!atom->q_flag) - error->all(FLERR,"Pair coul/wolf requires atom attribute q"); + if (!atom->q_flag) error->all(FLERR, "Pair coul/wolf requires atom attribute q"); - neighbor->request(this,instance_me); + neighbor->add_request(this, NeighConst::REQ_DEFAULT); - cut_coulsq = cut_coul*cut_coul; + cut_coulsq = cut_coul * cut_coul; } /* ---------------------------------------------------------------------- @@ -231,10 +227,9 @@ void PairCoulWolf::write_restart(FILE *fp) { write_restart_settings(fp); - int i,j; + int i, j; for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) - fwrite(&setflag[i][j],sizeof(int),1,fp); + for (j = i; j <= atom->ntypes; j++) fwrite(&setflag[i][j], sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -246,12 +241,12 @@ void PairCoulWolf::read_restart(FILE *fp) read_restart_settings(fp); allocate(); - int i,j; + int i, j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { - if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); } } @@ -261,10 +256,10 @@ void PairCoulWolf::read_restart(FILE *fp) void PairCoulWolf::write_restart_settings(FILE *fp) { - fwrite(&alf,sizeof(double),1,fp); - fwrite(&cut_coul,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); + fwrite(&alf, sizeof(double), 1, fp); + fwrite(&cut_coul, sizeof(double), 1, fp); + fwrite(&offset_flag, sizeof(int), 1, fp); + fwrite(&mix_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -274,15 +269,15 @@ void PairCoulWolf::write_restart_settings(FILE *fp) void PairCoulWolf::read_restart_settings(FILE *fp) { if (comm->me == 0) { - utils::sfread(FLERR,&alf,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&cut_coul,sizeof(double),1,fp,nullptr,error); - utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); - utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); + utils::sfread(FLERR, &alf, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut_coul, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); } - MPI_Bcast(&alf,1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&alf, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_coul, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- @@ -290,32 +285,31 @@ void PairCoulWolf::read_restart_settings(FILE *fp) ------------------------------------------------------------------------- */ double PairCoulWolf::single(int i, int j, int /*itype*/, int /*jtype*/, double rsq, - double factor_coul, double /*factor_lj*/, - double &fforce) + double factor_coul, double /*factor_lj*/, double &fforce) { - double r,prefactor; - double forcecoul,phicoul; - double e_shift,f_shift,dvdrr,erfcc,erfcd; + double r, prefactor; + double forcecoul, phicoul; + double e_shift, f_shift, dvdrr, erfcc, erfcd; - e_shift = erfc(alf*cut_coul) / cut_coul; - f_shift = -(e_shift+ 2.0*alf/MY_PIS * exp(-alf*alf*cut_coul*cut_coul)) / - cut_coul; + e_shift = erfc(alf * cut_coul) / cut_coul; + f_shift = -(e_shift + 2.0 * alf / MY_PIS * exp(-alf * alf * cut_coul * cut_coul)) / cut_coul; if (rsq < cut_coulsq) { r = sqrt(rsq); - prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; - erfcc = erfc(alf*r); - erfcd = exp(-alf*alf*r*r); - dvdrr = (erfcc/rsq + 2.0*alf/MY_PIS * erfcd/r) + f_shift; - forcecoul = dvdrr*rsq*prefactor; - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; - } else forcecoul = 0.0; + prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r; + erfcc = erfc(alf * r); + erfcd = exp(-alf * alf * r * r); + dvdrr = (erfcc / rsq + 2.0 * alf / MY_PIS * erfcd / r) + f_shift; + forcecoul = dvdrr * rsq * prefactor; + if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor; + } else + forcecoul = 0.0; fforce = forcecoul / rsq; double eng = 0.0; if (rsq < cut_coulsq) { - phicoul = prefactor * (erfcc-e_shift*r); - if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + phicoul = prefactor * (erfcc - e_shift * r); + if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor; eng += phicoul; } return eng; diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index 6287fb97db..5deb35b634 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -610,13 +610,12 @@ void PairHybrid::init_style() // create skip lists inside each pair neigh request // any kind of list can have its skip flag set in this loop - for (i = 0; i < neighbor->nrequest; i++) { - if (!neighbor->requests[i]->pair) continue; + for (auto &request : neighbor->get_pair_requests()) { // istyle = associated sub-style for the request for (istyle = 0; istyle < nstyles; istyle++) - if (styles[istyle] == neighbor->requests[i]->requestor) break; + if (styles[istyle] == request->get_requestor()) break; // allocate iskip and ijskip // initialize so as to skip all pair types @@ -663,9 +662,7 @@ void PairHybrid::init_style() if (ijskip[itype][jtype] == 1) skip = 1; if (skip) { - neighbor->requests[i]->skip = 1; - neighbor->requests[i]->iskip = iskip; - neighbor->requests[i]->ijskip = ijskip; + request->set_skip(iskip, ijskip); } else { delete[] iskip; memory->destroy(ijskip);