From 21079b3ac2d11af619518f27051d0c20a1a62123 Mon Sep 17 00:00:00 2001 From: "Dan S. Bolintineanu" Date: Fri, 2 Oct 2020 16:15:30 -0600 Subject: [PATCH 01/25] Added multi-stencil files from Ishan --- .../npair_half_size_multi_newtoff.cpp | 119 +++ .../npair_half_size_multi_newtoff.h | 43 + .../npair_half_size_multi_newton.cpp | 145 +++ src/from_ishan/npair_half_size_multi_newton.h | 43 + .../npair_half_size_multi_newton_tri.cpp | 126 +++ .../npair_half_size_multi_newton_tri.h | 43 + .../pair_gran_hooke_history_multi.cpp | 885 ++++++++++++++++++ .../pair_gran_hooke_history_multi.h | 109 +++ 8 files changed, 1513 insertions(+) create mode 100644 src/from_ishan/npair_half_size_multi_newtoff.cpp create mode 100644 src/from_ishan/npair_half_size_multi_newtoff.h create mode 100644 src/from_ishan/npair_half_size_multi_newton.cpp create mode 100644 src/from_ishan/npair_half_size_multi_newton.h create mode 100644 src/from_ishan/npair_half_size_multi_newton_tri.cpp create mode 100644 src/from_ishan/npair_half_size_multi_newton_tri.h create mode 100644 src/from_ishan/pair_gran_hooke_history_multi.cpp create mode 100644 src/from_ishan/pair_gran_hooke_history_multi.h diff --git a/src/from_ishan/npair_half_size_multi_newtoff.cpp b/src/from_ishan/npair_half_size_multi_newtoff.cpp new file mode 100644 index 0000000000..ec4ff805e2 --- /dev/null +++ b/src/from_ishan/npair_half_size_multi_newtoff.cpp @@ -0,0 +1,119 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include +#include "npair_half_size_multi_newtoff.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewtoff::NPairHalfSizeMultiNewtoff(LAMMPS *lmp) : NPair(lmp) {} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with partial Newton's 3rd law + each owned atom i checks own bin and other bins in stencil + multi-type stencil is itype dependent and is distance checked + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewtoff::build(NeighList *list) +{ + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *tag = atom->tag; + tagint *molecule = atom->molecule; + int nlocal = atom->nlocal; + if (includegroup) nlocal = atom->nfirst; + + int history = list->history; + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage *ipage = list->ipage; + + int mask_history = 3 << SBBITS; + + int inum = 0; + ipage->reset(); + + for (i = 0; i < nlocal; i++) { + n = 0; + neighptr = ipage->vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over all atoms in other bins in stencil including self + // only store pair if i < j + // skip if i,j neighbor cutoff is less than bin distance + // stores own/own pairs only once + // stores own/ghost pairs on both procs + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage->vgot(n); + if (ipage->status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + + list->inum = inum; +} diff --git a/src/from_ishan/npair_half_size_multi_newtoff.h b/src/from_ishan/npair_half_size_multi_newtoff.h new file mode 100644 index 0000000000..8abe5b5456 --- /dev/null +++ b/src/from_ishan/npair_half_size_multi_newtoff.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newtoff, + NPairHalfSizeMultiNewtoff, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTOFF | NP_ORTHO | NP_TRI) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewtoff : public NPair { + public: + NPairHalfSizeMultiNewtoff(class LAMMPS *); + ~NPairHalfSizeMultiNewtoff() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/from_ishan/npair_half_size_multi_newton.cpp b/src/from_ishan/npair_half_size_multi_newton.cpp new file mode 100644 index 0000000000..01aea66632 --- /dev/null +++ b/src/from_ishan/npair_half_size_multi_newton.cpp @@ -0,0 +1,145 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include +#include "npair_half_size_multi_newton.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewton::NPairHalfSizeMultiNewton(LAMMPS *lmp) : NPair(lmp) {} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with full Newton's 3rd law + each owned atom i checks its own bin and other bins in Newton stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewton::build(NeighList *list) +{ + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *tag = atom->tag; + tagint *molecule = atom->molecule; + int nlocal = atom->nlocal; + if (includegroup) nlocal = atom->nfirst; + + int history = list->history; + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage *ipage = list->ipage; + + int mask_history = 3 << SBBITS; + + int inum = 0; + ipage->reset(); + + for (i = 0; i < nlocal; i++) { + n = 0; + neighptr = ipage->vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over rest of atoms in i's bin, ghosts are at end of linked list + // if j is owned atom, store it, since j is beyond i in linked list + // if j is ghost, only store if j coords are "above and to the right" of i + + for (j = bins[i]; j >= 0; j = bins[j]) { + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + + // loop over all atoms in other bins in stencil, store every pair + // skip if i,j neighbor cutoff is less than bin distance + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage->vgot(n); + if (ipage->status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + + list->inum = inum; +} diff --git a/src/from_ishan/npair_half_size_multi_newton.h b/src/from_ishan/npair_half_size_multi_newton.h new file mode 100644 index 0000000000..109fe78c2b --- /dev/null +++ b/src/from_ishan/npair_half_size_multi_newton.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newton, + NPairHalfSizeMultiNewton, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_ORTHO) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewton : public NPair { + public: + NPairHalfSizeMultiNewton(class LAMMPS *); + ~NPairHalfSizeMultiNewton() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/from_ishan/npair_half_size_multi_newton_tri.cpp b/src/from_ishan/npair_half_size_multi_newton_tri.cpp new file mode 100644 index 0000000000..f7123a958d --- /dev/null +++ b/src/from_ishan/npair_half_size_multi_newton_tri.cpp @@ -0,0 +1,126 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include +#include "npair_half_size_multi_newton_tri.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewtonTri::NPairHalfSizeMultiNewtonTri(LAMMPS *lmp) : NPair(lmp) {} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with Newton's 3rd law for triclinic + each owned atom i checks its own bin and other bins in triclinic stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewtonTri::build(NeighList *list) +{ + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *tag = atom->tag; + tagint *molecule = atom->molecule; + int nlocal = atom->nlocal; + if (includegroup) nlocal = atom->nfirst; + + int history = list->history; + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage *ipage = list->ipage; + + int mask_history = 3 << SBBITS; + + int inum = 0; + ipage->reset(); + + for (i = 0; i < nlocal; i++) { + n = 0; + neighptr = ipage->vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + + // loop over all atoms in bins in stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and lower x) + // (equal zyx and j <= i) + // latter excludes self-self interaction but allows superposed atoms + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp) { + if (x[j][0] < xtmp) continue; + if (x[j][0] == xtmp && j <= i) continue; + } + } + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage->vgot(n); + if (ipage->status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + + list->inum = inum; +} diff --git a/src/from_ishan/npair_half_size_multi_newton_tri.h b/src/from_ishan/npair_half_size_multi_newton_tri.h new file mode 100644 index 0000000000..15fbaba0f7 --- /dev/null +++ b/src/from_ishan/npair_half_size_multi_newton_tri.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newton/tri, + NPairHalfSizeMultiNewtonTri, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_TRI) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewtonTri : public NPair { + public: + NPairHalfSizeMultiNewtonTri(class LAMMPS *); + ~NPairHalfSizeMultiNewtonTri() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/from_ishan/pair_gran_hooke_history_multi.cpp b/src/from_ishan/pair_gran_hooke_history_multi.cpp new file mode 100644 index 0000000000..83f5221741 --- /dev/null +++ b/src/from_ishan/pair_gran_hooke_history_multi.cpp @@ -0,0 +1,885 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Leo Silbert (SNL), Gary Grest (SNL) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_gran_hooke_history_multi.h" +#include "atom.h" +#include "atom_vec.h" +#include "domain.h" +#include "force.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "fix_neigh_history.h" +#include "comm.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairGranHookeHistoryMulti::PairGranHookeHistoryMulti(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 1; + no_virial_fdotr_compute = 1; + history = 1; + fix_history = NULL; + + single_extra = 10; + svector = new double[10]; + + neighprev = 0; + + nmax = 0; + mass_rigid = NULL; + + // set comm size needed by this Pair if used with fix rigid + + comm_forward = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairGranHookeHistoryMulti::~PairGranHookeHistoryMulti() +{ + delete [] svector; + if (fix_history) modify->delete_fix("NEIGH_HISTORY"); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(kn); + memory->destroy(kt); + memory->destroy(gamman); + memory->destroy(gammat); + memory->destroy(xmu); + memory->destroy(dampflag); + + delete [] onerad_dynamic; + delete [] onerad_frozen; + delete [] maxrad_dynamic; + delete [] maxrad_frozen; + } + + memory->destroy(mass_rigid); +} + +/* ---------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; + double radi,radj,radsum,rsq,r,rinv,rsqinv; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; + double wr1,wr2,wr3; + double vtr1,vtr2,vtr3,vrel; + double mi,mj,meff,damp,ccel,tor1,tor2,tor3; + double fn,fs,fs1,fs2,fs3; + double shrmag,rsht; + int *ilist,*jlist,*numneigh,**firstneigh; + int *touch,**firsttouch; + double *shear,*allshear,**firstshear; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + int shearupdate = 1; + if (update->setupflag) shearupdate = 0; + + // update rigid body info for owned & ghost atoms if using FixRigid masses + // body[i] = which body atom I is in, -1 if none + // mass_body = mass of each rigid body + + if (fix_rigid && neighbor->ago == 0) { + int tmp; + int *body = (int *) fix_rigid->extract("body",tmp); + double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); + if (atom->nmax > nmax) { + memory->destroy(mass_rigid); + nmax = atom->nmax; + memory->create(mass_rigid,nmax,"pair:mass_rigid"); + } + int nlocal = atom->nlocal; + for (i = 0; i < nlocal; i++) + if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; + else mass_rigid[i] = 0.0; + comm->forward_comm_pair(this); + } + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + int *type = atom->type; + double **omega = atom->omega; + double **torque = atom->torque; + double *radius = atom->radius; + double *rmass = atom->rmass; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + firsttouch = fix_history->firstflag; + firstshear = fix_history->firstvalue; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + radi = radius[i]; + touch = firsttouch[i]; + allshear = firstshear[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + radj = radius[j]; + radsum = radi + radj; + + if (rsq >= radsum*radsum) { + + // unset non-touching neighbors + + touch[jj] = 0; + shear = &allshear[3*jj]; + shear[0] = 0.0; + shear[1] = 0.0; + shear[2] = 0.0; + + } else { + r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; + + // relative translational velocity + + vr1 = v[i][0] - v[j][0]; + vr2 = v[i][1] - v[j][1]; + vr3 = v[i][2] - v[j][2]; + + // normal component + + vnnr = vr1*delx + vr2*dely + vr3*delz; + vn1 = delx*vnnr * rsqinv; + vn2 = dely*vnnr * rsqinv; + vn3 = delz*vnnr * rsqinv; + + // tangential component + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // relative rotational velocity + + wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv; + wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv; + wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv; + + // meff = effective mass of pair of particles + // if I or J part of rigid body, use body mass + // if I or J is frozen, meff is other particle + + mi = rmass[i]; + mj = rmass[j]; + if (fix_rigid) { + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; + } + + meff = mi*mj / (mi+mj); + if (mask[i] & freeze_group_bit) meff = mj; + if (mask[j] & freeze_group_bit) meff = mi; + + // normal forces = Hookian contact + normal velocity damping + + damp = meff*gamman[itype][jtype]*vnnr*rsqinv; + ccel = kn[itype][jtype]*(radsum-r)*rinv - damp; + + // relative velocities + + vtr1 = vt1 - (delz*wr2-dely*wr3); + vtr2 = vt2 - (delx*wr3-delz*wr1); + vtr3 = vt3 - (dely*wr1-delx*wr2); + vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; + vrel = sqrt(vrel); + + // shear history effects + + touch[jj] = 1; + shear = &allshear[3*jj]; + + if (shearupdate) { + shear[0] += vtr1*dt; + shear[1] += vtr2*dt; + shear[2] += vtr3*dt; + } + shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + + shear[2]*shear[2]); + + // rotate shear displacements + + rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; + rsht *= rsqinv; + if (shearupdate) { + shear[0] -= rsht*delx; + shear[1] -= rsht*dely; + shear[2] -= rsht*delz; + } + + // tangential forces = shear + tangential velocity damping + + fs1 = - (kt[itype][jtype]*shear[0] + meff*gammat[itype][jtype]*vtr1); + fs2 = - (kt[itype][jtype]*shear[1] + meff*gammat[itype][jtype]*vtr2); + fs3 = - (kt[itype][jtype]*shear[2] + meff*gammat[itype][jtype]*vtr3); + + // rescale frictional displacements and forces if needed + + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); + fn = xmu[itype][jtype] * fabs(ccel*r); + + if (fs > fn) { + if (shrmag != 0.0) { + shear[0] = (fn/fs) * (shear[0] + + meff*gammat[itype][jtype]*vtr1/kt[itype][jtype]) - + meff*gammat[itype][jtype]*vtr1/kt[itype][jtype]; + shear[1] = (fn/fs) * (shear[1] + + meff*gammat[itype][jtype]*vtr2/kt[itype][jtype]) - + meff*gammat[itype][jtype]*vtr2/kt[itype][jtype]; + shear[2] = (fn/fs) * (shear[2] + + meff*gammat[itype][jtype]*vtr3/kt[itype][jtype]) - + meff*gammat[itype][jtype]*vtr3/kt[itype][jtype]; + fs1 *= fn/fs; + fs2 *= fn/fs; + fs3 *= fn/fs; + } else fs1 = fs2 = fs3 = 0.0; + } + + // forces & torques + + fx = delx*ccel + fs1; + fy = dely*ccel + fs2; + fz = delz*ccel + fs3; + f[i][0] += fx; + f[i][1] += fy; + f[i][2] += fz; + + tor1 = rinv * (dely*fs3 - delz*fs2); + tor2 = rinv * (delz*fs1 - delx*fs3); + tor3 = rinv * (delx*fs2 - dely*fs1); + torque[i][0] -= radi*tor1; + torque[i][1] -= radi*tor2; + torque[i][2] -= radi*tor3; + + if (newton_pair || j < nlocal) { + f[j][0] -= fx; + f[j][1] -= fy; + f[j][2] -= fz; + torque[j][0] -= radj*tor1; + torque[j][1] -= radj*tor2; + torque[j][2] -= radj*tor3; + } + + if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, + 0.0,0.0,fx,fy,fz,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + 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; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(kn,n+1,n+1,"pair:kn"); + memory->create(kt,n+1,n+1,"pair:kt"); + memory->create(gamman,n+1,n+1,"pair:gamman"); + memory->create(gammat,n+1,n+1,"pair:gammat"); + memory->create(xmu,n+1,n+1,"pair:xmu"); + memory->create(dampflag,n+1,n+1,"pair:dampflag"); + + onerad_dynamic = new double[n+1]; + onerad_frozen = new double[n+1]; + maxrad_dynamic = new double[n+1]; + maxrad_frozen = new double[n+1]; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::settings(int narg, char **arg) +{ + if (narg != 1) error->all(FLERR,"Illegal pair_style command"); + + if (strcmp(arg[0],"NULL") == 0 ) cut_global = -1.0; + else cut_global = force->numeric(FLERR,arg[0]); + + // reset cutoffs that have been explicitly set + if (allocated) { + 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; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::coeff(int narg, char **arg) +{ + if (narg < 8 || narg > 9) + error->all(FLERR,"Incorrect args for pair coefficients"); + + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + double kn_one = force->numeric(FLERR,arg[2]); + double kt_one; + if (strcmp(arg[3],"NULL") == 0) kt_one = kn_one * 2.0/7.0; + else kt_one = force->numeric(FLERR,arg[3]); + + double gamman_one = force->numeric(FLERR,arg[4]); + double gammat_one; + if (strcmp(arg[5],"NULL") == 0) gammat_one = 0.5 * gamman_one; + else gammat_one = force->numeric(FLERR,arg[5]); + + double xmu_one = force->numeric(FLERR,arg[6]); + int dampflag_one = force->inumeric(FLERR,arg[7]); + if (dampflag_one == 0) gammat_one = 0.0; + + if (kn_one < 0.0 || kt_one < 0.0 || gamman_one < 0.0 || gammat_one < 0.0 || + xmu_one < 0.0 || xmu_one > 10000.0 || dampflag_one < 0 || dampflag_one > 1) + error->all(FLERR,"Illegal pair_style command"); + + // convert Kn and Kt from pressure units to force/distance^2 + kn_one /= force->nktv2p; + kt_one /= force->nktv2p; + + double cut_one = cut_global; + if (narg==9) { + if (strcmp(arg[8],"NULL") == 0) cut_one = -1.0; + else cut_one = force->numeric(FLERR,arg[8]); + } + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + kn[i][j] = kn_one; + kt[i][j] = kt_one; + gamman[i][j] = gamman_one; + gammat[i][j] = gammat_one; + xmu[i][j] = xmu_one; + dampflag[i][j] = dampflag_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::init_style() +{ + int i; + + // error and warning checks + + if (!atom->radius_flag || !atom->rmass_flag) + error->all(FLERR,"Pair granular requires atom attributes radius, rmass"); + if (comm->ghost_velocity == 0) + error->all(FLERR,"Pair granular requires ghost atoms store velocity"); + + // need a granular neigh list + + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->size = 1; + if (history) neighbor->requests[irequest]->history = 1; + + dt = update->dt; + + // if shear history is stored: + // if first init, create Fix needed for storing shear history + + if (history && fix_history == NULL) { + char dnumstr[16]; + sprintf(dnumstr,"%d",3); + char **fixarg = new char*[4]; + fixarg[0] = (char *) "NEIGH_HISTORY"; + fixarg[1] = (char *) "all"; + fixarg[2] = (char *) "NEIGH_HISTORY"; + fixarg[3] = dnumstr; + modify->add_fix(4,fixarg,1); + delete [] fixarg; + fix_history = (FixNeighHistory *) modify->fix[modify->nfix-1]; + fix_history->pair = this; + } + + // check for FixFreeze and set freeze_group_bit + + for (i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"freeze") == 0) break; + if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit; + else freeze_group_bit = 0; + + // check for FixRigid so can extract rigid body masses + + fix_rigid = NULL; + for (i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) break; + if (i < modify->nfix) fix_rigid = modify->fix[i]; + + // check for FixPour and FixDeposit so can extract particle radii + + int ipour; + for (ipour = 0; ipour < modify->nfix; ipour++) + if (strcmp(modify->fix[ipour]->style,"pour") == 0) break; + if (ipour == modify->nfix) ipour = -1; + + int idep; + for (idep = 0; idep < modify->nfix; idep++) + if (strcmp(modify->fix[idep]->style,"deposit") == 0) break; + if (idep == modify->nfix) idep = -1; + + // set maxrad_dynamic and maxrad_frozen for each type + // include future FixPour and FixDeposit particles as dynamic + + int itype; + for (i = 1; i <= atom->ntypes; i++) { + onerad_dynamic[i] = onerad_frozen[i] = 0.0; + if (ipour >= 0) { + itype = i; + onerad_dynamic[i] = + *((double *) modify->fix[ipour]->extract("radius",itype)); + } + if (idep >= 0) { + itype = i; + onerad_dynamic[i] = + *((double *) modify->fix[idep]->extract("radius",itype)); + } + } + + double *radius = atom->radius; + int *mask = atom->mask; + int *type = atom->type; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) + if (mask[i] & freeze_group_bit) + onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius[i]); + else + onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]); + + MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes, + MPI_DOUBLE,MPI_MAX,world); + MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes, + MPI_DOUBLE,MPI_MAX,world); + + // set fix which stores history info + + if (history) { + int ifix = modify->find_fix("NEIGH_HISTORY"); + if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID"); + fix_history = (FixNeighHistory *) modify->fix[ifix]; + } +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairGranHookeHistoryMulti::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + kn[i][j] = mix_stiffness(kn[i][i],kn[j][j]); + kt[i][j] = mix_stiffness(kt[i][i],kt[j][j]); + gamman[i][j] = mix_damping(gamman[i][i],gamman[j][j]); + gammat[i][j] = mix_damping(gammat[i][i],gammat[j][j]); + xmu[i][j] = mix_friction(xmu[i][i],xmu[j][j]); + + dampflag[i][j] = 0; + if (dampflag[i][i] || dampflag[j][j]) dampflag[i][j] = 1; + + } + + kn[j][i] = kn[i][j]; + kt[j][i] = kt[i][j]; + gamman[j][i] = gamman[i][j]; + gammat[j][i] = gammat[i][j]; + xmu[j][i] = xmu[i][j]; + dampflag[j][i] = dampflag[i][j]; + + // cutoff = sum of max I,J radii for + // dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen + double cutoff = cut[i][j]; + if (cut[i][j] < 0.0) cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; + if (cut[i][j] < 0.0) cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); + if (cut[i][j] < 0.0) cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); + + // cut[i][j] = MAX(maxrad_dynamic[i],maxrad_dynamic[j]); + // double cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; + // cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); + // cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); + return cutoff; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + 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); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::write_restart_settings(FILE *fp) +{ + fwrite(&kn,sizeof(double),1,fp); + fwrite(&kt,sizeof(double),1,fp); + fwrite(&gamman,sizeof(double),1,fp); + fwrite(&gammat,sizeof(double),1,fp); + fwrite(&xmu,sizeof(double),1,fp); + fwrite(&dampflag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&kn,sizeof(double),1,fp); + fread(&kt,sizeof(double),1,fp); + fread(&gamman,sizeof(double),1,fp); + fread(&gammat,sizeof(double),1,fp); + fread(&xmu,sizeof(double),1,fp); + fread(&dampflag,sizeof(int),1,fp); + } + MPI_Bcast(&kn,1,MPI_DOUBLE,0,world); + MPI_Bcast(&kt,1,MPI_DOUBLE,0,world); + MPI_Bcast(&gamman,1,MPI_DOUBLE,0,world); + MPI_Bcast(&gammat,1,MPI_DOUBLE,0,world); + MPI_Bcast(&xmu,1,MPI_DOUBLE,0,world); + MPI_Bcast(&dampflag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::reset_dt() +{ + dt = update->dt; +} + +/* ---------------------------------------------------------------------- */ + +double PairGranHookeHistoryMulti::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double radi,radj,radsum; + double r,rinv,rsqinv,delx,dely,delz; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; + double mi,mj,meff,damp,ccel; + double vtr1,vtr2,vtr3,vrel,shrmag,rsht; + double fs1,fs2,fs3,fs,fn; + + double *radius = atom->radius; + radi = radius[i]; + radj = radius[j]; + radsum = radi + radj; + + if (rsq >= radsum*radsum) { + fforce = 0.0; + for (int m = 0; m < single_extra; m++) svector[m] = 0.0; + return 0.0; + } + + r = sqrt(rsq); + rinv = 1.0/r; + rsqinv = 1.0/rsq; + + // relative translational velocity + + double **v = atom->v; + vr1 = v[i][0] - v[j][0]; + vr2 = v[i][1] - v[j][1]; + vr3 = v[i][2] - v[j][2]; + + // normal component + + double **x = atom->x; + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + + vnnr = vr1*delx + vr2*dely + vr3*delz; + vn1 = delx*vnnr * rsqinv; + vn2 = dely*vnnr * rsqinv; + vn3 = delz*vnnr * rsqinv; + + // tangential component + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // relative rotational velocity + + double **omega = atom->omega; + wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv; + wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv; + wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv; + + // meff = effective mass of pair of particles + // if I or J part of rigid body, use body mass + // if I or J is frozen, meff is other particle + + double *rmass = atom->rmass; + int *mask = atom->mask; + + mi = rmass[i]; + mj = rmass[j]; + if (fix_rigid) { + // NOTE: insure mass_rigid is current for owned+ghost atoms? + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; + } + + meff = mi*mj / (mi+mj); + if (mask[i] & freeze_group_bit) meff = mj; + if (mask[j] & freeze_group_bit) meff = mi; + + // normal forces = Hookian contact + normal velocity damping + + damp = meff*gamman[itype][jtype]*vnnr*rsqinv; + ccel = kn[itype][jtype]*(radsum-r)*rinv - damp; + + // relative velocities + + vtr1 = vt1 - (delz*wr2-dely*wr3); + vtr2 = vt2 - (delx*wr3-delz*wr1); + vtr3 = vt3 - (dely*wr1-delx*wr2); + vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; + vrel = sqrt(vrel); + + // shear history effects + // neighprev = index of found neigh on previous call + // search entire jnum list of neighbors of I for neighbor J + // start from neighprev, since will typically be next neighbor + // reset neighprev to 0 as necessary + + int jnum = list->numneigh[i]; + int *jlist = list->firstneigh[i]; + double *allshear = fix_history->firstvalue[i]; + + for (int jj = 0; jj < jnum; jj++) { + neighprev++; + if (neighprev >= jnum) neighprev = 0; + if (jlist[neighprev] == j) break; + } + + double *shear = &allshear[3*neighprev]; + shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + + shear[2]*shear[2]); + + // rotate shear displacements + + rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; + rsht *= rsqinv; + + // tangential forces = shear + tangential velocity damping + + fs1 = - (kt[itype][jtype]*shear[0] + meff*gammat[itype][jtype]*vtr1); + fs2 = - (kt[itype][jtype]*shear[1] + meff*gammat[itype][jtype]*vtr2); + fs3 = - (kt[itype][jtype]*shear[2] + meff*gammat[itype][jtype]*vtr3); + + // rescale frictional displacements and forces if needed + + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); + fn = xmu[itype][jtype] * fabs(ccel*r); + + if (fs > fn) { + if (shrmag != 0.0) { + fs1 *= fn/fs; + fs2 *= fn/fs; + fs3 *= fn/fs; + fs *= fn/fs; + } else fs1 = fs2 = fs3 = fs = 0.0; + } + + // set force and return no energy + + fforce = ccel; + + // set single_extra quantities + + svector[0] = fs1; + svector[1] = fs2; + svector[2] = fs3; + svector[3] = fs; + svector[4] = vn1; + svector[5] = vn2; + svector[6] = vn3; + svector[7] = vt1; + svector[8] = vt2; + svector[9] = vt3; + + return 0.0; +} + +/* ---------------------------------------------------------------------- */ + +int PairGranHookeHistoryMulti::pack_forward_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = mass_rigid[j]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void PairGranHookeHistoryMulti::unpack_forward_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) + mass_rigid[i] = buf[m++]; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairGranHookeHistoryMulti::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + mixing of stiffness +------------------------------------------------------------------------- */ + +double PairGranHookeHistoryMulti::mix_stiffness(double kii, double kjj) +{ + return kii*kjj/(kii + kjj); +} + +/* ---------------------------------------------------------------------- + mixing of damping +------------------------------------------------------------------------- */ + +double PairGranHookeHistoryMulti::mix_damping(double gammaii, double gammajj) +{ + return sqrt(gammaii*gammajj); +} + +/* ---------------------------------------------------------------------- + mixing of friction +------------------------------------------------------------------------- */ + +double PairGranHookeHistoryMulti::mix_friction(double xmuii, double xmujj) +{ + return MAX(xmuii,xmujj); +} + diff --git a/src/from_ishan/pair_gran_hooke_history_multi.h b/src/from_ishan/pair_gran_hooke_history_multi.h new file mode 100644 index 0000000000..f302ede96c --- /dev/null +++ b/src/from_ishan/pair_gran_hooke_history_multi.h @@ -0,0 +1,109 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(gran/hooke/history/multi,PairGranHookeHistoryMulti) + +#else + +#ifndef LMP_PAIR_GRAN_HOOKE_HISTORY_MULTI_H +#define LMP_PAIR_GRAN_HOOKE_HISTORY_MULTI_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairGranHookeHistoryMulti : public Pair { + public: + PairGranHookeHistoryMulti(class LAMMPS *); + virtual ~PairGranHookeHistoryMulti(); + virtual void compute(int, int); + virtual void settings(int, char **); + virtual void coeff(int, char **); // Made Virtual by IS Oct 7 2017 + void init_style(); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *); + void read_restart_settings(FILE *); + void reset_dt(); + virtual double single(int, int, int, int, double, double, double, double &); + int pack_forward_comm(int, int *, double *, int, int *); + void unpack_forward_comm(int, int, double *); + double memory_usage(); + + protected: + double cut_global; + double **kn,**kt,**gamman,**gammat,**xmu,**cut; + int **dampflag; + double dt; + int freeze_group_bit; + int history; + + int neighprev; + double *onerad_dynamic,*onerad_frozen; + double *maxrad_dynamic,*maxrad_frozen; + + class FixNeighHistory *fix_history; + + // storage of rigid body masses for use in granular interactions + + class Fix *fix_rigid; // ptr to rigid body fix, NULL if none + double *mass_rigid; // rigid mass for owned+ghost atoms + int nmax; // allocated size of mass_rigid + + virtual void allocate(); // Made Virtual by IS Oct 7 2017 + +private: + double mix_stiffness(double kii, double kjj); + double mix_damping(double gammaii, double gammajj); + double mix_friction(double xmuii, double xmujj); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair granular requires atom attributes radius, rmass + +The atom style defined does not have these attributes. + +E: Pair granular requires ghost atoms store velocity + +Use the comm_modify vel yes command to enable this. + +E: Pair granular with shear history requires newton pair off + +This is a current restriction of the implementation of pair +granular styles with history. + +E: Could not find pair fix ID + +A fix is created internally by the pair style to store shear +history information. You cannot delete it. + +*/ From dfb5cd3262bfa0a45eb4989391c1a62d2756232d Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Fri, 2 Oct 2020 16:37:01 -0600 Subject: [PATCH 02/25] merged with mastered, copied Ishan's files --- src/{from_ishan => }/npair_half_size_multi_newtoff.cpp | 0 src/{from_ishan => }/npair_half_size_multi_newtoff.h | 0 src/{from_ishan => }/npair_half_size_multi_newton.cpp | 0 src/{from_ishan => }/npair_half_size_multi_newton.h | 0 src/{from_ishan => }/npair_half_size_multi_newton_tri.cpp | 0 src/{from_ishan => }/npair_half_size_multi_newton_tri.h | 0 src/{from_ishan => }/pair_gran_hooke_history_multi.cpp | 0 src/{from_ishan => }/pair_gran_hooke_history_multi.h | 0 8 files changed, 0 insertions(+), 0 deletions(-) rename src/{from_ishan => }/npair_half_size_multi_newtoff.cpp (100%) rename src/{from_ishan => }/npair_half_size_multi_newtoff.h (100%) rename src/{from_ishan => }/npair_half_size_multi_newton.cpp (100%) rename src/{from_ishan => }/npair_half_size_multi_newton.h (100%) rename src/{from_ishan => }/npair_half_size_multi_newton_tri.cpp (100%) rename src/{from_ishan => }/npair_half_size_multi_newton_tri.h (100%) rename src/{from_ishan => }/pair_gran_hooke_history_multi.cpp (100%) rename src/{from_ishan => }/pair_gran_hooke_history_multi.h (100%) diff --git a/src/from_ishan/npair_half_size_multi_newtoff.cpp b/src/npair_half_size_multi_newtoff.cpp similarity index 100% rename from src/from_ishan/npair_half_size_multi_newtoff.cpp rename to src/npair_half_size_multi_newtoff.cpp diff --git a/src/from_ishan/npair_half_size_multi_newtoff.h b/src/npair_half_size_multi_newtoff.h similarity index 100% rename from src/from_ishan/npair_half_size_multi_newtoff.h rename to src/npair_half_size_multi_newtoff.h diff --git a/src/from_ishan/npair_half_size_multi_newton.cpp b/src/npair_half_size_multi_newton.cpp similarity index 100% rename from src/from_ishan/npair_half_size_multi_newton.cpp rename to src/npair_half_size_multi_newton.cpp diff --git a/src/from_ishan/npair_half_size_multi_newton.h b/src/npair_half_size_multi_newton.h similarity index 100% rename from src/from_ishan/npair_half_size_multi_newton.h rename to src/npair_half_size_multi_newton.h diff --git a/src/from_ishan/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp similarity index 100% rename from src/from_ishan/npair_half_size_multi_newton_tri.cpp rename to src/npair_half_size_multi_newton_tri.cpp diff --git a/src/from_ishan/npair_half_size_multi_newton_tri.h b/src/npair_half_size_multi_newton_tri.h similarity index 100% rename from src/from_ishan/npair_half_size_multi_newton_tri.h rename to src/npair_half_size_multi_newton_tri.h diff --git a/src/from_ishan/pair_gran_hooke_history_multi.cpp b/src/pair_gran_hooke_history_multi.cpp similarity index 100% rename from src/from_ishan/pair_gran_hooke_history_multi.cpp rename to src/pair_gran_hooke_history_multi.cpp diff --git a/src/from_ishan/pair_gran_hooke_history_multi.h b/src/pair_gran_hooke_history_multi.h similarity index 100% rename from src/from_ishan/pair_gran_hooke_history_multi.h rename to src/pair_gran_hooke_history_multi.h From bc1c16d1a632e4ec5a0b3efaa0b88a92c286e33e Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Mon, 5 Oct 2020 16:14:01 -0600 Subject: [PATCH 03/25] Removing gran/multi pairstyles --- src/pair_gran_hooke_history_multi.cpp | 885 -------------------------- src/pair_gran_hooke_history_multi.h | 109 ---- 2 files changed, 994 deletions(-) delete mode 100644 src/pair_gran_hooke_history_multi.cpp delete mode 100644 src/pair_gran_hooke_history_multi.h diff --git a/src/pair_gran_hooke_history_multi.cpp b/src/pair_gran_hooke_history_multi.cpp deleted file mode 100644 index 83f5221741..0000000000 --- a/src/pair_gran_hooke_history_multi.cpp +++ /dev/null @@ -1,885 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing authors: Leo Silbert (SNL), Gary Grest (SNL) -------------------------------------------------------------------------- */ - -#include -#include -#include -#include -#include "pair_gran_hooke_history_multi.h" -#include "atom.h" -#include "atom_vec.h" -#include "domain.h" -#include "force.h" -#include "update.h" -#include "modify.h" -#include "fix.h" -#include "fix_neigh_history.h" -#include "comm.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -PairGranHookeHistoryMulti::PairGranHookeHistoryMulti(LAMMPS *lmp) : Pair(lmp) -{ - single_enable = 1; - no_virial_fdotr_compute = 1; - history = 1; - fix_history = NULL; - - single_extra = 10; - svector = new double[10]; - - neighprev = 0; - - nmax = 0; - mass_rigid = NULL; - - // set comm size needed by this Pair if used with fix rigid - - comm_forward = 1; -} - -/* ---------------------------------------------------------------------- */ - -PairGranHookeHistoryMulti::~PairGranHookeHistoryMulti() -{ - delete [] svector; - if (fix_history) modify->delete_fix("NEIGH_HISTORY"); - - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - - memory->destroy(kn); - memory->destroy(kt); - memory->destroy(gamman); - memory->destroy(gammat); - memory->destroy(xmu); - memory->destroy(dampflag); - - delete [] onerad_dynamic; - delete [] onerad_frozen; - delete [] maxrad_dynamic; - delete [] maxrad_frozen; - } - - memory->destroy(mass_rigid); -} - -/* ---------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; - double radi,radj,radsum,rsq,r,rinv,rsqinv; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; - double wr1,wr2,wr3; - double vtr1,vtr2,vtr3,vrel; - double mi,mj,meff,damp,ccel,tor1,tor2,tor3; - double fn,fs,fs1,fs2,fs3; - double shrmag,rsht; - int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *shear,*allshear,**firstshear; - - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - int shearupdate = 1; - if (update->setupflag) shearupdate = 0; - - // update rigid body info for owned & ghost atoms if using FixRigid masses - // body[i] = which body atom I is in, -1 if none - // mass_body = mass of each rigid body - - if (fix_rigid && neighbor->ago == 0) { - int tmp; - int *body = (int *) fix_rigid->extract("body",tmp); - double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); - if (atom->nmax > nmax) { - memory->destroy(mass_rigid); - nmax = atom->nmax; - memory->create(mass_rigid,nmax,"pair:mass_rigid"); - } - int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) - if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; - else mass_rigid[i] = 0.0; - comm->forward_comm_pair(this); - } - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - int *type = atom->type; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - int *mask = atom->mask; - int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - firsttouch = fix_history->firstflag; - firstshear = fix_history->firstvalue; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - radi = radius[i]; - touch = firsttouch[i]; - allshear = firstshear[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - radj = radius[j]; - radsum = radi + radj; - - if (rsq >= radsum*radsum) { - - // unset non-touching neighbors - - touch[jj] = 0; - shear = &allshear[3*jj]; - shear[0] = 0.0; - shear[1] = 0.0; - shear[2] = 0.0; - - } else { - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - - // relative translational velocity - - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - vnnr = vr1*delx + vr2*dely + vr3*delz; - vn1 = delx*vnnr * rsqinv; - vn2 = dely*vnnr * rsqinv; - vn3 = delz*vnnr * rsqinv; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv; - wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv; - wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv; - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - // normal forces = Hookian contact + normal velocity damping - - damp = meff*gamman[itype][jtype]*vnnr*rsqinv; - ccel = kn[itype][jtype]*(radsum-r)*rinv - damp; - - // relative velocities - - vtr1 = vt1 - (delz*wr2-dely*wr3); - vtr2 = vt2 - (delx*wr3-delz*wr1); - vtr3 = vt3 - (dely*wr1-delx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - - touch[jj] = 1; - shear = &allshear[3*jj]; - - if (shearupdate) { - shear[0] += vtr1*dt; - shear[1] += vtr2*dt; - shear[2] += vtr3*dt; - } - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // rotate shear displacements - - rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; - rsht *= rsqinv; - if (shearupdate) { - shear[0] -= rsht*delx; - shear[1] -= rsht*dely; - shear[2] -= rsht*delz; - } - - // tangential forces = shear + tangential velocity damping - - fs1 = - (kt[itype][jtype]*shear[0] + meff*gammat[itype][jtype]*vtr1); - fs2 = - (kt[itype][jtype]*shear[1] + meff*gammat[itype][jtype]*vtr2); - fs3 = - (kt[itype][jtype]*shear[2] + meff*gammat[itype][jtype]*vtr3); - - // rescale frictional displacements and forces if needed - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - fn = xmu[itype][jtype] * fabs(ccel*r); - - if (fs > fn) { - if (shrmag != 0.0) { - shear[0] = (fn/fs) * (shear[0] + - meff*gammat[itype][jtype]*vtr1/kt[itype][jtype]) - - meff*gammat[itype][jtype]*vtr1/kt[itype][jtype]; - shear[1] = (fn/fs) * (shear[1] + - meff*gammat[itype][jtype]*vtr2/kt[itype][jtype]) - - meff*gammat[itype][jtype]*vtr2/kt[itype][jtype]; - shear[2] = (fn/fs) * (shear[2] + - meff*gammat[itype][jtype]*vtr3/kt[itype][jtype]) - - meff*gammat[itype][jtype]*vtr3/kt[itype][jtype]; - fs1 *= fn/fs; - fs2 *= fn/fs; - fs3 *= fn/fs; - } else fs1 = fs2 = fs3 = 0.0; - } - - // forces & torques - - fx = delx*ccel + fs1; - fy = dely*ccel + fs2; - fz = delz*ccel + fs3; - f[i][0] += fx; - f[i][1] += fy; - f[i][2] += fz; - - tor1 = rinv * (dely*fs3 - delz*fs2); - tor2 = rinv * (delz*fs1 - delx*fs3); - tor3 = rinv * (delx*fs2 - dely*fs1); - torque[i][0] -= radi*tor1; - torque[i][1] -= radi*tor2; - torque[i][2] -= radi*tor3; - - if (newton_pair || j < nlocal) { - f[j][0] -= fx; - f[j][1] -= fy; - f[j][2] -= fz; - torque[j][0] -= radj*tor1; - torque[j][1] -= radj*tor2; - torque[j][2] -= radj*tor3; - } - - if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, - 0.0,0.0,fx,fy,fz,delx,dely,delz); - } - } - } - - if (vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::allocate() -{ - allocated = 1; - int n = atom->ntypes; - - 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; - - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(kn,n+1,n+1,"pair:kn"); - memory->create(kt,n+1,n+1,"pair:kt"); - memory->create(gamman,n+1,n+1,"pair:gamman"); - memory->create(gammat,n+1,n+1,"pair:gammat"); - memory->create(xmu,n+1,n+1,"pair:xmu"); - memory->create(dampflag,n+1,n+1,"pair:dampflag"); - - onerad_dynamic = new double[n+1]; - onerad_frozen = new double[n+1]; - maxrad_dynamic = new double[n+1]; - maxrad_frozen = new double[n+1]; -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::settings(int narg, char **arg) -{ - if (narg != 1) error->all(FLERR,"Illegal pair_style command"); - - if (strcmp(arg[0],"NULL") == 0 ) cut_global = -1.0; - else cut_global = force->numeric(FLERR,arg[0]); - - // reset cutoffs that have been explicitly set - if (allocated) { - 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; - } -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::coeff(int narg, char **arg) -{ - if (narg < 8 || narg > 9) - error->all(FLERR,"Incorrect args for pair coefficients"); - - if (!allocated) allocate(); - - int ilo,ihi,jlo,jhi; - force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); - force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); - - double kn_one = force->numeric(FLERR,arg[2]); - double kt_one; - if (strcmp(arg[3],"NULL") == 0) kt_one = kn_one * 2.0/7.0; - else kt_one = force->numeric(FLERR,arg[3]); - - double gamman_one = force->numeric(FLERR,arg[4]); - double gammat_one; - if (strcmp(arg[5],"NULL") == 0) gammat_one = 0.5 * gamman_one; - else gammat_one = force->numeric(FLERR,arg[5]); - - double xmu_one = force->numeric(FLERR,arg[6]); - int dampflag_one = force->inumeric(FLERR,arg[7]); - if (dampflag_one == 0) gammat_one = 0.0; - - if (kn_one < 0.0 || kt_one < 0.0 || gamman_one < 0.0 || gammat_one < 0.0 || - xmu_one < 0.0 || xmu_one > 10000.0 || dampflag_one < 0 || dampflag_one > 1) - error->all(FLERR,"Illegal pair_style command"); - - // convert Kn and Kt from pressure units to force/distance^2 - kn_one /= force->nktv2p; - kt_one /= force->nktv2p; - - double cut_one = cut_global; - if (narg==9) { - if (strcmp(arg[8],"NULL") == 0) cut_one = -1.0; - else cut_one = force->numeric(FLERR,arg[8]); - } - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - kn[i][j] = kn_one; - kt[i][j] = kt_one; - gamman[i][j] = gamman_one; - gammat[i][j] = gammat_one; - xmu[i][j] = xmu_one; - dampflag[i][j] = dampflag_one; - cut[i][j] = cut_one; - setflag[i][j] = 1; - count++; - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::init_style() -{ - int i; - - // error and warning checks - - if (!atom->radius_flag || !atom->rmass_flag) - error->all(FLERR,"Pair granular requires atom attributes radius, rmass"); - if (comm->ghost_velocity == 0) - error->all(FLERR,"Pair granular requires ghost atoms store velocity"); - - // need a granular neigh list - - int irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->size = 1; - if (history) neighbor->requests[irequest]->history = 1; - - dt = update->dt; - - // if shear history is stored: - // if first init, create Fix needed for storing shear history - - if (history && fix_history == NULL) { - char dnumstr[16]; - sprintf(dnumstr,"%d",3); - char **fixarg = new char*[4]; - fixarg[0] = (char *) "NEIGH_HISTORY"; - fixarg[1] = (char *) "all"; - fixarg[2] = (char *) "NEIGH_HISTORY"; - fixarg[3] = dnumstr; - modify->add_fix(4,fixarg,1); - delete [] fixarg; - fix_history = (FixNeighHistory *) modify->fix[modify->nfix-1]; - fix_history->pair = this; - } - - // check for FixFreeze and set freeze_group_bit - - for (i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"freeze") == 0) break; - if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit; - else freeze_group_bit = 0; - - // check for FixRigid so can extract rigid body masses - - fix_rigid = NULL; - for (i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) break; - if (i < modify->nfix) fix_rigid = modify->fix[i]; - - // check for FixPour and FixDeposit so can extract particle radii - - int ipour; - for (ipour = 0; ipour < modify->nfix; ipour++) - if (strcmp(modify->fix[ipour]->style,"pour") == 0) break; - if (ipour == modify->nfix) ipour = -1; - - int idep; - for (idep = 0; idep < modify->nfix; idep++) - if (strcmp(modify->fix[idep]->style,"deposit") == 0) break; - if (idep == modify->nfix) idep = -1; - - // set maxrad_dynamic and maxrad_frozen for each type - // include future FixPour and FixDeposit particles as dynamic - - int itype; - for (i = 1; i <= atom->ntypes; i++) { - onerad_dynamic[i] = onerad_frozen[i] = 0.0; - if (ipour >= 0) { - itype = i; - onerad_dynamic[i] = - *((double *) modify->fix[ipour]->extract("radius",itype)); - } - if (idep >= 0) { - itype = i; - onerad_dynamic[i] = - *((double *) modify->fix[idep]->extract("radius",itype)); - } - } - - double *radius = atom->radius; - int *mask = atom->mask; - int *type = atom->type; - int nlocal = atom->nlocal; - - for (i = 0; i < nlocal; i++) - if (mask[i] & freeze_group_bit) - onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius[i]); - else - onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]); - - MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes, - MPI_DOUBLE,MPI_MAX,world); - MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes, - MPI_DOUBLE,MPI_MAX,world); - - // set fix which stores history info - - if (history) { - int ifix = modify->find_fix("NEIGH_HISTORY"); - if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID"); - fix_history = (FixNeighHistory *) modify->fix[ifix]; - } -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::init_one(int i, int j) -{ - if (setflag[i][j] == 0) { - kn[i][j] = mix_stiffness(kn[i][i],kn[j][j]); - kt[i][j] = mix_stiffness(kt[i][i],kt[j][j]); - gamman[i][j] = mix_damping(gamman[i][i],gamman[j][j]); - gammat[i][j] = mix_damping(gammat[i][i],gammat[j][j]); - xmu[i][j] = mix_friction(xmu[i][i],xmu[j][j]); - - dampflag[i][j] = 0; - if (dampflag[i][i] || dampflag[j][j]) dampflag[i][j] = 1; - - } - - kn[j][i] = kn[i][j]; - kt[j][i] = kt[i][j]; - gamman[j][i] = gamman[i][j]; - gammat[j][i] = gammat[i][j]; - xmu[j][i] = xmu[i][j]; - dampflag[j][i] = dampflag[i][j]; - - // cutoff = sum of max I,J radii for - // dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen - double cutoff = cut[i][j]; - if (cut[i][j] < 0.0) cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; - if (cut[i][j] < 0.0) cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); - if (cut[i][j] < 0.0) cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); - - // cut[i][j] = MAX(maxrad_dynamic[i],maxrad_dynamic[j]); - // double cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; - // cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); - // cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); - return cutoff; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::write_restart(FILE *fp) -{ - write_restart_settings(fp); - - 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); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::read_restart(FILE *fp) -{ - read_restart_settings(fp); - allocate(); - - int i,j; - int me = comm->me; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) { - if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::write_restart_settings(FILE *fp) -{ - fwrite(&kn,sizeof(double),1,fp); - fwrite(&kt,sizeof(double),1,fp); - fwrite(&gamman,sizeof(double),1,fp); - fwrite(&gammat,sizeof(double),1,fp); - fwrite(&xmu,sizeof(double),1,fp); - fwrite(&dampflag,sizeof(int),1,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::read_restart_settings(FILE *fp) -{ - if (comm->me == 0) { - fread(&kn,sizeof(double),1,fp); - fread(&kt,sizeof(double),1,fp); - fread(&gamman,sizeof(double),1,fp); - fread(&gammat,sizeof(double),1,fp); - fread(&xmu,sizeof(double),1,fp); - fread(&dampflag,sizeof(int),1,fp); - } - MPI_Bcast(&kn,1,MPI_DOUBLE,0,world); - MPI_Bcast(&kt,1,MPI_DOUBLE,0,world); - MPI_Bcast(&gamman,1,MPI_DOUBLE,0,world); - MPI_Bcast(&gammat,1,MPI_DOUBLE,0,world); - MPI_Bcast(&xmu,1,MPI_DOUBLE,0,world); - MPI_Bcast(&dampflag,1,MPI_INT,0,world); -} - -/* ---------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::reset_dt() -{ - dt = update->dt; -} - -/* ---------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::single(int i, int j, int itype, int jtype, - double rsq, - double factor_coul, double factor_lj, - double &fforce) -{ - double radi,radj,radsum; - double r,rinv,rsqinv,delx,dely,delz; - double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; - double mi,mj,meff,damp,ccel; - double vtr1,vtr2,vtr3,vrel,shrmag,rsht; - double fs1,fs2,fs3,fs,fn; - - double *radius = atom->radius; - radi = radius[i]; - radj = radius[j]; - radsum = radi + radj; - - if (rsq >= radsum*radsum) { - fforce = 0.0; - for (int m = 0; m < single_extra; m++) svector[m] = 0.0; - return 0.0; - } - - r = sqrt(rsq); - rinv = 1.0/r; - rsqinv = 1.0/rsq; - - // relative translational velocity - - double **v = atom->v; - vr1 = v[i][0] - v[j][0]; - vr2 = v[i][1] - v[j][1]; - vr3 = v[i][2] - v[j][2]; - - // normal component - - double **x = atom->x; - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - - vnnr = vr1*delx + vr2*dely + vr3*delz; - vn1 = delx*vnnr * rsqinv; - vn2 = dely*vnnr * rsqinv; - vn3 = delz*vnnr * rsqinv; - - // tangential component - - vt1 = vr1 - vn1; - vt2 = vr2 - vn2; - vt3 = vr3 - vn3; - - // relative rotational velocity - - double **omega = atom->omega; - wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv; - wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv; - wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv; - - // meff = effective mass of pair of particles - // if I or J part of rigid body, use body mass - // if I or J is frozen, meff is other particle - - double *rmass = atom->rmass; - int *mask = atom->mask; - - mi = rmass[i]; - mj = rmass[j]; - if (fix_rigid) { - // NOTE: insure mass_rigid is current for owned+ghost atoms? - if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; - if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; - } - - meff = mi*mj / (mi+mj); - if (mask[i] & freeze_group_bit) meff = mj; - if (mask[j] & freeze_group_bit) meff = mi; - - // normal forces = Hookian contact + normal velocity damping - - damp = meff*gamman[itype][jtype]*vnnr*rsqinv; - ccel = kn[itype][jtype]*(radsum-r)*rinv - damp; - - // relative velocities - - vtr1 = vt1 - (delz*wr2-dely*wr3); - vtr2 = vt2 - (delx*wr3-delz*wr1); - vtr3 = vt3 - (dely*wr1-delx*wr2); - vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; - vrel = sqrt(vrel); - - // shear history effects - // neighprev = index of found neigh on previous call - // search entire jnum list of neighbors of I for neighbor J - // start from neighprev, since will typically be next neighbor - // reset neighprev to 0 as necessary - - int jnum = list->numneigh[i]; - int *jlist = list->firstneigh[i]; - double *allshear = fix_history->firstvalue[i]; - - for (int jj = 0; jj < jnum; jj++) { - neighprev++; - if (neighprev >= jnum) neighprev = 0; - if (jlist[neighprev] == j) break; - } - - double *shear = &allshear[3*neighprev]; - shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + - shear[2]*shear[2]); - - // rotate shear displacements - - rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz; - rsht *= rsqinv; - - // tangential forces = shear + tangential velocity damping - - fs1 = - (kt[itype][jtype]*shear[0] + meff*gammat[itype][jtype]*vtr1); - fs2 = - (kt[itype][jtype]*shear[1] + meff*gammat[itype][jtype]*vtr2); - fs3 = - (kt[itype][jtype]*shear[2] + meff*gammat[itype][jtype]*vtr3); - - // rescale frictional displacements and forces if needed - - fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); - fn = xmu[itype][jtype] * fabs(ccel*r); - - if (fs > fn) { - if (shrmag != 0.0) { - fs1 *= fn/fs; - fs2 *= fn/fs; - fs3 *= fn/fs; - fs *= fn/fs; - } else fs1 = fs2 = fs3 = fs = 0.0; - } - - // set force and return no energy - - fforce = ccel; - - // set single_extra quantities - - svector[0] = fs1; - svector[1] = fs2; - svector[2] = fs3; - svector[3] = fs; - svector[4] = vn1; - svector[5] = vn2; - svector[6] = vn3; - svector[7] = vt1; - svector[8] = vt2; - svector[9] = vt3; - - return 0.0; -} - -/* ---------------------------------------------------------------------- */ - -int PairGranHookeHistoryMulti::pack_forward_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) -{ - int i,j,m; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - buf[m++] = mass_rigid[j]; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void PairGranHookeHistoryMulti::unpack_forward_comm(int n, int first, double *buf) -{ - int i,m,last; - - m = 0; - last = first + n; - for (i = first; i < last; i++) - mass_rigid[i] = buf[m++]; -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based arrays -------------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::memory_usage() -{ - double bytes = nmax * sizeof(double); - return bytes; -} - -/* ---------------------------------------------------------------------- - mixing of stiffness -------------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::mix_stiffness(double kii, double kjj) -{ - return kii*kjj/(kii + kjj); -} - -/* ---------------------------------------------------------------------- - mixing of damping -------------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::mix_damping(double gammaii, double gammajj) -{ - return sqrt(gammaii*gammajj); -} - -/* ---------------------------------------------------------------------- - mixing of friction -------------------------------------------------------------------------- */ - -double PairGranHookeHistoryMulti::mix_friction(double xmuii, double xmujj) -{ - return MAX(xmuii,xmujj); -} - diff --git a/src/pair_gran_hooke_history_multi.h b/src/pair_gran_hooke_history_multi.h deleted file mode 100644 index f302ede96c..0000000000 --- a/src/pair_gran_hooke_history_multi.h +++ /dev/null @@ -1,109 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef PAIR_CLASS - -PairStyle(gran/hooke/history/multi,PairGranHookeHistoryMulti) - -#else - -#ifndef LMP_PAIR_GRAN_HOOKE_HISTORY_MULTI_H -#define LMP_PAIR_GRAN_HOOKE_HISTORY_MULTI_H - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairGranHookeHistoryMulti : public Pair { - public: - PairGranHookeHistoryMulti(class LAMMPS *); - virtual ~PairGranHookeHistoryMulti(); - virtual void compute(int, int); - virtual void settings(int, char **); - virtual void coeff(int, char **); // Made Virtual by IS Oct 7 2017 - void init_style(); - double init_one(int, int); - void write_restart(FILE *); - void read_restart(FILE *); - void write_restart_settings(FILE *); - void read_restart_settings(FILE *); - void reset_dt(); - virtual double single(int, int, int, int, double, double, double, double &); - int pack_forward_comm(int, int *, double *, int, int *); - void unpack_forward_comm(int, int, double *); - double memory_usage(); - - protected: - double cut_global; - double **kn,**kt,**gamman,**gammat,**xmu,**cut; - int **dampflag; - double dt; - int freeze_group_bit; - int history; - - int neighprev; - double *onerad_dynamic,*onerad_frozen; - double *maxrad_dynamic,*maxrad_frozen; - - class FixNeighHistory *fix_history; - - // storage of rigid body masses for use in granular interactions - - class Fix *fix_rigid; // ptr to rigid body fix, NULL if none - double *mass_rigid; // rigid mass for owned+ghost atoms - int nmax; // allocated size of mass_rigid - - virtual void allocate(); // Made Virtual by IS Oct 7 2017 - -private: - double mix_stiffness(double kii, double kjj); - double mix_damping(double gammaii, double gammajj); - double mix_friction(double xmuii, double xmujj); -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: Incorrect args for pair coefficients - -Self-explanatory. Check the input script or data file. - -E: Pair granular requires atom attributes radius, rmass - -The atom style defined does not have these attributes. - -E: Pair granular requires ghost atoms store velocity - -Use the comm_modify vel yes command to enable this. - -E: Pair granular with shear history requires newton pair off - -This is a current restriction of the implementation of pair -granular styles with history. - -E: Could not find pair fix ID - -A fix is created internally by the pair style to store shear -history information. You cannot delete it. - -*/ From dc86c37e2380e516a09e7569cc8564295f21e5bb Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Fri, 9 Oct 2020 11:38:18 -0600 Subject: [PATCH 04/25] Minor updates and documentation --- doc/src/comm_modify.rst | 6 ++++-- src/npair_half_size_multi_newtoff.cpp | 6 ++---- src/npair_half_size_multi_newtoff.h | 3 +++ src/npair_half_size_multi_newton.cpp | 4 +--- src/npair_half_size_multi_newton.h | 3 +++ src/npair_half_size_multi_newton_tri.cpp | 5 +++-- src/npair_half_size_multi_newton_tri.h | 3 +++ 7 files changed, 19 insertions(+), 11 deletions(-) diff --git a/doc/src/comm_modify.rst b/doc/src/comm_modify.rst index 9a2ae60f1e..68d40281c1 100644 --- a/doc/src/comm_modify.rst +++ b/doc/src/comm_modify.rst @@ -84,13 +84,15 @@ information is available, then also a heuristic based on that bond length is computed. It is used as communication cutoff, if there is no pair style present and no *comm_modify cutoff* command used. Otherwise a warning is printed, if this bond based estimate is larger than the -communication cutoff used. A +communication cutoff used. The *cutoff/multi* option is equivalent to *cutoff*\ , but applies to communication mode *multi* instead. Since in this case the communication cutoffs are determined per atom type, a type specifier is needed and cutoff for one or multiple types can be extended. Also ranges of types -using the usual asterisk notation can be given. +using the usual asterisk notation can be given. For granular pairstyles, +the default cutoff is set to the sum of the current maximum atomic radii +for each type. These are simulation scenarios in which it may be useful or even necessary to set a ghost cutoff > neighbor cutoff: diff --git a/src/npair_half_size_multi_newtoff.cpp b/src/npair_half_size_multi_newtoff.cpp index ec4ff805e2..813a642b96 100644 --- a/src/npair_half_size_multi_newtoff.cpp +++ b/src/npair_half_size_multi_newtoff.cpp @@ -11,9 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include #include "npair_half_size_multi_newtoff.h" -#include "neighbor.h" #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" @@ -27,6 +25,7 @@ using namespace LAMMPS_NS; NPairHalfSizeMultiNewtoff::NPairHalfSizeMultiNewtoff(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- + size particles binned neighbor list construction with partial Newton's 3rd law each owned atom i checks own bin and other bins in stencil multi-type stencil is itype dependent and is distance checked @@ -46,7 +45,6 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; - tagint *tag = atom->tag; tagint *molecule = atom->molecule; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; @@ -99,7 +97,7 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) + if (history && rsq < radsum*radsum) neighptr[n++] = j ^ mask_history; else neighptr[n++] = j; diff --git a/src/npair_half_size_multi_newtoff.h b/src/npair_half_size_multi_newtoff.h index 8abe5b5456..f255f9a17d 100644 --- a/src/npair_half_size_multi_newtoff.h +++ b/src/npair_half_size_multi_newtoff.h @@ -40,4 +40,7 @@ class NPairHalfSizeMultiNewtoff : public NPair { /* ERROR/WARNING messages: +E: Neighbor list overflow, boost neigh_modify one + +UNDOCUMENTED */ diff --git a/src/npair_half_size_multi_newton.cpp b/src/npair_half_size_multi_newton.cpp index 01aea66632..0d12924629 100644 --- a/src/npair_half_size_multi_newton.cpp +++ b/src/npair_half_size_multi_newton.cpp @@ -11,9 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include #include "npair_half_size_multi_newton.h" -#include "neighbor.h" #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" @@ -27,6 +25,7 @@ using namespace LAMMPS_NS; NPairHalfSizeMultiNewton::NPairHalfSizeMultiNewton(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- + size particles binned neighbor list construction with full Newton's 3rd law each owned atom i checks its own bin and other bins in Newton stencil multi-type stencil is itype dependent and is distance checked @@ -45,7 +44,6 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; - tagint *tag = atom->tag; tagint *molecule = atom->molecule; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; diff --git a/src/npair_half_size_multi_newton.h b/src/npair_half_size_multi_newton.h index 109fe78c2b..3e3d6f4180 100644 --- a/src/npair_half_size_multi_newton.h +++ b/src/npair_half_size_multi_newton.h @@ -40,4 +40,7 @@ class NPairHalfSizeMultiNewton : public NPair { /* ERROR/WARNING messages: +E: Neighbor list overflow, boost neigh_modify one + +UNDOCUMENTED */ diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp index f7123a958d..64b9659c8f 100644 --- a/src/npair_half_size_multi_newton_tri.cpp +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -45,7 +45,6 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; - tagint *tag = atom->tag; tagint *molecule = atom->molecule; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; @@ -72,7 +71,9 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) radi = radius[i]; - // loop over all atoms in bins in stencil + // loop over all atoms in bins, including self, in stencil + // skip if i,j neighbor cutoff is less than bin distance + // bins below self are excluded from stencil // pairs for atoms j "below" i are excluded // below = lower z or (equal z and lower y) or (equal zy and lower x) // (equal zyx and j <= i) diff --git a/src/npair_half_size_multi_newton_tri.h b/src/npair_half_size_multi_newton_tri.h index 15fbaba0f7..6afe8201a7 100644 --- a/src/npair_half_size_multi_newton_tri.h +++ b/src/npair_half_size_multi_newton_tri.h @@ -40,4 +40,7 @@ class NPairHalfSizeMultiNewtonTri : public NPair { /* ERROR/WARNING messages: +E: Neighbor list overflow, boost neigh_modify one + +UNDOCUMENTED */ From d0981db66a6e244fb176929831c4a3fe2c44b077 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Fri, 9 Oct 2020 13:51:35 -0600 Subject: [PATCH 05/25] Minor edits --- doc/src/comm_modify.rst | 2 +- src/npair_half_size_multi_newton_tri.cpp | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/doc/src/comm_modify.rst b/doc/src/comm_modify.rst index 68d40281c1..1ccf77e995 100644 --- a/doc/src/comm_modify.rst +++ b/doc/src/comm_modify.rst @@ -84,7 +84,7 @@ information is available, then also a heuristic based on that bond length is computed. It is used as communication cutoff, if there is no pair style present and no *comm_modify cutoff* command used. Otherwise a warning is printed, if this bond based estimate is larger than the -communication cutoff used. +communication cutoff used. The *cutoff/multi* option is equivalent to *cutoff*\ , but applies to communication mode *multi* instead. Since in this case the communication diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp index 64b9659c8f..adc4c080ec 100644 --- a/src/npair_half_size_multi_newton_tri.cpp +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -11,9 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include #include "npair_half_size_multi_newton_tri.h" -#include "neighbor.h" #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" From 0ae09c0f3b426ca7c82066359e0b79c4075e2523 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Wed, 14 Oct 2020 15:52:20 -0600 Subject: [PATCH 06/25] Adding OMP classes --- .../npair_half_size_multi_newtoff_omp.cpp | 128 +++++++++++++++ .../npair_half_size_multi_newtoff_omp.h | 44 +++++ .../npair_half_size_multi_newton_omp.cpp | 151 ++++++++++++++++++ .../npair_half_size_multi_newton_omp.h | 43 +++++ .../npair_half_size_multi_newton_tri_omp.cpp | 136 ++++++++++++++++ .../npair_half_size_multi_newton_tri_omp.h | 43 +++++ 6 files changed, 545 insertions(+) create mode 100644 src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp create mode 100644 src/USER-OMP/npair_half_size_multi_newtoff_omp.h create mode 100644 src/USER-OMP/npair_half_size_multi_newton_omp.cpp create mode 100644 src/USER-OMP/npair_half_size_multi_newton_omp.h create mode 100644 src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp create mode 100644 src/USER-OMP/npair_half_size_multi_newton_tri_omp.h diff --git a/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp b/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp new file mode 100644 index 0000000000..a72be3f982 --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp @@ -0,0 +1,128 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "omp_compat.h" +#include "npair_half_size_multi_newtoff_omp.h" +#include "npair_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewtoffOmp::NPairHalfSizeMultiNewtoffOmp(LAMMPS *lmp) : + NPair(lmp) {} + +/* ---------------------------------------------------------------------- + size particles + binned neighbor list construction with partial Newton's 3rd law + each owned atom i checks own bin and other bins in stencil + multi-type stencil is itype dependent and is distance checked + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int history = list->history; + const int mask_history = 3 << SBBITS; + + NPAIR_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(list) +#endif + NPAIR_OMP_SETUP(nlocal); + + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + // loop over each atom, storing neighbors + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + // each thread has its own page allocator + MyPage &ipage = list->ipage[tid]; + ipage.reset(); + + for (i = ifrom; i < ito; i++) { + + n = 0; + neighptr = ipage.vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over all atoms in other bins in stencil including self + // only store pair if i < j + // skip if i,j neighbor cutoff is less than bin distance + // stores own/own pairs only once + // stores own/ghost pairs on both procs + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage.vgot(n); + if (ipage.status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + NPAIR_OMP_CLOSE; + list->inum = nlocal; +} diff --git a/src/USER-OMP/npair_half_size_multi_newtoff_omp.h b/src/USER-OMP/npair_half_size_multi_newtoff_omp.h new file mode 100644 index 0000000000..dd883d7f3c --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newtoff_omp.h @@ -0,0 +1,44 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newtoff/omp, + NPairHalfSizeMultiNewtoffOmp, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTOFF | NP_OMP | + NP_ORTHO | NP_TRI) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_OMP_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_OMP_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewtoffOmp : public NPair { + public: + NPairHalfSizeMultiNewtoffOmp(class LAMMPS *); + ~NPairHalfSizeMultiNewtoffOmp() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/USER-OMP/npair_half_size_multi_newton_omp.cpp b/src/USER-OMP/npair_half_size_multi_newton_omp.cpp new file mode 100644 index 0000000000..751805dce1 --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newton_omp.cpp @@ -0,0 +1,151 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "omp_compat.h" +#include "npair_half_size_multi_newton_omp.h" +#include "npair_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewtonOmp::NPairHalfSizeMultiNewtonOmp(LAMMPS *lmp) : + NPair(lmp) {} + +/* ---------------------------------------------------------------------- + size particles + binned neighbor list construction with full Newton's 3rd law + each owned atom i checks its own bin and other bins in Newton stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int history = list->history; + const int mask_history = 3 << SBBITS; + + NPAIR_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(list) +#endif + NPAIR_OMP_SETUP(nlocal); + + int i,j,k,n,m,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + // each thread has its own page allocator + MyPage &ipage = list->ipage[tid]; + ipage.reset(); + + for (i = ifrom; i < ito; i++) { + + n = 0; + neighptr = ipage.vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over rest of atoms in i's bin, ghosts are at end of linked list + // if j is owned atom, store it, since j is beyond i in linked list + // if j is ghost, only store if j coords are "above and to the right" of i + + for (j = bins[i]; j >= 0; j = bins[j]) { + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + + // loop over all atoms in other bins in stencil, store every pair + // skip if i,j neighbor cutoff is less than bin distance + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage.vgot(n); + if (ipage.status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + NPAIR_OMP_CLOSE; + list->inum = nlocal; +} diff --git a/src/USER-OMP/npair_half_size_multi_newton_omp.h b/src/USER-OMP/npair_half_size_multi_newton_omp.h new file mode 100644 index 0000000000..15cc828629 --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newton_omp.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/bin/newton/omp, + NPairHalfSizeMultiNewtonOmp, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_OMP | NP_ORTHO) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_OMP_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_OMP_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewtonOmp : public NPair { + public: + NPairHalfSizeMultiNewtonOmp(class LAMMPS *); + ~NPairHalfSizeMultiNewtonOmp() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp new file mode 100644 index 0000000000..29832fba8d --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp @@ -0,0 +1,136 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "omp_compat.h" +#include "npair_half_size_multi_newton_tri_omp.h" +#include "npair_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfSizeMultiNewtonTriOmp::NPairHalfSizeMultiNewtonTriOmp(LAMMPS *lmp) : + NPair(lmp) {} + +/* ---------------------------------------------------------------------- + size particles + binned neighbor list construction with Newton's 3rd law for triclinic + each owned atom i checks its own bin and other bins in triclinic stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int history = list->history; + const int mask_history = 3 << SBBITS; + + NPAIR_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(list) +#endif + NPAIR_OMP_SETUP(nlocal); + + int i,j,k,m,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + double *cutsq,*distsq; + + // loop over each atom, storing neighbors + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + // each thread has its own page allocator + MyPage &ipage = list->ipage[tid]; + ipage.reset(); + + for (i = ifrom; i < ito; i++) { + + n = 0; + neighptr = ipage.vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over all atoms in bins, including self, in stencil + // skip if i,j neighbor cutoff is less than bin distance + // bins below self are excluded from stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and lower x) + // (equal zyx and j <= i) + // latter excludes self-self interaction but allows superposed atoms + + ibin = atom2bin[i]; + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp) { + if (x[j][0] < xtmp) continue; + if (x[j][0] == xtmp && j <= i) continue; + } + } + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage.vgot(n); + if (ipage.status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + NPAIR_OMP_CLOSE; + list->inum = nlocal; +} diff --git a/src/USER-OMP/npair_half_size_multi_newton_tri_omp.h b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.h new file mode 100644 index 0000000000..6e936c8da4 --- /dev/null +++ b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.h @@ -0,0 +1,43 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/size/multi/newton/tri/omp, + NPairHalfSizeMultiNewtonTriOmp, + NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_TRI | NP_OMP) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_OMP_H +#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_OMP_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeMultiNewtonTriOmp : public NPair { + public: + NPairHalfSizeMultiNewtonTriOmp(class LAMMPS *); + ~NPairHalfSizeMultiNewtonTriOmp() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ From d519f4fd4f0807c969c11e6c90ef52fcd8542459 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Wed, 14 Oct 2020 16:06:14 -0600 Subject: [PATCH 07/25] Missed reference to bin, minor uniform style changes --- src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp | 2 -- src/USER-OMP/npair_half_size_multi_newton_omp.cpp | 2 +- src/USER-OMP/npair_half_size_multi_newton_omp.h | 2 +- src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp | 2 -- 4 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp b/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp index a72be3f982..33722eb9d3 100644 --- a/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp +++ b/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp @@ -54,8 +54,6 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list) int *neighptr,*s; double *cutsq,*distsq; - // loop over each atom, storing neighbors - double **x = atom->x; double *radius = atom->radius; int *type = atom->type; diff --git a/src/USER-OMP/npair_half_size_multi_newton_omp.cpp b/src/USER-OMP/npair_half_size_multi_newton_omp.cpp index 751805dce1..0be87457e3 100644 --- a/src/USER-OMP/npair_half_size_multi_newton_omp.cpp +++ b/src/USER-OMP/npair_half_size_multi_newton_omp.cpp @@ -47,7 +47,7 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,m,itype,jtype,ibin,ns; + int i,j,k,m,n,itype,jtype,ibin,ns; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; diff --git a/src/USER-OMP/npair_half_size_multi_newton_omp.h b/src/USER-OMP/npair_half_size_multi_newton_omp.h index 15cc828629..03d712145f 100644 --- a/src/USER-OMP/npair_half_size_multi_newton_omp.h +++ b/src/USER-OMP/npair_half_size_multi_newton_omp.h @@ -13,7 +13,7 @@ #ifdef NPAIR_CLASS -NPairStyle(half/size/bin/newton/omp, +NPairStyle(half/size/multi/newton/omp, NPairHalfSizeMultiNewtonOmp, NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_OMP | NP_ORTHO) diff --git a/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp index 29832fba8d..ef26a5f2bd 100644 --- a/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp +++ b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp @@ -53,8 +53,6 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) int *neighptr,*s; double *cutsq,*distsq; - // loop over each atom, storing neighbors - double **x = atom->x; double *radius = atom->radius; int *type = atom->type; From 5a00d3c157186ff4e1163429dbfc541ec1d58f15 Mon Sep 17 00:00:00 2001 From: Julien Devemy Date: Tue, 1 Dec 2020 09:21:12 +0100 Subject: [PATCH 08/25] Fix invalid result when dihedral quadratic angle is > pi or < -pi --- src/USER-MISC/dihedral_quadratic.cpp | 2 ++ src/USER-OMP/dihedral_quadratic_omp.cpp | 4 ++++ 2 files changed, 6 insertions(+) diff --git a/src/USER-MISC/dihedral_quadratic.cpp b/src/USER-MISC/dihedral_quadratic.cpp index bfe13d13d1..ef5711215c 100644 --- a/src/USER-MISC/dihedral_quadratic.cpp +++ b/src/USER-MISC/dihedral_quadratic.cpp @@ -195,6 +195,8 @@ void DihedralQuadratic::compute(int eflag, int vflag) siinv = 1.0/si; double dphi = phi-phi0[type]; + if (dphi > MY_PI) dphi -= 2*MY_PI; + else if (dphi < -MY_PI) dphi += 2*MY_PI; p = k[type]*dphi; pd = - 2.0 * p * siinv; p = p * dphi; diff --git a/src/USER-OMP/dihedral_quadratic_omp.cpp b/src/USER-OMP/dihedral_quadratic_omp.cpp index 2e25258b13..d183500182 100644 --- a/src/USER-OMP/dihedral_quadratic_omp.cpp +++ b/src/USER-OMP/dihedral_quadratic_omp.cpp @@ -23,11 +23,13 @@ #include "neighbor.h" #include "force.h" #include "update.h" +#include "math_const.h" #include "error.h" #include "suffix.h" using namespace LAMMPS_NS; +using namespace MathConst; #define TOLERANCE 0.05 #define SMALL 0.001 @@ -218,6 +220,8 @@ void DihedralQuadraticOMP::eval(int nfrom, int nto, ThrData * const thr) siinv = 1.0/si; double dphi = phi-phi0[type]; + if (dphi > MY_PI) dphi -= 2*MY_PI; + else if (dphi < -MY_PI) dphi += 2*MY_PI; p = k[type]*dphi; pd = - 2.0 * p * siinv; p = p * dphi; From ba64e7c75c142ab77674f9d3b363f3d386541f6c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 07:20:53 -0500 Subject: [PATCH 09/25] simplify/improve multi-partition test --- unittest/c-library/test_library_mpi.cpp | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/unittest/c-library/test_library_mpi.cpp b/unittest/c-library/test_library_mpi.cpp index cf33ed13b5..7502da767a 100644 --- a/unittest/c-library/test_library_mpi.cpp +++ b/unittest/c-library/test_library_mpi.cpp @@ -173,7 +173,7 @@ TEST(MPI, multi_partition) MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &me); - const char *args[] = {"LAMMPS_test", "-log", "none", "-partition", "2x2", + const char *args[] = {"LAMMPS_test", "-log", "none", "-partition", "4x1", "-echo", "screen", "-nocite", "-in", "none"}; char **argv = (char **)args; int argc = sizeof(args) / sizeof(char *); @@ -183,19 +183,15 @@ TEST(MPI, multi_partition) lammps_command(lmp, "atom_style atomic"); lammps_command(lmp, "region box block 0 2 0 2 0 2"); lammps_command(lmp, "create_box 1 box"); - lammps_command(lmp, "variable partition universe 1 2"); + lammps_command(lmp, "variable partition universe 1 2 3 4"); EXPECT_EQ(lammps_extract_setting(lmp, "universe_size"), nprocs); EXPECT_EQ(lammps_extract_setting(lmp, "universe_rank"), me); - EXPECT_EQ(lammps_extract_setting(lmp, "world_size"), nprocs / 2); - EXPECT_EQ(lammps_extract_setting(lmp, "world_rank"), me % 2); + EXPECT_EQ(lammps_extract_setting(lmp, "world_size"), 1); + EXPECT_EQ(lammps_extract_setting(lmp, "world_rank"), 0); char *part_id = (char *)lammps_extract_variable(lmp, "partition", nullptr); - if (me < 2) { - ASSERT_THAT(part_id, StrEq("1")); - } else { - ASSERT_THAT(part_id, StrEq("2")); - } + ASSERT_THAT(part_id, StrEq(std::to_string(me+1))); lammps_close(lmp); }; From 181b18beeb675a20237767bca3547cd3e1bbc53f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 07:21:17 -0500 Subject: [PATCH 10/25] tweak epsilon for better compiler compatibility --- unittest/force-styles/tests/fix-timestep-oneway.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unittest/force-styles/tests/fix-timestep-oneway.yaml b/unittest/force-styles/tests/fix-timestep-oneway.yaml index 4fcef41f5b..01131a3803 100644 --- a/unittest/force-styles/tests/fix-timestep-oneway.yaml +++ b/unittest/force-styles/tests/fix-timestep-oneway.yaml @@ -1,7 +1,7 @@ --- lammps_version: 24 Aug 2020 date_generated: Tue Sep 15 09:44:41 202 -epsilon: 2e-14 +epsilon: 4e-14 prerequisites: ! | atom full fix oneway From e3c3106795a3f8fc4906fcfb5f56e0992caa9cce Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 09:20:38 -0500 Subject: [PATCH 11/25] fix incorrect use of MPI_Gather() --- src/fix_temp_csld.cpp | 2 +- src/fix_temp_csvr.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fix_temp_csld.cpp b/src/fix_temp_csld.cpp index f01de7c2fa..48c18c4e4f 100644 --- a/src/fix_temp_csld.cpp +++ b/src/fix_temp_csld.cpp @@ -314,7 +314,7 @@ void FixTempCSLD::write_restart(FILE *fp) } double state[103]; random->get_state(state); - MPI_Gather(state,103,MPI_DOUBLE,list+2,103*comm->nprocs,MPI_DOUBLE,0,world); + MPI_Gather(state,103,MPI_DOUBLE,list+2,103,MPI_DOUBLE,0,world); if (comm->me == 0) { int size = nsize * sizeof(double); diff --git a/src/fix_temp_csvr.cpp b/src/fix_temp_csvr.cpp index 0130631172..e3fdb7fc80 100644 --- a/src/fix_temp_csvr.cpp +++ b/src/fix_temp_csvr.cpp @@ -347,7 +347,7 @@ void FixTempCSVR::write_restart(FILE *fp) } double state[103]; random->get_state(state); - MPI_Gather(state,103,MPI_DOUBLE,list+2,103*comm->nprocs,MPI_DOUBLE,0,world); + MPI_Gather(state,103,MPI_DOUBLE,list+2,103,MPI_DOUBLE,0,world); if (comm->me == 0) { int size = nsize * sizeof(double); From 5ddeb45c0a27ae7e263b8323a5357a1dbc5ca24d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 09:25:54 -0500 Subject: [PATCH 12/25] whitespace --- doc/src/Build_extras.rst | 8 ++++---- doc/src/fix_temp_csvr.rst | 18 +++++++++--------- src/KIM/pair_kim.cpp | 2 +- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/doc/src/Build_extras.rst b/doc/src/Build_extras.rst index 6847b8c030..4dd17d71f8 100644 --- a/doc/src/Build_extras.rst +++ b/doc/src/Build_extras.rst @@ -352,7 +352,7 @@ KIM Extra unit tests (CMake only) During development, testing, or debugging, if :doc:`unit testing ` is enabled in LAMMPS, one can also enable extra tests on :doc:`KIM commands ` by setting the -``KIM_EXTRA_UNITTESTS`` to *yes* (or *on*). +``KIM_EXTRA_UNITTESTS`` to *yes* (or *on*). Enabling the extra unit tests have some requirements, @@ -367,9 +367,9 @@ Enabling the extra unit tests have some requirements, *conda-forge* channel as ``conda install kim-property`` if LAMMPS is built in Conda. More detailed information is available at: `kim-property installation `_. -* It is also necessary to install - ``EAM_Dynamo_Mendelev_2007_Zr__MO_848899341753_000``, and - ``EAM_Dynamo_ErcolessiAdams_1994_Al__MO_123629422045_005`` KIM models. +* It is also necessary to install + ``EAM_Dynamo_Mendelev_2007_Zr__MO_848899341753_000``, and + ``EAM_Dynamo_ErcolessiAdams_1994_Al__MO_123629422045_005`` KIM models. See `Obtaining KIM Models `_ to learn how to install a pre-build binary of the OpenKIM Repository of Models or see diff --git a/doc/src/fix_temp_csvr.rst b/doc/src/fix_temp_csvr.rst index 5e6074bc66..440a264ed5 100644 --- a/doc/src/fix_temp_csvr.rst +++ b/doc/src/fix_temp_csvr.rst @@ -124,16 +124,16 @@ temperature is calculated taking the bias into account, bias is removed from each atom, thermostatting is performed on the remaining thermal degrees of freedom, and the bias is added back in. -An important feature of these thermostats is that they have an +An important feature of these thermostats is that they have an associated effective energy that is a constant of motion. -The effective energy is the total energy (kinetic + potential) plus -the accumulated kinetic energy changes due to the thermostat. The -latter quantity is the global scalar computed by these fixes. This -feature is useful to check the integration of the equations of motion -against discretization errors. In other words, the conservation of -the effective energy can be used to choose an appropriate integration -:doc:`timestep `. This is similar to the usual paradigm of -checking the conservation of the total energy in the microcanonical +The effective energy is the total energy (kinetic + potential) plus +the accumulated kinetic energy changes due to the thermostat. The +latter quantity is the global scalar computed by these fixes. This +feature is useful to check the integration of the equations of motion +against discretization errors. In other words, the conservation of +the effective energy can be used to choose an appropriate integration +:doc:`timestep `. This is similar to the usual paradigm of +checking the conservation of the total energy in the microcanonical ensemble. diff --git a/src/KIM/pair_kim.cpp b/src/KIM/pair_kim.cpp index ad671de5b3..ee6f5f98a5 100644 --- a/src/KIM/pair_kim.cpp +++ b/src/KIM/pair_kim.cpp @@ -941,7 +941,7 @@ void PairKIM::set_argument_pointers() // Set KIM pointer appropriately for particleVirial if (KIM_SupportStatus_Equal(kim_model_support_for_particleVirial, - KIM_SUPPORT_STATUS_required) + KIM_SUPPORT_STATUS_required) && (vflag_atom == 0)) { // reallocate per-atom virial array if necessary if (atom->nmax > maxvatom) { From 527ffa79dc11f41cbe3a8a8cd9d369b9103e520a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 09:34:56 -0500 Subject: [PATCH 13/25] refactor for clarity --- src/fix_temp_csld.cpp | 7 ++++--- src/fix_temp_csvr.cpp | 7 ++++--- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/fix_temp_csld.cpp b/src/fix_temp_csld.cpp index 48c18c4e4f..3b522c185f 100644 --- a/src/fix_temp_csld.cpp +++ b/src/fix_temp_csld.cpp @@ -305,16 +305,17 @@ double FixTempCSLD::compute_scalar() void FixTempCSLD::write_restart(FILE *fp) { - int nsize = (98+2+3)*comm->nprocs+2; // pRNG state per proc + nprocs + energy + const int PRNGSIZE = 98+2+3; + int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy double *list = nullptr; if (comm->me == 0) { list = new double[nsize]; list[0] = energy; list[1] = comm->nprocs; } - double state[103]; + double state[PRNGSIZE]; random->get_state(state); - MPI_Gather(state,103,MPI_DOUBLE,list+2,103,MPI_DOUBLE,0,world); + MPI_Gather(state,PRNGSIZE,MPI_DOUBLE,list+2,PRNGSIZE,MPI_DOUBLE,0,world); if (comm->me == 0) { int size = nsize * sizeof(double); diff --git a/src/fix_temp_csvr.cpp b/src/fix_temp_csvr.cpp index e3fdb7fc80..674c436e9d 100644 --- a/src/fix_temp_csvr.cpp +++ b/src/fix_temp_csvr.cpp @@ -338,16 +338,17 @@ double FixTempCSVR::compute_scalar() void FixTempCSVR::write_restart(FILE *fp) { - int nsize = (98+2+3)*comm->nprocs+2; // pRNG state per proc + nprocs + energy + const int PRNGSIZE = 98+2+3; + int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy double *list = nullptr; if (comm->me == 0) { list = new double[nsize]; list[0] = energy; list[1] = comm->nprocs; } - double state[103]; + double state[PRNGSIZE]; random->get_state(state); - MPI_Gather(state,103,MPI_DOUBLE,list+2,103,MPI_DOUBLE,0,world); + MPI_Gather(state,PRNGSIZE,MPI_DOUBLE,list+2,PRNGSIZE,MPI_DOUBLE,0,world); if (comm->me == 0) { int size = nsize * sizeof(double); From 47cdafe651a70fe228a2159fd41af226dd0d5b91 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 12:06:44 -0500 Subject: [PATCH 14/25] reduce size of executable/library by replacing static const int with enum --- src/fix.h | 51 +++++++++++++++-------------- src/neighbor.h | 88 +++++++++++++++++++++++++++----------------------- 2 files changed, 73 insertions(+), 66 deletions(-) diff --git a/src/fix.h b/src/fix.h index 75c644e87f..ec951e3bb2 100644 --- a/src/fix.h +++ b/src/fix.h @@ -251,31 +251,32 @@ class Fix : protected Pointers { }; namespace FixConst { - static const int INITIAL_INTEGRATE = 1<<0; - static const int POST_INTEGRATE = 1<<1; - static const int PRE_EXCHANGE = 1<<2; - static const int PRE_NEIGHBOR = 1<<3; - static const int POST_NEIGHBOR = 1<<4; - static const int PRE_FORCE = 1<<5; - static const int PRE_REVERSE = 1<<6; - static const int POST_FORCE = 1<<7; - static const int FINAL_INTEGRATE = 1<<8; - static const int END_OF_STEP = 1<<9; - static const int POST_RUN = 1<<10; - static const int THERMO_ENERGY = 1<<11; - static const int INITIAL_INTEGRATE_RESPA = 1<<12; - static const int POST_INTEGRATE_RESPA = 1<<13; - static const int PRE_FORCE_RESPA = 1<<14; - static const int POST_FORCE_RESPA = 1<<15; - static const int FINAL_INTEGRATE_RESPA = 1<<16; - static const int MIN_PRE_EXCHANGE = 1<<17; - static const int MIN_PRE_NEIGHBOR = 1<<18; - static const int MIN_POST_NEIGHBOR = 1<<19; - static const int MIN_PRE_FORCE = 1<<20; - static const int MIN_PRE_REVERSE = 1<<21; - static const int MIN_POST_FORCE = 1<<22; - static const int MIN_ENERGY = 1<<23; - static const int FIX_CONST_LAST = 1<<24; + enum { + INITIAL_INTEGRATE = 1<<0, + POST_INTEGRATE = 1<<1, + PRE_EXCHANGE = 1<<2, + PRE_NEIGHBOR = 1<<3, + POST_NEIGHBOR = 1<<4, + PRE_FORCE = 1<<5, + PRE_REVERSE = 1<<6, + POST_FORCE = 1<<7, + FINAL_INTEGRATE = 1<<8, + END_OF_STEP = 1<<9, + POST_RUN = 1<<10, + THERMO_ENERGY = 1<<11, + INITIAL_INTEGRATE_RESPA = 1<<12, + POST_INTEGRATE_RESPA = 1<<13, + PRE_FORCE_RESPA = 1<<14, + POST_FORCE_RESPA = 1<<15, + FINAL_INTEGRATE_RESPA = 1<<16, + MIN_PRE_EXCHANGE = 1<<17, + MIN_PRE_NEIGHBOR = 1<<18, + MIN_POST_NEIGHBOR = 1<<19, + MIN_PRE_FORCE = 1<<20, + MIN_PRE_REVERSE = 1<<21, + MIN_POST_FORCE = 1<<22, + MIN_ENERGY = 1<<23 + }; } } diff --git a/src/neighbor.h b/src/neighbor.h index 9ee2af9c75..b9b40bcf1a 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -235,49 +235,55 @@ class Neighbor : protected Pointers { }; namespace NeighConst { - static const int NB_INTEL = 1<<0; - static const int NB_KOKKOS_DEVICE = 1<<1; - static const int NB_KOKKOS_HOST = 1<<2; - static const int NB_SSA = 1<<3; + enum { + NB_INTEL = 1<<0, + NB_KOKKOS_DEVICE = 1<<1, + NB_KOKKOS_HOST = 1<<2, + NB_SSA = 1<<3 + }; - static const int NS_BIN = 1<<0; - static const int NS_MULTI = 1<<1; - static const int NS_HALF = 1<<2; - static const int NS_FULL = 1<<3; - static const int NS_2D = 1<<4; - static const int NS_3D = 1<<5; - static const int NS_NEWTON = 1<<6; - static const int NS_NEWTOFF = 1<<7; - static const int NS_ORTHO = 1<<8; - static const int NS_TRI = 1<<9; - static const int NS_GHOST = 1<<10; - static const int NS_SSA = 1<<11; + enum { + NS_BIN = 1<<0, + NS_MULTI = 1<<1, + NS_HALF = 1<<2, + NS_FULL = 1<<3, + NS_2D = 1<<4, + NS_3D = 1<<5, + NS_NEWTON = 1<<6, + NS_NEWTOFF = 1<<7, + NS_ORTHO = 1<<8, + NS_TRI = 1<<9, + NS_GHOST = 1<<10, + NS_SSA = 1<<11 + }; - static const int NP_NSQ = 1<<0; - static const int NP_BIN = 1<<1; - static const int NP_MULTI = 1<<2; - static const int NP_HALF = 1<<3; - static const int NP_FULL = 1<<4; - static const int NP_ORTHO = 1<<5; - static const int NP_TRI = 1<<6; - static const int NP_ATOMONLY = 1<<7; - static const int NP_MOLONLY = 1<<8; - static const int NP_NEWTON = 1<<9; - static const int NP_NEWTOFF = 1<<10; - static const int NP_GHOST = 1<<11; - static const int NP_SIZE = 1<<12; - static const int NP_ONESIDE = 1<<13; - static const int NP_RESPA = 1<<14; - static const int NP_BOND = 1<<15; - static const int NP_OMP = 1<<16; - static const int NP_INTEL = 1<<17; - static const int NP_KOKKOS_DEVICE = 1<<18; - static const int NP_KOKKOS_HOST = 1<<19; - static const int NP_SSA = 1<<20; - static const int NP_COPY = 1<<21; - static const int NP_SKIP = 1<<22; - static const int NP_HALF_FULL = 1<<23; - static const int NP_OFF2ON = 1<<24; + enum { + NP_NSQ = 1<<0, + NP_BIN = 1<<1, + NP_MULTI = 1<<2, + NP_HALF = 1<<3, + NP_FULL = 1<<4, + NP_ORTHO = 1<<5, + NP_TRI = 1<<6, + NP_ATOMONLY = 1<<7, + NP_MOLONLY = 1<<8, + NP_NEWTON = 1<<9, + NP_NEWTOFF = 1<<10, + NP_GHOST = 1<<11, + NP_SIZE = 1<<12, + NP_ONESIDE = 1<<13, + NP_RESPA = 1<<14, + NP_BOND = 1<<15, + NP_OMP = 1<<16, + NP_INTEL = 1<<17, + NP_KOKKOS_DEVICE = 1<<18, + NP_KOKKOS_HOST = 1<<19, + NP_SSA = 1<<20, + NP_COPY = 1<<21, + NP_SKIP = 1<<22, + NP_HALF_FULL = 1<<23, + NP_OFF2ON = 1<<24 + }; } } From 72f68f3d56d01cf888e268703057491276d8359a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 12:47:04 -0500 Subject: [PATCH 15/25] consolidate multiple factorial() function definitions into MathSpecial::factorial() --- .../compute_orientorder_atom_kokkos.cpp | 2 + src/SNAP/sna.cpp | 192 +---------------- src/SNAP/sna.h | 4 - src/USER-SMTBQ/pair_smtbq.cpp | 11 - src/compute_orientorder_atom.cpp | 191 +---------------- src/compute_orientorder_atom.h | 3 - src/math_special.cpp | 194 ++++++++++++++++++ src/math_special.h | 4 + 8 files changed, 205 insertions(+), 396 deletions(-) diff --git a/src/KOKKOS/compute_orientorder_atom_kokkos.cpp b/src/KOKKOS/compute_orientorder_atom_kokkos.cpp index 1323e589d3..3703769e90 100644 --- a/src/KOKKOS/compute_orientorder_atom_kokkos.cpp +++ b/src/KOKKOS/compute_orientorder_atom_kokkos.cpp @@ -21,6 +21,7 @@ #include "atom_masks.h" #include "kokkos.h" #include "math_const.h" +#include "math_special.h" #include "memory_kokkos.h" #include "neigh_list.h" #include "neigh_request.h" @@ -32,6 +33,7 @@ using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; using namespace std; #ifdef DBL_EPSILON diff --git a/src/SNAP/sna.cpp b/src/SNAP/sna.cpp index 6d4197b59b..3c7a40c0dc 100644 --- a/src/SNAP/sna.cpp +++ b/src/SNAP/sna.cpp @@ -18,6 +18,7 @@ #include "sna.h" #include #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "error.h" #include "comm.h" @@ -25,6 +26,7 @@ using namespace std; using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; /* ---------------------------------------------------------------------- @@ -1363,196 +1365,6 @@ void SNA::destroy_twojmax_arrays() } -/* ---------------------------------------------------------------------- - factorial n, wrapper for precomputed table -------------------------------------------------------------------------- */ - -double SNA::factorial(int n) -{ - if (n < 0 || n > nmaxfactorial) { - char str[128]; - sprintf(str, "Invalid argument to factorial %d", n); - error->all(FLERR, str); - } - - return nfac_table[n]; -} - -/* ---------------------------------------------------------------------- - factorial n table, size SNA::nmaxfactorial+1 -------------------------------------------------------------------------- */ - -const double SNA::nfac_table[] = { - 1, - 1, - 2, - 6, - 24, - 120, - 720, - 5040, - 40320, - 362880, - 3628800, - 39916800, - 479001600, - 6227020800, - 87178291200, - 1307674368000, - 20922789888000, - 355687428096000, - 6.402373705728e+15, - 1.21645100408832e+17, - 2.43290200817664e+18, - 5.10909421717094e+19, - 1.12400072777761e+21, - 2.5852016738885e+22, - 6.20448401733239e+23, - 1.5511210043331e+25, - 4.03291461126606e+26, - 1.08888694504184e+28, - 3.04888344611714e+29, - 8.8417619937397e+30, - 2.65252859812191e+32, - 8.22283865417792e+33, - 2.63130836933694e+35, - 8.68331761881189e+36, - 2.95232799039604e+38, - 1.03331479663861e+40, - 3.71993326789901e+41, - 1.37637530912263e+43, - 5.23022617466601e+44, - 2.03978820811974e+46, - 8.15915283247898e+47, - 3.34525266131638e+49, - 1.40500611775288e+51, - 6.04152630633738e+52, - 2.65827157478845e+54, - 1.1962222086548e+56, - 5.50262215981209e+57, - 2.58623241511168e+59, - 1.24139155925361e+61, - 6.08281864034268e+62, - 3.04140932017134e+64, - 1.55111875328738e+66, - 8.06581751709439e+67, - 4.27488328406003e+69, - 2.30843697339241e+71, - 1.26964033536583e+73, - 7.10998587804863e+74, - 4.05269195048772e+76, - 2.35056133128288e+78, - 1.3868311854569e+80, - 8.32098711274139e+81, - 5.07580213877225e+83, - 3.14699732603879e+85, - 1.98260831540444e+87, - 1.26886932185884e+89, - 8.24765059208247e+90, - 5.44344939077443e+92, - 3.64711109181887e+94, - 2.48003554243683e+96, - 1.71122452428141e+98, - 1.19785716699699e+100, - 8.50478588567862e+101, - 6.12344583768861e+103, - 4.47011546151268e+105, - 3.30788544151939e+107, - 2.48091408113954e+109, - 1.88549470166605e+111, - 1.45183092028286e+113, - 1.13242811782063e+115, - 8.94618213078297e+116, - 7.15694570462638e+118, - 5.79712602074737e+120, - 4.75364333701284e+122, - 3.94552396972066e+124, - 3.31424013456535e+126, - 2.81710411438055e+128, - 2.42270953836727e+130, - 2.10775729837953e+132, - 1.85482642257398e+134, - 1.65079551609085e+136, - 1.48571596448176e+138, - 1.3520015276784e+140, - 1.24384140546413e+142, - 1.15677250708164e+144, - 1.08736615665674e+146, - 1.03299784882391e+148, - 9.91677934870949e+149, - 9.61927596824821e+151, - 9.42689044888324e+153, - 9.33262154439441e+155, - 9.33262154439441e+157, - 9.42594775983835e+159, - 9.61446671503512e+161, - 9.90290071648618e+163, - 1.02990167451456e+166, - 1.08139675824029e+168, - 1.14628056373471e+170, - 1.22652020319614e+172, - 1.32464181945183e+174, - 1.44385958320249e+176, - 1.58824554152274e+178, - 1.76295255109024e+180, - 1.97450685722107e+182, - 2.23119274865981e+184, - 2.54355973347219e+186, - 2.92509369349301e+188, - 3.3931086844519e+190, - 3.96993716080872e+192, - 4.68452584975429e+194, - 5.5745857612076e+196, - 6.68950291344912e+198, - 8.09429852527344e+200, - 9.8750442008336e+202, - 1.21463043670253e+205, - 1.50614174151114e+207, - 1.88267717688893e+209, - 2.37217324288005e+211, - 3.01266001845766e+213, - 3.8562048236258e+215, - 4.97450422247729e+217, - 6.46685548922047e+219, - 8.47158069087882e+221, - 1.118248651196e+224, - 1.48727070609069e+226, - 1.99294274616152e+228, - 2.69047270731805e+230, - 3.65904288195255e+232, - 5.01288874827499e+234, - 6.91778647261949e+236, - 9.61572319694109e+238, - 1.34620124757175e+241, - 1.89814375907617e+243, - 2.69536413788816e+245, - 3.85437071718007e+247, - 5.5502938327393e+249, - 8.04792605747199e+251, - 1.17499720439091e+254, - 1.72724589045464e+256, - 2.55632391787286e+258, - 3.80892263763057e+260, - 5.71338395644585e+262, - 8.62720977423323e+264, - 1.31133588568345e+267, - 2.00634390509568e+269, - 3.08976961384735e+271, - 4.78914290146339e+273, - 7.47106292628289e+275, - 1.17295687942641e+278, - 1.85327186949373e+280, - 2.94670227249504e+282, - 4.71472363599206e+284, - 7.59070505394721e+286, - 1.22969421873945e+289, - 2.0044015765453e+291, - 3.28721858553429e+293, - 5.42391066613159e+295, - 9.00369170577843e+297, - 1.503616514865e+300, // nmaxfactorial = 167 -}; - /* ---------------------------------------------------------------------- the function delta given by VMK Eq. 8.2(1) ------------------------------------------------------------------------- */ diff --git a/src/SNAP/sna.h b/src/SNAP/sna.h index 30c4aa6bbc..746b55fb70 100644 --- a/src/SNAP/sna.h +++ b/src/SNAP/sna.h @@ -100,10 +100,6 @@ private: double** dulist_r, ** dulist_i; int elem_duarray; // element of j in derivative - static const int nmaxfactorial = 167; - static const double nfac_table[]; - double factorial(int); - void create_twojmax_arrays(); void destroy_twojmax_arrays(); void init_clebsch_gordan(); diff --git a/src/USER-SMTBQ/pair_smtbq.cpp b/src/USER-SMTBQ/pair_smtbq.cpp index 94647cb552..aedf90da07 100644 --- a/src/USER-SMTBQ/pair_smtbq.cpp +++ b/src/USER-SMTBQ/pair_smtbq.cpp @@ -72,17 +72,6 @@ using namespace MathSpecial; #define PGDELTA 1 #define MAXNEIGH 24 -/* ------------------------------------------------------------------------------------ - - Calculates the factorial of an integer n via recursion. - - ------------------------------------------------------------------------------------ */ -static double factorial(int n) -{ - if (n <= 1) return 1.0; - else return static_cast(n)*factorial(n-1); -} - /* ---------------------------------------------------------------------- */ PairSMTBQ::PairSMTBQ(LAMMPS *lmp) : Pair(lmp) diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp index db8a2f404b..b334654e82 100644 --- a/src/compute_orientorder_atom.cpp +++ b/src/compute_orientorder_atom.cpp @@ -23,6 +23,7 @@ #include "error.h" #include "force.h" #include "math_const.h" +#include "math_special.h" #include "memory.h" #include "modify.h" #include "neigh_list.h" @@ -36,6 +37,8 @@ using namespace LAMMPS_NS; using namespace MathConst; +using namespace MathSpecial; + #ifdef DBL_EPSILON #define MY_EPSILON (10.0*DBL_EPSILON) @@ -708,191 +711,3 @@ void ComputeOrientOrderAtom::init_clebsch_gordan() } } } - -/* ---------------------------------------------------------------------- - factorial n, wrapper for precomputed table -------------------------------------------------------------------------- */ - -double ComputeOrientOrderAtom::factorial(int n) -{ - if (n < 0 || n > nmaxfactorial) - error->all(FLERR,fmt::format("Invalid argument to factorial {}", n)); - - return nfac_table[n]; -} - -/* ---------------------------------------------------------------------- - factorial n table, size SNA::nmaxfactorial+1 -------------------------------------------------------------------------- */ - -const double ComputeOrientOrderAtom::nfac_table[] = { - 1, - 1, - 2, - 6, - 24, - 120, - 720, - 5040, - 40320, - 362880, - 3628800, - 39916800, - 479001600, - 6227020800, - 87178291200, - 1307674368000, - 20922789888000, - 355687428096000, - 6.402373705728e+15, - 1.21645100408832e+17, - 2.43290200817664e+18, - 5.10909421717094e+19, - 1.12400072777761e+21, - 2.5852016738885e+22, - 6.20448401733239e+23, - 1.5511210043331e+25, - 4.03291461126606e+26, - 1.08888694504184e+28, - 3.04888344611714e+29, - 8.8417619937397e+30, - 2.65252859812191e+32, - 8.22283865417792e+33, - 2.63130836933694e+35, - 8.68331761881189e+36, - 2.95232799039604e+38, - 1.03331479663861e+40, - 3.71993326789901e+41, - 1.37637530912263e+43, - 5.23022617466601e+44, - 2.03978820811974e+46, - 8.15915283247898e+47, - 3.34525266131638e+49, - 1.40500611775288e+51, - 6.04152630633738e+52, - 2.65827157478845e+54, - 1.1962222086548e+56, - 5.50262215981209e+57, - 2.58623241511168e+59, - 1.24139155925361e+61, - 6.08281864034268e+62, - 3.04140932017134e+64, - 1.55111875328738e+66, - 8.06581751709439e+67, - 4.27488328406003e+69, - 2.30843697339241e+71, - 1.26964033536583e+73, - 7.10998587804863e+74, - 4.05269195048772e+76, - 2.35056133128288e+78, - 1.3868311854569e+80, - 8.32098711274139e+81, - 5.07580213877225e+83, - 3.14699732603879e+85, - 1.98260831540444e+87, - 1.26886932185884e+89, - 8.24765059208247e+90, - 5.44344939077443e+92, - 3.64711109181887e+94, - 2.48003554243683e+96, - 1.71122452428141e+98, - 1.19785716699699e+100, - 8.50478588567862e+101, - 6.12344583768861e+103, - 4.47011546151268e+105, - 3.30788544151939e+107, - 2.48091408113954e+109, - 1.88549470166605e+111, - 1.45183092028286e+113, - 1.13242811782063e+115, - 8.94618213078297e+116, - 7.15694570462638e+118, - 5.79712602074737e+120, - 4.75364333701284e+122, - 3.94552396972066e+124, - 3.31424013456535e+126, - 2.81710411438055e+128, - 2.42270953836727e+130, - 2.10775729837953e+132, - 1.85482642257398e+134, - 1.65079551609085e+136, - 1.48571596448176e+138, - 1.3520015276784e+140, - 1.24384140546413e+142, - 1.15677250708164e+144, - 1.08736615665674e+146, - 1.03299784882391e+148, - 9.91677934870949e+149, - 9.61927596824821e+151, - 9.42689044888324e+153, - 9.33262154439441e+155, - 9.33262154439441e+157, - 9.42594775983835e+159, - 9.61446671503512e+161, - 9.90290071648618e+163, - 1.02990167451456e+166, - 1.08139675824029e+168, - 1.14628056373471e+170, - 1.22652020319614e+172, - 1.32464181945183e+174, - 1.44385958320249e+176, - 1.58824554152274e+178, - 1.76295255109024e+180, - 1.97450685722107e+182, - 2.23119274865981e+184, - 2.54355973347219e+186, - 2.92509369349301e+188, - 3.3931086844519e+190, - 3.96993716080872e+192, - 4.68452584975429e+194, - 5.5745857612076e+196, - 6.68950291344912e+198, - 8.09429852527344e+200, - 9.8750442008336e+202, - 1.21463043670253e+205, - 1.50614174151114e+207, - 1.88267717688893e+209, - 2.37217324288005e+211, - 3.01266001845766e+213, - 3.8562048236258e+215, - 4.97450422247729e+217, - 6.46685548922047e+219, - 8.47158069087882e+221, - 1.118248651196e+224, - 1.48727070609069e+226, - 1.99294274616152e+228, - 2.69047270731805e+230, - 3.65904288195255e+232, - 5.01288874827499e+234, - 6.91778647261949e+236, - 9.61572319694109e+238, - 1.34620124757175e+241, - 1.89814375907617e+243, - 2.69536413788816e+245, - 3.85437071718007e+247, - 5.5502938327393e+249, - 8.04792605747199e+251, - 1.17499720439091e+254, - 1.72724589045464e+256, - 2.55632391787286e+258, - 3.80892263763057e+260, - 5.71338395644585e+262, - 8.62720977423323e+264, - 1.31133588568345e+267, - 2.00634390509568e+269, - 3.08976961384735e+271, - 4.78914290146339e+273, - 7.47106292628289e+275, - 1.17295687942641e+278, - 1.85327186949373e+280, - 2.94670227249504e+282, - 4.71472363599206e+284, - 7.59070505394721e+286, - 1.22969421873945e+289, - 2.0044015765453e+291, - 3.28721858553429e+293, - 5.42391066613159e+295, - 9.00369170577843e+297, - 1.503616514865e+300, // nmaxfactorial = 167 -}; - diff --git a/src/compute_orientorder_atom.h b/src/compute_orientorder_atom.h index 4e5084bfcd..4fc122c5c8 100644 --- a/src/compute_orientorder_atom.h +++ b/src/compute_orientorder_atom.h @@ -56,9 +56,6 @@ class ComputeOrientOrderAtom : public Compute { double polar_prefactor(int, int, double); double associated_legendre(int, int, double); - static const int nmaxfactorial = 167; - static const double nfac_table[]; - double factorial(int); virtual void init_clebsch_gordan(); double *cglist; // Clebsch-Gordan coeffs int idxcg_max; diff --git a/src/math_special.cpp b/src/math_special.cpp index 5ec0fc3f61..4ed91ba51f 100644 --- a/src/math_special.cpp +++ b/src/math_special.cpp @@ -1,9 +1,203 @@ #include "math_special.h" + #include #include // IWYU pragma: keep +#include + +#include "error.h" using namespace LAMMPS_NS; +static constexpr int nmaxfactorial = 167; + +/* ---------------------------------------------------------------------- + factorial n table, size nmaxfactorial+1 +------------------------------------------------------------------------- */ + +static const double nfac_table[] = { + 1, + 1, + 2, + 6, + 24, + 120, + 720, + 5040, + 40320, + 362880, + 3628800, + 39916800, + 479001600, + 6227020800, + 87178291200, + 1307674368000, + 20922789888000, + 355687428096000, + 6.402373705728e+15, + 1.21645100408832e+17, + 2.43290200817664e+18, + 5.10909421717094e+19, + 1.12400072777761e+21, + 2.5852016738885e+22, + 6.20448401733239e+23, + 1.5511210043331e+25, + 4.03291461126606e+26, + 1.08888694504184e+28, + 3.04888344611714e+29, + 8.8417619937397e+30, + 2.65252859812191e+32, + 8.22283865417792e+33, + 2.63130836933694e+35, + 8.68331761881189e+36, + 2.95232799039604e+38, + 1.03331479663861e+40, + 3.71993326789901e+41, + 1.37637530912263e+43, + 5.23022617466601e+44, + 2.03978820811974e+46, + 8.15915283247898e+47, + 3.34525266131638e+49, + 1.40500611775288e+51, + 6.04152630633738e+52, + 2.65827157478845e+54, + 1.1962222086548e+56, + 5.50262215981209e+57, + 2.58623241511168e+59, + 1.24139155925361e+61, + 6.08281864034268e+62, + 3.04140932017134e+64, + 1.55111875328738e+66, + 8.06581751709439e+67, + 4.27488328406003e+69, + 2.30843697339241e+71, + 1.26964033536583e+73, + 7.10998587804863e+74, + 4.05269195048772e+76, + 2.35056133128288e+78, + 1.3868311854569e+80, + 8.32098711274139e+81, + 5.07580213877225e+83, + 3.14699732603879e+85, + 1.98260831540444e+87, + 1.26886932185884e+89, + 8.24765059208247e+90, + 5.44344939077443e+92, + 3.64711109181887e+94, + 2.48003554243683e+96, + 1.71122452428141e+98, + 1.19785716699699e+100, + 8.50478588567862e+101, + 6.12344583768861e+103, + 4.47011546151268e+105, + 3.30788544151939e+107, + 2.48091408113954e+109, + 1.88549470166605e+111, + 1.45183092028286e+113, + 1.13242811782063e+115, + 8.94618213078297e+116, + 7.15694570462638e+118, + 5.79712602074737e+120, + 4.75364333701284e+122, + 3.94552396972066e+124, + 3.31424013456535e+126, + 2.81710411438055e+128, + 2.42270953836727e+130, + 2.10775729837953e+132, + 1.85482642257398e+134, + 1.65079551609085e+136, + 1.48571596448176e+138, + 1.3520015276784e+140, + 1.24384140546413e+142, + 1.15677250708164e+144, + 1.08736615665674e+146, + 1.03299784882391e+148, + 9.91677934870949e+149, + 9.61927596824821e+151, + 9.42689044888324e+153, + 9.33262154439441e+155, + 9.33262154439441e+157, + 9.42594775983835e+159, + 9.61446671503512e+161, + 9.90290071648618e+163, + 1.02990167451456e+166, + 1.08139675824029e+168, + 1.14628056373471e+170, + 1.22652020319614e+172, + 1.32464181945183e+174, + 1.44385958320249e+176, + 1.58824554152274e+178, + 1.76295255109024e+180, + 1.97450685722107e+182, + 2.23119274865981e+184, + 2.54355973347219e+186, + 2.92509369349301e+188, + 3.3931086844519e+190, + 3.96993716080872e+192, + 4.68452584975429e+194, + 5.5745857612076e+196, + 6.68950291344912e+198, + 8.09429852527344e+200, + 9.8750442008336e+202, + 1.21463043670253e+205, + 1.50614174151114e+207, + 1.88267717688893e+209, + 2.37217324288005e+211, + 3.01266001845766e+213, + 3.8562048236258e+215, + 4.97450422247729e+217, + 6.46685548922047e+219, + 8.47158069087882e+221, + 1.118248651196e+224, + 1.48727070609069e+226, + 1.99294274616152e+228, + 2.69047270731805e+230, + 3.65904288195255e+232, + 5.01288874827499e+234, + 6.91778647261949e+236, + 9.61572319694109e+238, + 1.34620124757175e+241, + 1.89814375907617e+243, + 2.69536413788816e+245, + 3.85437071718007e+247, + 5.5502938327393e+249, + 8.04792605747199e+251, + 1.17499720439091e+254, + 1.72724589045464e+256, + 2.55632391787286e+258, + 3.80892263763057e+260, + 5.71338395644585e+262, + 8.62720977423323e+264, + 1.31133588568345e+267, + 2.00634390509568e+269, + 3.08976961384735e+271, + 4.78914290146339e+273, + 7.47106292628289e+275, + 1.17295687942641e+278, + 1.85327186949373e+280, + 2.94670227249504e+282, + 4.71472363599206e+284, + 7.59070505394721e+286, + 1.22969421873945e+289, + 2.0044015765453e+291, + 3.28721858553429e+293, + 5.42391066613159e+295, + 9.00369170577843e+297, + 1.503616514865e+300, // nmaxfactorial = 167 +}; + +/* ---------------------------------------------------------------------- + factorial n vial lookup from precomputed table +------------------------------------------------------------------------- */ + +double MathSpecial::factorial(const int n) +{ + if (n < 0 || n > nmaxfactorial) + return std::numeric_limits::quiet_NaN(); + + return nfac_table[n]; +} + + /* Library libcerf: * Compute complex error functions, based on a new implementation of * Faddeeva's w_of_z. Also provide Dawson and Voigt functions. diff --git a/src/math_special.h b/src/math_special.h index 1e7b10d9fd..59517a2f76 100644 --- a/src/math_special.h +++ b/src/math_special.h @@ -20,6 +20,10 @@ namespace LAMMPS_NS { namespace MathSpecial { + // tabulated factorial function + + extern double factorial(const int); + // support function for scaled error function complement extern double erfcx_y100(const double y100); From a403e209d38db8b06a74401e8ed0109cea310771 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 13:00:51 -0500 Subject: [PATCH 16/25] use C++11 constexpr instead of C++98 const. --- src/math_const.h | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/math_const.h b/src/math_const.h index fefb442f68..53a485bf99 100644 --- a/src/math_const.h +++ b/src/math_const.h @@ -17,18 +17,18 @@ namespace LAMMPS_NS { namespace MathConst { - static const double THIRD = 1.0/3.0; - static const double TWOTHIRDS = 2.0/3.0; - static const double MY_PI = 3.14159265358979323846; // pi - static const double MY_2PI = 6.28318530717958647692; // 2pi - static const double MY_3PI = 9.42477796076937971538; // 3pi - static const double MY_4PI = 12.56637061435917295384; // 4pi - static const double MY_PI2 = 1.57079632679489661923; // pi/2 - static const double MY_PI4 = 0.78539816339744830962; // pi/4 - static const double MY_PIS = 1.77245385090551602729; // sqrt(pi) - static const double MY_ISPI4 = 1.12837916709551257389; // 1/sqrt(pi/4) - static const double MY_SQRT2 = 1.41421356237309504880; // sqrt(2) - static const double MY_CBRT2 = 1.25992104989487316476; // 2*(1/3) + static constexpr double THIRD = 1.0/3.0; + static constexpr double TWOTHIRDS = 2.0/3.0; + static constexpr double MY_PI = 3.14159265358979323846; // pi + static constexpr double MY_2PI = 6.28318530717958647692; // 2pi + static constexpr double MY_3PI = 9.42477796076937971538; // 3pi + static constexpr double MY_4PI = 12.56637061435917295384; // 4pi + static constexpr double MY_PI2 = 1.57079632679489661923; // pi/2 + static constexpr double MY_PI4 = 0.78539816339744830962; // pi/4 + static constexpr double MY_PIS = 1.77245385090551602729; // sqrt(pi) + static constexpr double MY_ISPI4 = 1.12837916709551257389; // 1/sqrt(pi/4) + static constexpr double MY_SQRT2 = 1.41421356237309504880; // sqrt(2) + static constexpr double MY_CBRT2 = 1.25992104989487316476; // 2*(1/3) } } From 2426e6245db130237a645230aaec564fa1871a91 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 14:02:56 -0500 Subject: [PATCH 17/25] move static constants in pair style out of the globally included header --- src/KOKKOS/pair_zbl_kokkos.cpp | 2 ++ src/USER-OMP/pair_zbl_omp.cpp | 1 - src/pair_zbl.cpp | 7 +++++-- src/pair_zbl.h | 17 ----------------- src/pair_zbl_const.h | 34 ++++++++++++++++++++++++++++++++++ 5 files changed, 41 insertions(+), 20 deletions(-) create mode 100644 src/pair_zbl_const.h diff --git a/src/KOKKOS/pair_zbl_kokkos.cpp b/src/KOKKOS/pair_zbl_kokkos.cpp index 2cbe65dcf7..4e7e67a074 100644 --- a/src/KOKKOS/pair_zbl_kokkos.cpp +++ b/src/KOKKOS/pair_zbl_kokkos.cpp @@ -30,6 +30,8 @@ #include "atom_masks.h" #include "kokkos.h" +#include "pair_zbl_const.h" + // From J.F. Zeigler, J. P. Biersack and U. Littmark, // "The Stopping and Range of Ions in Matter" volume 1, Pergamon, 1985. diff --git a/src/USER-OMP/pair_zbl_omp.cpp b/src/USER-OMP/pair_zbl_omp.cpp index 993e75f499..6794789f25 100644 --- a/src/USER-OMP/pair_zbl_omp.cpp +++ b/src/USER-OMP/pair_zbl_omp.cpp @@ -24,7 +24,6 @@ #include "omp_compat.h" using namespace LAMMPS_NS; -using namespace PairZBLConstants; /* ---------------------------------------------------------------------- */ diff --git a/src/pair_zbl.cpp b/src/pair_zbl.cpp index 6dff702d2c..80b623ff35 100644 --- a/src/pair_zbl.cpp +++ b/src/pair_zbl.cpp @@ -16,15 +16,18 @@ ------------------------------------------------------------------------- */ #include "pair_zbl.h" + #include + #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" +#include "memory.h" #include "neighbor.h" #include "neigh_list.h" -#include "memory.h" -#include "error.h" +#include "pair_zbl_const.h" // From J.F. Zeigler, J. P. Biersack and U. Littmark, // "The Stopping and Range of Ions in Matter" volume 1, Pergamon, 1985. diff --git a/src/pair_zbl.h b/src/pair_zbl.h index 55a5dc5fa4..6a16bc7419 100644 --- a/src/pair_zbl.h +++ b/src/pair_zbl.h @@ -54,23 +54,6 @@ class PairZBL : public Pair { double d2zbldr2(double, int, int); void set_coeff(int, int, double, double); }; - -namespace PairZBLConstants { - - // ZBL constants - - static const double pzbl = 0.23; - static const double a0 = 0.46850; - static const double c1 = 0.02817; - static const double c2 = 0.28022; - static const double c3 = 0.50986; - static const double c4 = 0.18175; - static const double d1 = 0.20162; - static const double d2 = 0.40290; - static const double d3 = 0.94229; - static const double d4 = 3.19980; -} - } #endif diff --git a/src/pair_zbl_const.h b/src/pair_zbl_const.h new file mode 100644 index 0000000000..385657693f --- /dev/null +++ b/src/pair_zbl_const.h @@ -0,0 +1,34 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_PAIR_ZBL_CONST_H +#define LMP_PAIR_ZBL_CONST_H + +namespace LAMMPS_NS { +namespace PairZBLConstants { + + // ZBL constants + + static constexpr double pzbl = 0.23; + static constexpr double a0 = 0.46850; + static constexpr double c1 = 0.02817; + static constexpr double c2 = 0.28022; + static constexpr double c3 = 0.50986; + static constexpr double c4 = 0.18175; + static constexpr double d1 = 0.20162; + static constexpr double d2 = 0.40290; + static constexpr double d3 = 0.94229; + static constexpr double d4 = 3.19980; +} +} +#endif From 3e2b004d216e9b6dfceab78afa889657a4f3f754 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 16:01:26 -0500 Subject: [PATCH 18/25] more use of constexpr --- src/MANYBODY/pair_comb.h | 2 +- src/MANYBODY/pair_comb3.h | 2 +- src/MANYBODY/pair_gw.h | 2 +- src/MANYBODY/pair_gw_zbl.h | 2 +- src/MANYBODY/pair_nb3b_harmonic.h | 2 +- src/MANYBODY/pair_sw.h | 2 +- src/MANYBODY/pair_tersoff.h | 2 +- src/MANYBODY/pair_tersoff_mod.h | 2 +- src/MANYBODY/pair_tersoff_mod_c.h | 2 +- src/MANYBODY/pair_tersoff_zbl.h | 2 +- src/MANYBODY/pair_vashishta.h | 2 +- src/RIGID/rigid_const.h | 20 ++++++++++---------- src/USER-MISC/pair_tersoff_table.h | 2 +- src/pair_coul_streitz.h | 2 +- src/text_file_reader.h | 2 +- 15 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/MANYBODY/pair_comb.h b/src/MANYBODY/pair_comb.h index 7a3d279033..f809dfaf0e 100644 --- a/src/MANYBODY/pair_comb.h +++ b/src/MANYBODY/pair_comb.h @@ -38,7 +38,7 @@ class PairComb : public Pair { virtual double yasu_char(double *, int &); double enegtot; - static const int NPARAMS_PER_LINE = 49; + static constexpr int NPARAMS_PER_LINE = 49; protected: struct Param { diff --git a/src/MANYBODY/pair_comb3.h b/src/MANYBODY/pair_comb3.h index 567859127b..b1c216fee6 100644 --- a/src/MANYBODY/pair_comb3.h +++ b/src/MANYBODY/pair_comb3.h @@ -37,7 +37,7 @@ class PairComb3 : public Pair { virtual double combqeq(double *, int &); double enegtot; - static const int NPARAMS_PER_LINE = 74; + static constexpr int NPARAMS_PER_LINE = 74; protected: // general potential parameters diff --git a/src/MANYBODY/pair_gw.h b/src/MANYBODY/pair_gw.h index 77ff5c0399..2dbf3dcf84 100644 --- a/src/MANYBODY/pair_gw.h +++ b/src/MANYBODY/pair_gw.h @@ -34,7 +34,7 @@ class PairGW : public Pair { void init_style(); double init_one(int, int); - static const int NPARAMS_PER_LINE = 17; + static constexpr int NPARAMS_PER_LINE = 17; protected: struct Param { diff --git a/src/MANYBODY/pair_gw_zbl.h b/src/MANYBODY/pair_gw_zbl.h index 95129c4eb8..de344e94f7 100644 --- a/src/MANYBODY/pair_gw_zbl.h +++ b/src/MANYBODY/pair_gw_zbl.h @@ -29,7 +29,7 @@ class PairGWZBL : public PairGW { PairGWZBL(class LAMMPS *); ~PairGWZBL() {} - static const int NPARAMS_PER_LINE = 21; + static constexpr int NPARAMS_PER_LINE = 21; private: double global_a_0; // Bohr radius for Coulomb repulsion diff --git a/src/MANYBODY/pair_nb3b_harmonic.h b/src/MANYBODY/pair_nb3b_harmonic.h index 0267b6d0e2..b68974e15f 100644 --- a/src/MANYBODY/pair_nb3b_harmonic.h +++ b/src/MANYBODY/pair_nb3b_harmonic.h @@ -34,7 +34,7 @@ class PairNb3bHarmonic : public Pair { double init_one(int, int); void init_style(); - static const int NPARAMS_PER_LINE = 6; + static constexpr int NPARAMS_PER_LINE = 6; protected: struct Param { diff --git a/src/MANYBODY/pair_sw.h b/src/MANYBODY/pair_sw.h index 441b1c0303..4fdc538430 100644 --- a/src/MANYBODY/pair_sw.h +++ b/src/MANYBODY/pair_sw.h @@ -34,7 +34,7 @@ class PairSW : public Pair { virtual double init_one(int, int); virtual void init_style(); - static const int NPARAMS_PER_LINE = 14; + static constexpr int NPARAMS_PER_LINE = 14; struct Param { double epsilon,sigma; diff --git a/src/MANYBODY/pair_tersoff.h b/src/MANYBODY/pair_tersoff.h index 8ebbc83544..eb29cc227d 100644 --- a/src/MANYBODY/pair_tersoff.h +++ b/src/MANYBODY/pair_tersoff.h @@ -34,7 +34,7 @@ class PairTersoff : public Pair { virtual void init_style(); double init_one(int, int); - static const int NPARAMS_PER_LINE = 17; + static constexpr int NPARAMS_PER_LINE = 17; protected: diff --git a/src/MANYBODY/pair_tersoff_mod.h b/src/MANYBODY/pair_tersoff_mod.h index 2cdbae9518..4f1d7a12a6 100644 --- a/src/MANYBODY/pair_tersoff_mod.h +++ b/src/MANYBODY/pair_tersoff_mod.h @@ -30,7 +30,7 @@ class PairTersoffMOD : public PairTersoff { PairTersoffMOD(class LAMMPS *); ~PairTersoffMOD() {} - static const int NPARAMS_PER_LINE = 20; + static constexpr int NPARAMS_PER_LINE = 20; protected: virtual void read_file(char *); diff --git a/src/MANYBODY/pair_tersoff_mod_c.h b/src/MANYBODY/pair_tersoff_mod_c.h index bad5f988be..4949cb76db 100644 --- a/src/MANYBODY/pair_tersoff_mod_c.h +++ b/src/MANYBODY/pair_tersoff_mod_c.h @@ -29,7 +29,7 @@ class PairTersoffMODC : public PairTersoffMOD { PairTersoffMODC(class LAMMPS *lmp) : PairTersoffMOD(lmp) {}; ~PairTersoffMODC() {} - static const int NPARAMS_PER_LINE = 21; + static constexpr int NPARAMS_PER_LINE = 21; protected: void read_file(char *); diff --git a/src/MANYBODY/pair_tersoff_zbl.h b/src/MANYBODY/pair_tersoff_zbl.h index be3f496b62..42483b7d89 100644 --- a/src/MANYBODY/pair_tersoff_zbl.h +++ b/src/MANYBODY/pair_tersoff_zbl.h @@ -29,7 +29,7 @@ class PairTersoffZBL : public PairTersoff { PairTersoffZBL(class LAMMPS *); ~PairTersoffZBL() {} - static const int NPARAMS_PER_LINE = 21; + static constexpr int NPARAMS_PER_LINE = 21; protected: double global_a_0; // Bohr radius for Coulomb repulsion diff --git a/src/MANYBODY/pair_vashishta.h b/src/MANYBODY/pair_vashishta.h index bd250d328a..19a8e01cd9 100644 --- a/src/MANYBODY/pair_vashishta.h +++ b/src/MANYBODY/pair_vashishta.h @@ -34,7 +34,7 @@ class PairVashishta : public Pair { double init_one(int, int); void init_style(); - static const int NPARAMS_PER_LINE = 17; + static constexpr int NPARAMS_PER_LINE = 17; struct Param { double bigb,gamma,r0,bigc,costheta; diff --git a/src/RIGID/rigid_const.h b/src/RIGID/rigid_const.h index 3aae988197..b345a6b9dc 100644 --- a/src/RIGID/rigid_const.h +++ b/src/RIGID/rigid_const.h @@ -33,21 +33,21 @@ namespace LAMMPS_NS { TORQUE = 1<<8 }; - static const double TOLERANCE = 1.0e-6; - static const double EPSILON = 1.0e-7; - static const double BIG = 1.0e20; + static constexpr double TOLERANCE = 1.0e-6; + static constexpr double EPSILON = 1.0e-7; + static constexpr double BIG = 1.0e20; // moment of inertia prefactor for sphere - static const double SINERTIA = 0.4; + static constexpr double SINERTIA = 0.4; // moment of inertia prefactor for ellipsoid - static const double EINERTIA = 0.2; + static constexpr double EINERTIA = 0.2; // moment of inertia prefactor for line segment - static const double LINERTIA = 1.0/12.0; + static constexpr double LINERTIA = 1.0/12.0; - static const int MAXLINE = 1024; - static const int CHUNK = 1024; - static const int DELTA_BODY = 10000; - static const int ATTRIBUTE_PERBODY = 20; + static constexpr int MAXLINE = 1024; + static constexpr int CHUNK = 1024; + static constexpr int DELTA_BODY = 10000; + static constexpr int ATTRIBUTE_PERBODY = 20; } } diff --git a/src/USER-MISC/pair_tersoff_table.h b/src/USER-MISC/pair_tersoff_table.h index faa704191c..d63a5be5a6 100644 --- a/src/USER-MISC/pair_tersoff_table.h +++ b/src/USER-MISC/pair_tersoff_table.h @@ -43,7 +43,7 @@ class PairTersoffTable : public Pair { void init_style(); double init_one(int, int); - static const int NPARAMS_PER_LINE = 17; + static constexpr int NPARAMS_PER_LINE = 17; protected: struct Param { diff --git a/src/pair_coul_streitz.h b/src/pair_coul_streitz.h index d5a4a2de5c..2f62846212 100644 --- a/src/pair_coul_streitz.h +++ b/src/pair_coul_streitz.h @@ -36,7 +36,7 @@ class PairCoulStreitz : public Pair { double memory_usage(); virtual void *extract(const char *, int &); - static const int NPARAMS_PER_LINE = 6; + static constexpr int NPARAMS_PER_LINE = 6; protected: struct Param { diff --git a/src/text_file_reader.h b/src/text_file_reader.h index 0b90304911..0da21e4581 100644 --- a/src/text_file_reader.h +++ b/src/text_file_reader.h @@ -27,7 +27,7 @@ namespace LAMMPS_NS class TextFileReader { std::string filename; std::string filetype; - static const int MAXLINE = 1024; + static constexpr int MAXLINE = 1024; char line[MAXLINE]; FILE *fp; From f6fa564ef2b587a48aee0d0f2395b527be009236 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 17:43:04 -0500 Subject: [PATCH 19/25] whitespace --- src/math_special.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/math_special.cpp b/src/math_special.cpp index 4ed91ba51f..243d5a05f3 100644 --- a/src/math_special.cpp +++ b/src/math_special.cpp @@ -183,7 +183,7 @@ static const double nfac_table[] = { 5.42391066613159e+295, 9.00369170577843e+297, 1.503616514865e+300, // nmaxfactorial = 167 -}; +}; /* ---------------------------------------------------------------------- factorial n vial lookup from precomputed table From de94e28c8bcf07510f1f6426914a48c4231b6dd6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 17 Dec 2020 21:48:04 -0500 Subject: [PATCH 20/25] correct path to find liblammps.dll with Windows installer package --- python/lammps/core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/lammps/core.py b/python/lammps/core.py index 161583b78c..c0cbaac533 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -77,7 +77,7 @@ class lammps(object): modpath = dirname(abspath(getsourcefile(lambda:0))) # for windows installers the shared library is in a different folder - winpath = abspath(os.path.join(modpath,'..','bin')) + winpath = abspath(os.path.join(modpath,'..','..','bin')) self.lib = None self.lmp = None From 71e340a792607f0ede5874088c895583ed9acf2a Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Fri, 18 Dec 2020 17:56:51 -0700 Subject: [PATCH 21/25] Adding multi+granular note in neighbor documentation --- doc/src/neighbor.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index 1bb591587c..6671e7e1ac 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -49,7 +49,9 @@ sometimes be faster. Either style should give the same answers. The *multi* style is a modified binning algorithm that is useful for systems with a wide range of cutoff distances, e.g. due to different -size particles. For the *bin* style, the bin size is set to 1/2 of +size particles. For granular pairstyles, cutoffs are set to the +sum of the maximum atomic radii for each atom type. +For the *bin* style, the bin size is set to 1/2 of the largest cutoff distance between any pair of atom types and a single set of bins is defined to search over for all atom types. This can be inefficient if one pair of types has a very long cutoff, but From 384376bc51cc37f64edc2d612bd47205e7eca9c7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Dec 2020 22:33:56 -0500 Subject: [PATCH 22/25] increase precision of lj/cubic constants and document how to compute them --- src/GPU/pair_lj_cubic_gpu.cpp | 5 +- src/USER-OMP/pair_lj_cubic_omp.cpp | 2 + src/pair_lj_cubic.cpp | 1 + src/pair_lj_cubic.h | 14 -- src/pair_lj_cubic_const.h | 48 +++++++ .../force-styles/tests/mol-pair-lj_cubic.yaml | 128 +++++++++--------- 6 files changed, 119 insertions(+), 79 deletions(-) create mode 100644 src/pair_lj_cubic_const.h diff --git a/src/GPU/pair_lj_cubic_gpu.cpp b/src/GPU/pair_lj_cubic_gpu.cpp index a669d52a19..35062a5d71 100644 --- a/src/GPU/pair_lj_cubic_gpu.cpp +++ b/src/GPU/pair_lj_cubic_gpu.cpp @@ -16,10 +16,11 @@ ------------------------------------------------------------------------- */ #include "pair_lj_cubic_gpu.h" + #include #include - #include + #include "atom.h" #include "atom_vec.h" #include "comm.h" @@ -36,6 +37,8 @@ #include "gpu_extra.h" #include "suffix.h" +#include "pair_lj_cubic_const.h" + using namespace LAMMPS_NS; using namespace PairLJCubicConstants; diff --git a/src/USER-OMP/pair_lj_cubic_omp.cpp b/src/USER-OMP/pair_lj_cubic_omp.cpp index 872574a3b2..e96b3d0a79 100644 --- a/src/USER-OMP/pair_lj_cubic_omp.cpp +++ b/src/USER-OMP/pair_lj_cubic_omp.cpp @@ -23,6 +23,8 @@ #include #include "omp_compat.h" +#include "pair_lj_cubic_const.h" + using namespace LAMMPS_NS; using namespace PairLJCubicConstants; diff --git a/src/pair_lj_cubic.cpp b/src/pair_lj_cubic.cpp index 7fe4757358..c5f1f5b257 100644 --- a/src/pair_lj_cubic.cpp +++ b/src/pair_lj_cubic.cpp @@ -26,6 +26,7 @@ #include "memory.h" #include "error.h" +#include "pair_lj_cubic_const.h" using namespace LAMMPS_NS; using namespace PairLJCubicConstants; diff --git a/src/pair_lj_cubic.h b/src/pair_lj_cubic.h index d0af06c0a0..4b5edf7e97 100644 --- a/src/pair_lj_cubic.h +++ b/src/pair_lj_cubic.h @@ -46,20 +46,6 @@ class PairLJCubic : public Pair { void allocate(); }; - -namespace PairLJCubicConstants { - - // LJ quantities scaled by epsilon and rmin = sigma*2^1/6 - - static const double RT6TWO = 1.1224621; // 2^1/6 - static const double SS = 1.1086834; // inflection point (13/7)^1/6 - static const double PHIS = -0.7869823; // energy at s - static const double DPHIDS = 2.6899009; // gradient at s - static const double A3 = 27.93357; // cubic coefficient - static const double SM = 1.5475375; // cubic cutoff = s*67/48 - -} - } #endif diff --git a/src/pair_lj_cubic_const.h b/src/pair_lj_cubic_const.h new file mode 100644 index 0000000000..5fc773ac3c --- /dev/null +++ b/src/pair_lj_cubic_const.h @@ -0,0 +1,48 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_PAIR_LJ_CUBIC_CONST_H +#define LMP_PAIR_LJ_CUBIC_CONST_H + +namespace LAMMPS_NS { +namespace PairLJCubicConstants { + + // LJ quantities scaled by epsilon and rmin = sigma*2^1/6 + + static constexpr double RT6TWO = 1.1224620483093730; // 2^1/6 + static constexpr double SS = 1.1086834179687215; // inflection point (13/7)^1/6 + static constexpr double PHIS = -0.7869822485207097; // energy at s + static constexpr double DPHIDS = 2.6899008972047196; // gradient at s + static constexpr double A3 = 27.9335700460986445; // cubic coefficient + static constexpr double SM = 1.5475372709146737; // cubic cutoff = s*67/48} +} +} +#endif + +// python script to compute the constants +// +// sixth = 1.0/6.0 +// rmin = pow(2.0,sixth) +// rs = pow(26.0/7.0,sixth) +// ss = rs/rmin +// pow6 = pow(1.0/rs,6.0) +// phis = 4.0*pow6*(pow6-1.0) +// dpds = -24.0*pow6*(2.0*pow6-1.0)/rs*rmin +// a3 = 8.0*pow(dpds,3)/(9.0*phis*phis) +// sm = 67.0/48.0*ss +// print("static constexpr double RT6TWO = %19.16f; // 2^1/6" % rmin) +// print("static constexpr double SS = %19.16f; // inflection point (13/7)^1/6" % ss) +// print("static constexpr double PHIS = %19.16f; // energy at s" % phis) +// print("static constexpr double DPHIDS = %19.16f; // gradient at s" % dpds) +// print("static constexpr double A3 = %19.16f; // cubic coefficient" % a3) +// print("static constexpr double SM = %19.16f; // cubic cutoff = s*67/48" % sm) diff --git a/unittest/force-styles/tests/mol-pair-lj_cubic.yaml b/unittest/force-styles/tests/mol-pair-lj_cubic.yaml index ad7f0666ca..6d23998ad5 100644 --- a/unittest/force-styles/tests/mol-pair-lj_cubic.yaml +++ b/unittest/force-styles/tests/mol-pair-lj_cubic.yaml @@ -1,7 +1,7 @@ --- -lammps_version: 24 Aug 2020 -date_generated: Tue Sep 15 09:44:15 202 -epsilon: 2e-13 +lammps_version: 30 Nov 2020 +date_generated: Fri Dec 18 22:30:42 202 +epsilon: 1e-13 prerequisites: ! | atom full pair lj/cubic @@ -19,72 +19,72 @@ pair_coeff: ! | 5 5 0.015 3.1 extract: ! "" natoms: 29 -init_vdwl: 166.005266308562 +init_vdwl: 166.005266338502 init_coul: 0 init_stress: ! |2- - 4.5860676757908186e+02 4.8091912919212928e+02 1.0767204080701006e+03 -2.1005546139122362e+02 -2.9491286717936713e+00 1.6145675857120941e+02 + 4.5860676760613336e+02 4.8091912923656065e+02 1.0767204080957010e+03 -2.1005546139899724e+02 -2.9491286649299440e+00 1.6145675855773089e+02 init_forces: ! |2 - 1 9.1849370411551270e+00 7.6268937957720553e+01 6.1726872441625311e+01 - 2 2.2858712118514426e+01 1.8809274242266209e+01 -2.6905829837199740e+01 - 3 -3.2016987482543328e+01 -9.4135849525427091e+01 -3.4799279593035926e+01 - 4 -5.5341015869901478e-01 1.5206999898436971e-01 -3.9418368928369890e-01 - 5 -1.8042057425348118e-01 -3.0459951056385326e-01 8.7068483241007189e-01 - 6 -2.0038994438822397e+02 2.3344446299945159e+02 2.8487343926572851e+02 - 7 8.0909912172413883e+00 -7.8410849891085633e+01 -4.3214084684451740e+02 - 8 4.7943581255133857e+01 -2.1287511456246008e+01 1.4094503445180061e+02 - 9 1.1447552368270737e+01 1.2328709806786962e+01 5.0656476982000299e+01 - 10 1.3071496571967870e+02 -1.4589264560693914e+02 -4.4748155922123622e+01 - 11 -1.6551880116149281e-01 -4.1534332040572380e-01 -6.8284765241715795e-01 - 12 1.7721533626133388e+00 6.3456329073685158e-01 -8.2372301448028962e-01 - 13 5.6789360334118277e-01 -2.2634410312439054e-01 -9.7536738055328392e-03 - 14 -2.4337021468262635e-01 4.6659433642728905e-02 -6.1110664501270184e-01 - 15 -2.1936997101927893e-02 5.9238263972968364e-01 2.1493099548264527e-01 - 16 1.1121534968449923e+02 -7.8056927924992834e+01 -2.9249212971206231e+02 - 17 -1.1020604609843586e+02 7.6481296254913858e+01 2.9430701446263464e+02 - 18 -1.6570656719723909e-02 -2.7996966177077785e-02 2.6456326954440619e-02 - 19 7.4243353115058947e-04 6.3524893127716046e-04 1.8675586277048476e-04 - 20 -7.4243353115058947e-04 -6.3524893127716046e-04 -1.8675586277048476e-04 + 1 9.1849370400954342e+00 7.6268937960282486e+01 6.1726872440953251e+01 + 2 2.2858712118557477e+01 1.8809274242339992e+01 -2.6905829837209875e+01 + 3 -3.2016987482586380e+01 -9.4135849525500873e+01 -3.4799279593025787e+01 + 4 -5.5341015917141478e-01 1.5206999930354809e-01 -3.9418368927471259e-01 + 5 -1.8042057490785210e-01 -3.0459951013385317e-01 8.7068483295449139e-01 + 6 -2.0038994438983289e+02 2.3344446299845976e+02 2.8487343926562755e+02 + 7 8.0909912152126484e+00 -7.8410849888970290e+01 -4.3214084684693944e+02 + 8 4.7943581253518161e+01 -2.1287511458762772e+01 1.4094503445043824e+02 + 9 1.1447552368334570e+01 1.2328709806695397e+01 5.0656476981938759e+01 + 10 1.3071496571895125e+02 -1.4589264560750766e+02 -4.4748155921681018e+01 + 11 -1.6551880130741767e-01 -4.1534332035668403e-01 -6.8284765271893910e-01 + 12 1.7721533633976141e+00 6.3456328808910523e-01 -8.2372301296999062e-01 + 13 5.6789360517170440e-01 -2.2634410322310763e-01 -9.7536744156455080e-03 + 14 -2.4337021462140446e-01 4.6659433737428486e-02 -6.1110664517858926e-01 + 15 -2.1936996885396808e-02 5.9238264018028652e-01 2.1493099542010694e-01 + 16 1.1121534968641315e+02 -7.8056927926703892e+01 -2.9249212971001145e+02 + 17 -1.1020604609471918e+02 7.6481296251709367e+01 2.9430701446464468e+02 + 18 -1.6570656172995385e-02 -2.7996961634495443e-02 2.6456324093123054e-02 + 19 7.4243314752471426e-04 6.3524860303507949e-04 1.8675576627109685e-04 + 20 -7.4243314752471426e-04 -6.3524860303507949e-04 -1.8675576627109685e-04 21 -1.1415041486189516e+01 -1.3016363071591645e+01 3.6007276733401099e+01 - 22 -1.7227422089792942e+01 -4.1746638094950628e+00 -2.7029162034499002e+01 - 23 2.8642463575982458e+01 1.7191026881086707e+01 -8.9781146989020968e+00 - 24 5.8150644491939154e+00 -3.3774314134628064e+01 1.7867788752379695e+01 - 25 -2.3666545027773044e+01 3.8106021846559952e+00 -1.9896269873584632e+01 - 26 1.7843812244577855e+01 2.9960339884741117e+01 2.0167430316952100e+00 - 27 8.2825859209946024e+00 -3.6194570066818969e+01 1.4492694351988913e+01 - 28 -2.8773892796642542e+01 1.2366374307374247e+01 -1.9468877181285176e+01 - 29 2.0497044211022661e+01 2.3831279505404666e+01 4.9748677441078746e+00 -run_vdwl: 164.176727313193 + 22 -1.7227422090167991e+01 -4.1746638096674396e+00 -2.7029162034658786e+01 + 23 2.8642463576357507e+01 1.7191026881259084e+01 -8.9781146987423170e+00 + 24 5.8150644495067958e+00 -3.3774314132751599e+01 1.7867788754173183e+01 + 25 -2.3666545028156815e+01 3.8106021844353424e+00 -1.9896269873795571e+01 + 26 1.7843812244961626e+01 2.9960339884961773e+01 2.0167430319061479e+00 + 27 8.2825859198612442e+00 -3.6194570067428131e+01 1.4492694352248694e+01 + 28 -2.8773892797077618e+01 1.2366374307246014e+01 -1.9468877181493664e+01 + 29 2.0497044211457737e+01 2.3831279505532898e+01 4.9748677443163629e+00 +run_vdwl: 164.176727343123 run_coul: 0 run_stress: ! |2- - 4.5696669905868288e+02 4.7904134871830234e+02 1.0582153129928076e+03 -2.0794233475912671e+02 -2.0206236780739086e+00 1.5948889983617320e+02 + 4.5696669908578690e+02 4.7904134876274537e+02 1.0582153130184231e+03 -2.0794233476691440e+02 -2.0206236712071677e+00 1.5948889982266843e+02 run_forces: ! |2 - 1 9.2483453892728740e+00 7.5945844239974676e+01 6.1367289784738311e+01 - 2 2.2732852134554197e+01 1.8737493288759875e+01 -2.6669600068587577e+01 - 3 -3.1964986246503269e+01 -9.3734416340601200e+01 -3.4671255449722295e+01 - 4 -5.5001585694065913e-01 1.5070776418957152e-01 -3.9275226164990551e-01 - 5 -1.7969915375876547e-01 -3.0362292678324765e-01 8.6793365916706400e-01 - 6 -1.9796806590003533e+02 2.2990162655369829e+02 2.7501790745035373e+02 - 7 7.9109888171861886e+00 -7.6524641847864288e+01 -4.2032627555539375e+02 - 8 4.5992175211879918e+01 -1.9873148355500923e+01 1.3922798515578887e+02 - 9 1.1403303962503362e+01 1.2231165268449960e+01 5.0409090529604853e+01 - 10 1.3047784096912977e+02 -1.4556513307073578e+02 -4.4756481385998420e+01 - 11 -1.6346238349194633e-01 -4.0922088697872150e-01 -6.7303941060221573e-01 - 12 1.7689668574627495e+00 6.3166692380469824e-01 -8.3233385524693426e-01 - 13 5.6480104331071312e-01 -2.2395151872256039e-01 -9.5790993973179691e-03 - 14 -2.4030540495346994e-01 4.5179188229734012e-02 -6.0304369153488313e-01 - 15 -2.3220111317111478e-02 5.9408611078423390e-01 2.1676726911830960e-01 - 16 1.0963039832302356e+02 -7.7096357855469549e+01 -2.8842624961188653e+02 - 17 -1.0862142117090467e+02 7.5521002117836630e+01 2.9024023117219969e+02 - 18 -1.6565212275867415e-02 -2.7990268691876326e-02 2.6458602067006932e-02 - 19 7.1473709241289744e-04 6.1248437925700357e-04 1.7889258733572757e-04 - 20 -7.1473709241289744e-04 -6.1248437925700357e-04 -1.7889258733572757e-04 - 21 -1.1536904971247665e+01 -1.3021625993962397e+01 3.6108894191673429e+01 - 22 -1.7333879764643559e+01 -4.2314763344327275e+00 -2.7103019136756011e+01 - 23 2.8870784735891224e+01 1.7253102328395126e+01 -9.0058750549174214e+00 - 24 6.1437425316795213e+00 -3.4297207023632204e+01 1.8296742414004438e+01 - 25 -2.4276461284621075e+01 3.8560435260643189e+00 -2.0415720860228767e+01 - 26 1.8125049871956215e+01 3.0437790982988908e+01 2.1072387594169975e+00 - 27 8.4078124309265192e+00 -3.6323119973255714e+01 1.4505938075037919e+01 - 28 -2.8937319168272772e+01 1.2421253477627801e+01 -1.9540501416079319e+01 - 29 2.0535244350189391e+01 2.3904950625827414e+01 5.0332497948309163e+00 + 1 9.2483453882111135e+00 7.5945844242532630e+01 6.1367289784064468e+01 + 2 2.2732852134599622e+01 1.8737493288836678e+01 -2.6669600068598680e+01 + 3 -3.1964986246546172e+01 -9.3734416340673434e+01 -3.4671255449709932e+01 + 4 -5.5001585741105219e-01 1.5070776450703974e-01 -3.9275226164162935e-01 + 5 -1.7969915441332354e-01 -3.0362292635308630e-01 8.6793365971015668e-01 + 6 -1.9796806590165264e+02 2.2990162655270728e+02 2.7501790745023669e+02 + 7 7.9109888151537930e+00 -7.6524641845739424e+01 -4.2032627555780203e+02 + 8 4.5992175210254615e+01 -1.9873148358014198e+01 1.3922798515443196e+02 + 9 1.1403303962567387e+01 1.2231165268357689e+01 5.0409090529541324e+01 + 10 1.3047784096840681e+02 -1.4556513307130629e+02 -4.4756481385556341e+01 + 11 -1.6346238363922402e-01 -4.0922088693014880e-01 -6.7303941090576869e-01 + 12 1.7689668582530251e+00 6.3166692115802814e-01 -8.3233385373913271e-01 + 13 5.6480104514758445e-01 -2.2395151881989347e-01 -9.5791000096411370e-03 + 14 -2.4030540489083232e-01 4.5179188324553629e-02 -6.0304369170302041e-01 + 15 -2.3220111102749796e-02 5.9408611123108124e-01 2.1676726905655996e-01 + 16 1.0963039832494319e+02 -7.7096357857185083e+01 -2.8842624960984045e+02 + 17 -1.0862142116718823e+02 7.5521002114629724e+01 2.9024023117422536e+02 + 18 -1.6565211729771653e-02 -2.7990264151000276e-02 2.6458599205916818e-02 + 19 7.1473670261989000e-04 6.1248404522910675e-04 1.7889248977386897e-04 + 20 -7.1473670261989000e-04 -6.1248404522910675e-04 -1.7889248977386897e-04 + 21 -1.1536904971247079e+01 -1.3021625993961788e+01 3.6108894191671638e+01 + 22 -1.7333879765020438e+01 -4.2314763346057394e+00 -2.7103019136915339e+01 + 23 2.8870784736267517e+01 1.7253102328567529e+01 -9.0058750547563005e+00 + 24 6.1437425319896413e+00 -3.4297207021755547e+01 1.8296742415795482e+01 + 25 -2.4276461285000558e+01 3.8560435258438814e+00 -2.0415720860437396e+01 + 26 1.8125049872338021e+01 3.0437790983209180e+01 2.1072387596271578e+00 + 27 8.4078124297937791e+00 -3.6323119973862774e+01 1.4505938075297303e+01 + 28 -2.8937319168706676e+01 1.2421253477499173e+01 -1.9540501416287395e+01 + 29 2.0535244350622666e+01 2.3904950625953887e+01 5.0332497950390769e+00 ... From 2022dc0aa9b0fee7c0150b6887bf6940c43fecf6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Dec 2020 22:43:12 -0500 Subject: [PATCH 23/25] whitespace fixes --- doc/src/comm_modify.rst | 6 +- doc/src/neighbor.rst | 2 +- src/balance.cpp | 92 ++++++++++++------------ src/npair_half_size_multi_newton.cpp | 6 +- src/npair_half_size_multi_newton_tri.cpp | 2 +- 5 files changed, 54 insertions(+), 54 deletions(-) diff --git a/doc/src/comm_modify.rst b/doc/src/comm_modify.rst index 1ccf77e995..56d3c08a57 100644 --- a/doc/src/comm_modify.rst +++ b/doc/src/comm_modify.rst @@ -84,15 +84,15 @@ information is available, then also a heuristic based on that bond length is computed. It is used as communication cutoff, if there is no pair style present and no *comm_modify cutoff* command used. Otherwise a warning is printed, if this bond based estimate is larger than the -communication cutoff used. +communication cutoff used. The *cutoff/multi* option is equivalent to *cutoff*\ , but applies to communication mode *multi* instead. Since in this case the communication cutoffs are determined per atom type, a type specifier is needed and cutoff for one or multiple types can be extended. Also ranges of types -using the usual asterisk notation can be given. For granular pairstyles, +using the usual asterisk notation can be given. For granular pairstyles, the default cutoff is set to the sum of the current maximum atomic radii -for each type. +for each type. These are simulation scenarios in which it may be useful or even necessary to set a ghost cutoff > neighbor cutoff: diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index 6671e7e1ac..878fa53f8c 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -49,7 +49,7 @@ sometimes be faster. Either style should give the same answers. The *multi* style is a modified binning algorithm that is useful for systems with a wide range of cutoff distances, e.g. due to different -size particles. For granular pairstyles, cutoffs are set to the +size particles. For granular pairstyles, cutoffs are set to the sum of the maximum atomic radii for each atom type. For the *bin* style, the bin size is set to 1/2 of the largest cutoff distance between any pair of atom types and a diff --git a/src/balance.cpp b/src/balance.cpp index c47412cc9e..6294e023b3 100644 --- a/src/balance.cpp +++ b/src/balance.cpp @@ -924,66 +924,66 @@ int Balance::shift() i = 0; while (i < np) { if (split[i+1] - split[i] < close) { - j = i+1; + j = i+1; - // I,J = set of consecutive splits that are collectively too close - // if can expand set and not become too close to splits I-1 or J+1, do it - // else add split I-1 or J+1 to set and try again - // delta = size of expanded split set that will satisy criterion + // I,J = set of consecutive splits that are collectively too close + // if can expand set and not become too close to splits I-1 or J+1, do it + // else add split I-1 or J+1 to set and try again + // delta = size of expanded split set that will satisy criterion - while (1) { - delta = (j-i) * close; - midpt = 0.5 * (split[i]+split[j]); - start = midpt - 0.5*delta; - stop = midpt + 0.5*delta; + while (1) { + delta = (j-i) * close; + midpt = 0.5 * (split[i]+split[j]); + start = midpt - 0.5*delta; + stop = midpt + 0.5*delta; - if (i > 0) lbound = split[i-1] + close; - else lbound = 0.0; - if (j < np) ubound = split[j+1] - close; - else ubound = 1.0; + if (i > 0) lbound = split[i-1] + close; + else lbound = 0.0; + if (j < np) ubound = split[j+1] - close; + else ubound = 1.0; - // start/stop are within bounds, reset the splits + // start/stop are within bounds, reset the splits - if (start >= lbound && stop <= ubound) break; + if (start >= lbound && stop <= ubound) break; - // try a shift to either bound, reset the splits if delta fits - // these tests change start/stop + // try a shift to either bound, reset the splits if delta fits + // these tests change start/stop - if (start < lbound) { - start = lbound; - stop = start + delta; - if (stop <= ubound) break; - } else if (stop > ubound) { - stop = ubound; - start = stop - delta; - if (start >= lbound) break; - } + if (start < lbound) { + start = lbound; + stop = start + delta; + if (stop <= ubound) break; + } else if (stop > ubound) { + stop = ubound; + start = stop - delta; + if (start >= lbound) break; + } - // delta does not fit between lbound and ubound - // exit if can't expand set, else expand set - // if can expand in either direction, - // pick new split closest to current midpt of set + // delta does not fit between lbound and ubound + // exit if can't expand set, else expand set + // if can expand in either direction, + // pick new split closest to current midpt of set - if (i == 0 && j == np) { - start = 0.0; stop = 1.0; - break; - } - if (i == 0) j++; - else if (j == np) i--; - else if (midpt-lbound < ubound-midpt) i--; - else j++; - } + if (i == 0 && j == np) { + start = 0.0; stop = 1.0; + break; + } + if (i == 0) j++; + else if (j == np) i--; + else if (midpt-lbound < ubound-midpt) i--; + else j++; + } - // reset all splits between I,J inclusive to be equi-spaced + // reset all splits between I,J inclusive to be equi-spaced - spacing = (stop-start) / (j-i); - for (m = i; m <= j; m++) - split[m] = start + (m-i)*spacing; + spacing = (stop-start) / (j-i); + for (m = i; m <= j; m++) + split[m] = start + (m-i)*spacing; if (j == np) split[np] = 1.0; - // continue testing beyond the J split + // continue testing beyond the J split - i = j+1; + i = j+1; } else i++; } diff --git a/src/npair_half_size_multi_newton.cpp b/src/npair_half_size_multi_newton.cpp index 0d12924629..943965ab63 100644 --- a/src/npair_half_size_multi_newton.cpp +++ b/src/npair_half_size_multi_newton.cpp @@ -93,9 +93,9 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) + if (history && rsq < radsum*radsum) neighptr[n++] = j ^ mask_history; - else + else neighptr[n++] = j; } } @@ -123,7 +123,7 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) + if (history && rsq < radsum*radsum) neighptr[n++] = j ^ mask_history; else neighptr[n++] = j; diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp index adc4c080ec..e59a3f088a 100644 --- a/src/npair_half_size_multi_newton_tri.cpp +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -105,7 +105,7 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) + if (history && rsq < radsum*radsum) neighptr[n++] = j ^ mask_history; else neighptr[n++] = j; From 151110f07f59ebba7b4b75470797e09a9680b22c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Dec 2020 22:55:15 -0500 Subject: [PATCH 24/25] refactor constant vs define in edip pair styles --- src/USER-MISC/pair_edip.cpp | 2 +- src/USER-MISC/pair_edip_multi.cpp | 5 ++--- src/USER-MISC/pair_edip_multi.h | 26 +++++++++++--------------- src/USER-OMP/pair_edip_omp.cpp | 2 +- 4 files changed, 15 insertions(+), 20 deletions(-) diff --git a/src/USER-MISC/pair_edip.cpp b/src/USER-MISC/pair_edip.cpp index b29f499de8..12390ba69d 100644 --- a/src/USER-MISC/pair_edip.cpp +++ b/src/USER-MISC/pair_edip.cpp @@ -46,7 +46,7 @@ using namespace LAMMPS_NS; // max number of interaction per atom for f(Z) environment potential -#define leadDimInteractionList 64 +static constexpr int leadDimInteractionList = 64; /* ---------------------------------------------------------------------- */ diff --git a/src/USER-MISC/pair_edip_multi.cpp b/src/USER-MISC/pair_edip_multi.cpp index 811f80abfc..9d36db7041 100644 --- a/src/USER-MISC/pair_edip_multi.cpp +++ b/src/USER-MISC/pair_edip_multi.cpp @@ -32,13 +32,11 @@ #include "error.h" #include "citeme.h" - using namespace LAMMPS_NS; #define MAXLINE 1024 #define DELTA 4 - static const char cite_pair_edip[] = "@article{cjiang2012\n" " author = {Jian, Chao and Morgan, Dane, and Szlufarska, Izabella},\n" @@ -56,7 +54,9 @@ static const char cite_pair_edip[] = " year = {2010},\n" "}\n\n"; +// max number of interaction per atom for f(Z) environment potential +static constexpr int leadDimInteractionList = 64; /* ---------------------------------------------------------------------- */ @@ -94,7 +94,6 @@ PairEDIPMulti::~PairEDIPMulti() memory->destroy(cutsq); delete [] map; -//XXX deallocateGrids(); deallocatePreLoops(); } } diff --git a/src/USER-MISC/pair_edip_multi.h b/src/USER-MISC/pair_edip_multi.h index 26433a66fa..5120ece5d2 100644 --- a/src/USER-MISC/pair_edip_multi.h +++ b/src/USER-MISC/pair_edip_multi.h @@ -36,17 +36,17 @@ class PairEDIPMulti : public Pair { protected: struct Param { - double A, B;//coefficients for pair interaction I-J - double cutoffA;//cut-off distance for pair interaction I-J - double cutoffC;//lower cut-off distance for calculating Z_I - double alpha;//coefficient for calculating Z_I - double beta;//attractive term for pair I-J - double sigma;//cut-off coefficient for pair I-J - double rho;//pair I-J - double gamma;//coefficient for three-body interaction I-J-K - double eta, lambda;//coefficients for function h(l,Z) - double mu, Q0;//coefficients for function Q(Z) - double u1, u2, u3, u4;//coefficients for function tau(Z) + double A, B; // coefficients for pair interaction I-J + double cutoffA; // cut-off distance for pair interaction I-J + double cutoffC; // lower cut-off distance for calculating Z_I + double alpha; // coefficient for calculating Z_I + double beta; // attractive term for pair I-J + double sigma; // cut-off coefficient for pair I-J + double rho; // pair I-J + double gamma; // coefficient for three-body interaction I-J-K + double eta, lambda; // coefficients for function h(l,Z) + double mu, Q0; // coefficients for function Q(Z) + double u1, u2, u3, u4; // coefficients for function tau(Z) double cutsq; int ielement,jelement,kelement; }; @@ -62,10 +62,6 @@ class PairEDIPMulti : public Pair { int maxparam; // max # of parameter sets Param *params; // parameter set for an I-J-K interaction - // max number of interaction per atom for f(Z) environment potential - - static const int leadDimInteractionList = 64; - void allocate(); void allocatePreLoops(void); void deallocatePreLoops(void); diff --git a/src/USER-OMP/pair_edip_omp.cpp b/src/USER-OMP/pair_edip_omp.cpp index 2eb8384113..106d038743 100644 --- a/src/USER-OMP/pair_edip_omp.cpp +++ b/src/USER-OMP/pair_edip_omp.cpp @@ -28,7 +28,7 @@ using namespace LAMMPS_NS; // max number of interaction per atom for f(Z) environment potential -#define leadDimInteractionList 64 +static constexpr int leadDimInteractionList = 64; /* ---------------------------------------------------------------------- */ From a13f70790a279a629bb2082aa8874b6854ee8152 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Dec 2020 23:02:28 -0500 Subject: [PATCH 25/25] spelling --- doc/src/comm_modify.rst | 2 +- doc/src/neighbor.rst | 8 +++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/doc/src/comm_modify.rst b/doc/src/comm_modify.rst index 56d3c08a57..864c544fed 100644 --- a/doc/src/comm_modify.rst +++ b/doc/src/comm_modify.rst @@ -90,7 +90,7 @@ The *cutoff/multi* option is equivalent to *cutoff*\ , but applies to communication mode *multi* instead. Since in this case the communication cutoffs are determined per atom type, a type specifier is needed and cutoff for one or multiple types can be extended. Also ranges of types -using the usual asterisk notation can be given. For granular pairstyles, +using the usual asterisk notation can be given. For granular pair styles, the default cutoff is set to the sum of the current maximum atomic radii for each type. diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index 878fa53f8c..98ee0d0b6a 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -49,7 +49,7 @@ sometimes be faster. Either style should give the same answers. The *multi* style is a modified binning algorithm that is useful for systems with a wide range of cutoff distances, e.g. due to different -size particles. For granular pairstyles, cutoffs are set to the +size particles. For granular pair styles, cutoffs are set to the sum of the maximum atomic radii for each atom type. For the *bin* style, the bin size is set to 1/2 of the largest cutoff distance between any pair of atom types and a @@ -59,8 +59,10 @@ other type pairs have a much shorter cutoff. For style *multi* the bin size is set to 1/2 of the shortest cutoff distance and multiple sets of bins are defined to search over for different atom types. This imposes some extra setup overhead, but the searches themselves -may be much faster for the short-cutoff cases. See the :doc:`comm_modify mode multi ` command for a communication option -that may also be beneficial for simulations of this kind. +may be much faster for the short-cutoff cases. +See the :doc:`comm_modify mode multi ` command for a +communication option that may also be beneficial for simulations of +this kind. The :doc:`neigh_modify ` command has additional options that control how often neighbor lists are built and which pairs are