diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index e52aac10db..9723531625 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -24,6 +24,7 @@ #include "update.h" #include "force.h" #include "fix.h" +#include "fix_neigh_history.h" #include "neighbor.h" #include "neigh_list.h" #include "comm.h" @@ -95,8 +96,8 @@ void PairGranHertzHistory::compute(int eflag, int vflag) ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - firsttouch = list->listhistory->firstneigh; - firstshear = list->listhistory->firstdouble; + firsttouch = fix_history->firstflag; + firstshear = fix_history->firstvalue; // loop over neighbors of my atoms @@ -407,7 +408,7 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype, int jnum = list->numneigh[i]; int *jlist = list->firstneigh[i]; - double *allshear = list->listhistory->firstdouble[i]; + double *allshear = fix_history->firstvalue[i]; for (int jj = 0; jj < jnum; jj++) { neighprev++; diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index e9662c9e73..4f120150de 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -27,7 +27,7 @@ #include "update.h" #include "modify.h" #include "fix.h" -#include "fix_shear_history.h" +#include "fix_neigh_history.h" #include "comm.h" #include "neighbor.h" #include "neigh_list.h" @@ -64,7 +64,7 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp) PairGranHookeHistory::~PairGranHookeHistory() { delete [] svector; - if (fix_history) modify->delete_fix("SHEAR_HISTORY"); + if (fix_history) modify->delete_fix("NEIGH_HISTORY"); if (allocated) { memory->destroy(setflag); @@ -137,8 +137,8 @@ void PairGranHookeHistory::compute(int eflag, int vflag) ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; + firsttouch = fix_history->firstflag; + firstshear = fix_history->firstvalue; // loop over neighbors of my atoms @@ -400,35 +400,28 @@ void PairGranHookeHistory::init_style() if (comm->ghost_velocity == 0) error->all(FLERR,"Pair granular requires ghost atoms store velocity"); - // need a granular neigh list and optionally a granular history neigh list + // need a granular neigh list int irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->size = 1; - if (history) { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->history = 1; - neighbor->requests[irequest]->dnum = 3; - } + 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 *) "SHEAR_HISTORY"; + fixarg[0] = (char *) "NEIGH_HISTORY"; fixarg[1] = (char *) "all"; - fixarg[2] = (char *) "SHEAR_HISTORY"; + fixarg[2] = (char *) "NEIGH_HISTORY"; fixarg[3] = dnumstr; - modify->add_fix(4,fixarg); + modify->add_fix(4,fixarg,1); delete [] fixarg; - fix_history = (FixShearHistory *) modify->fix[modify->nfix-1]; + fix_history = (FixNeighHistory *) modify->fix[modify->nfix-1]; fix_history->pair = this; - neighbor->requests[irequest]->fix_history = fix_history; } // check for FixFreeze and set freeze_group_bit @@ -494,23 +487,12 @@ void PairGranHookeHistory::init_style() // set fix which stores history info if (history) { - int ifix = modify->find_fix("SHEAR_HISTORY"); - if (ifix < 0) error->all(FLERR,"Could not find pair fix ID"); - fix_history = (FixShearHistory *) modify->fix[ifix]; + 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]; } } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - optional granular history list -------------------------------------------------------------------------- */ - -void PairGranHookeHistory::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listhistory = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ @@ -704,7 +686,7 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype, int jnum = list->numneigh[i]; int *jlist = list->firstneigh[i]; - double *allshear = list->listhistory->firstdouble[i]; + double *allshear = fix_history->firstvalue[i]; for (int jj = 0; jj < jnum; jj++) { neighprev++; @@ -797,14 +779,3 @@ double PairGranHookeHistory::memory_usage() double bytes = nmax * sizeof(double); return bytes; } - -/* ---------------------------------------------------------------------- - return ptr to FixShearHistory class - called by Neighbor when setting up neighbor lists -------------------------------------------------------------------------- */ - -void *PairGranHookeHistory::extract(const char *str, int &dim) -{ - if (strcmp(str,"history") == 0) return (void *) fix_history; - return NULL; -} diff --git a/src/GRANULAR/pair_gran_hooke_history.h b/src/GRANULAR/pair_gran_hooke_history.h index afeab93413..f02cccd55e 100644 --- a/src/GRANULAR/pair_gran_hooke_history.h +++ b/src/GRANULAR/pair_gran_hooke_history.h @@ -32,7 +32,6 @@ class PairGranHookeHistory : public Pair { virtual void settings(int, char **); void coeff(int, char **); void init_style(); - void init_list(int, class NeighList *); double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); @@ -43,7 +42,6 @@ class PairGranHookeHistory : public Pair { int pack_forward_comm(int, int *, double *, int, int *); void unpack_forward_comm(int, int, double *); double memory_usage(); - void *extract(const char *, int &); protected: double kn,kt,gamman,gammat,xmu; @@ -56,7 +54,7 @@ class PairGranHookeHistory : public Pair { double *onerad_dynamic,*onerad_frozen; double *maxrad_dynamic,*maxrad_frozen; - class FixShearHistory *fix_history; + class FixNeighHistory *fix_history; // storage of rigid body masses for use in granular interactions diff --git a/src/fix.h b/src/fix.h index 3f32895309..21dfc955a8 100644 --- a/src/fix.h +++ b/src/fix.h @@ -113,6 +113,7 @@ class Fix : protected Pointers { virtual void setup(int) {} virtual void setup_pre_exchange() {} virtual void setup_pre_neighbor() {} + virtual void setup_post_neighbor() {} virtual void setup_pre_force(int) {} virtual void setup_pre_reverse(int, int) {} virtual void min_setup(int) {} @@ -120,6 +121,7 @@ class Fix : protected Pointers { virtual void post_integrate() {} virtual void pre_exchange() {} virtual void pre_neighbor() {} + virtual void post_neighbor() {} virtual void pre_force(int) {} virtual void pre_reverse(int,int) {} virtual void post_force(int) {} @@ -155,6 +157,7 @@ class Fix : protected Pointers { virtual void min_pre_exchange() {} virtual void min_pre_neighbor() {} + virtual void min_post_neighbor() {} virtual void min_pre_force(int) {} virtual void min_pre_reverse(int,int) {} virtual void min_post_force(int) {} @@ -244,25 +247,27 @@ namespace FixConst { static const int POST_INTEGRATE = 1<<1; static const int PRE_EXCHANGE = 1<<2; static const int PRE_NEIGHBOR = 1<<3; - static const int PRE_FORCE = 1<<4; - static const int PRE_REVERSE = 1<<5; - static const int POST_FORCE = 1<<6; - static const int FINAL_INTEGRATE = 1<<7; - static const int END_OF_STEP = 1<<8; - static const int POST_RUN = 1<<9; - static const int THERMO_ENERGY = 1<<10; - static const int INITIAL_INTEGRATE_RESPA = 1<<11; - static const int POST_INTEGRATE_RESPA = 1<<12; - static const int PRE_FORCE_RESPA = 1<<13; - static const int POST_FORCE_RESPA = 1<<14; - static const int FINAL_INTEGRATE_RESPA = 1<<15; - static const int MIN_PRE_EXCHANGE = 1<<16; - static const int MIN_PRE_NEIGHBOR = 1<<17; - static const int MIN_PRE_FORCE = 1<<18; - static const int MIN_PRE_REVERSE = 1<<19; - static const int MIN_POST_FORCE = 1<<20; - static const int MIN_ENERGY = 1<<21; - static const int FIX_CONST_LAST = 1<<22; + 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; } } diff --git a/src/fix_shear_history.cpp b/src/fix_neigh_history.cpp similarity index 55% rename from src/fix_shear_history.cpp rename to src/fix_neigh_history.cpp index 17e78830f4..de94fbecef 100644 --- a/src/fix_shear_history.cpp +++ b/src/fix_neigh_history.cpp @@ -14,7 +14,7 @@ #include #include #include -#include "fix_shear_history.h" +#include "fix_neigh_history.h" #include "atom.h" #include "comm.h" #include "neighbor.h" @@ -33,12 +33,12 @@ enum{DEFAULT,NPARTNER,PERPARTNER}; /* ---------------------------------------------------------------------- */ -FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : +FixNeighHistory::FixNeighHistory(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - npartner(NULL), partner(NULL), shearpartner(NULL), pair(NULL), - ipage(NULL), dpage(NULL) + npartner(NULL), partner(NULL), valuepartner(NULL), pair(NULL), + ipage_atom(NULL), dpage_atom(NULL), ipage_neigh(NULL), dpage_neigh(NULL) { - if (narg != 4) error->all(FLERR,"Illegal fix SHEAR_HISTORY command"); + if (narg != 4) error->all(FLERR,"Illegal fix NEIGH_HISTORY command"); restart_peratom = 1; create_attribute = 1; @@ -48,9 +48,12 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : dnum = force->inumeric(FLERR,arg[3]); dnumbytes = dnum * sizeof(double); + zeroes = new double[dnum]; + for (int i = 0; i < dnum; i++) zeroes[i] = 0.0; + onesided = 0; - if (strcmp(id,"LINE_SHEAR_HISTORY") == 0) onesided = 1; - if (strcmp(id,"TRI_SHEAR_HISTORY") == 0) onesided = 1; + if (strcmp(id,"LINE_NEIGH_HISTORY") == 0) onesided = 1; + if (strcmp(id,"TRI_NEIGH_HISTORY") == 0) onesided = 1; if (newton_pair) comm_reverse = 1; // just for single npartner value // variable-size history communicated via @@ -65,11 +68,24 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : pgsize = oneatom = 0; + // other per-atom vectors + + firstflag = NULL; + firstvalue = NULL; + maxatom = 0; + + // per-atom and per-neighbor data structs + + ipage_atom = NULL; + dpage_atom = NULL; + ipage_neigh = NULL; + dpage_neigh = NULL; + // initialize npartner to 0 so neighbor list creation is OK the 1st time int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) npartner[i] = 0; - maxtouch = 0; + maxpartner = 0; nlocal_neigh = nall_neigh = 0; commflag = DEFAULT; @@ -77,7 +93,7 @@ FixShearHistory::FixShearHistory(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -FixShearHistory::~FixShearHistory() +FixNeighHistory::~FixNeighHistory() { // unregister this fix so atom class doesn't invoke it any more @@ -86,86 +102,111 @@ FixShearHistory::~FixShearHistory() // delete locally stored arrays + delete [] zeroes; + + memory->sfree(firstflag); + memory->sfree(firstvalue); + memory->destroy(npartner); memory->sfree(partner); - memory->sfree(shearpartner); + memory->sfree(valuepartner); + + delete [] ipage_atom; + delete [] dpage_atom; + delete [] ipage_neigh; + delete [] dpage_neigh; // to better detect use-after-delete errors + firstflag = NULL; + firstvalue = NULL; + pair = NULL; npartner = NULL; partner = NULL; - shearpartner = NULL; - - delete [] ipage; - delete [] dpage; + valuepartner = NULL; } /* ---------------------------------------------------------------------- */ -int FixShearHistory::setmask() +int FixNeighHistory::setmask() { int mask = 0; mask |= PRE_EXCHANGE; mask |= MIN_PRE_EXCHANGE; + mask |= POST_NEIGHBOR; + mask |= MIN_POST_NEIGHBOR; return mask; } /* ---------------------------------------------------------------------- */ -void FixShearHistory::init() +void FixNeighHistory::init() { if (atom->tag_enable == 0) - error->all(FLERR,"Granular shear history requires atoms have IDs"); + error->all(FLERR,"Neighbor history requires atoms have IDs"); allocate_pages(); } /* ---------------------------------------------------------------------- create pages if first time or if neighbor pgsize/oneatom has changed - note that latter could cause shear history info to be discarded + note that latter could cause neighbor history info to be discarded ------------------------------------------------------------------------- */ -void FixShearHistory::allocate_pages() +void FixNeighHistory::allocate_pages() { int create = 0; - if (ipage == NULL) create = 1; + if (ipage_atom == NULL) create = 1; if (pgsize != neighbor->pgsize) create = 1; if (oneatom != neighbor->oneatom) create = 1; if (create) { - delete [] ipage; - delete [] dpage; + delete [] ipage_atom; + delete [] dpage_atom; + delete [] ipage_neigh; + delete [] dpage_neigh; pgsize = neighbor->pgsize; oneatom = neighbor->oneatom; int nmypage = comm->nthreads; - ipage = new MyPage[nmypage]; - dpage = new MyPage[nmypage]; + ipage_atom = new MyPage[nmypage]; + dpage_atom = new MyPage[nmypage]; + ipage_neigh = new MyPage[nmypage]; + dpage_neigh = new MyPage[nmypage]; for (int i = 0; i < nmypage; i++) { - ipage[i].init(oneatom,pgsize); - dpage[i].init(dnum*oneatom,dnum*pgsize); + ipage_atom[i].init(oneatom,pgsize); + dpage_atom[i].init(dnum*oneatom,dnum*pgsize); + ipage_neigh[i].init(oneatom,pgsize); + dpage_neigh[i].init(dnum*oneatom,dnum*pgsize); } } } +/* ---------------------------------------------------------------------- */ + +void FixNeighHistory::setup_post_neighbor() +{ + post_neighbor(); +} + /* ---------------------------------------------------------------------- - copy shear partner info from neighbor lists to atom arrays - should be called whenever neighbor list stores current history info - and need to store the info with owned atoms - e.g. so atoms can migrate to new procs or between runs - when atoms may be added or deleted (neighbor list becomes out-of-date) - the next granular neigh list build will put this info back into neigh list + copy partner info from neighbor data structs (NDS) to atom arrays + should be called whenever NDS store current history info + and need to transfer the info to owned atoms + e.g. when atoms migrate to new procs, new neigh list built, or between runs + when atoms may be added or deleted (NDS becomes out-of-date) + the next post_neighbor() will put this info back into new NDS called during run before atom exchanges, including for restart files called at end of run via post_run() do not call during setup of run (setup_pre_exchange) - b/c there is no guarantee of a current neigh list (even on continued run) + b/c there is no guarantee of a current NDS (even on continued run) if run command does a 2nd run with pre = no, then no neigh list will be built, but old neigh list will still have the info onesided and newton on and newton off versions ------------------------------------------------------------------------- */ -void FixShearHistory::pre_exchange() +void FixNeighHistory::pre_exchange() { if (onesided) pre_exchange_onesided(); else if (newton_pair) pre_exchange_newton(); @@ -178,60 +219,57 @@ void FixShearHistory::pre_exchange() only store history info with spheres ------------------------------------------------------------------------- */ -void FixShearHistory::pre_exchange_onesided() +void FixNeighHistory::pre_exchange_onesided() { int i,j,ii,jj,m,n,inum,jnum; int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *shear,*allshear,**firstshear; + int *allflags; + double *allvalues,*onevalues; // NOTE: all operations until very end are with nlocal_neigh <= current nlocal // b/c previous neigh list was built with nlocal_neigh // nlocal can be larger if other fixes added atoms at this pre_exchange() - // zero npartner for owned atoms - // clear 2 page data structures + // clear two paged data structures - for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0; - - ipage->reset(); - dpage->reset(); + ipage_atom->reset(); + dpage_atom->reset(); // 1st loop over neighbor list, I = sphere, J = tri // only calculate npartner for owned spheres + for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0; + tagint *tag = atom->tag; NeighList *list = pair->list; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - firsttouch = list->listhistory->firstneigh; - firstshear = list->listhistory->firstdouble; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; jnum = numneigh[i]; - touch = firsttouch[i]; + allflags = firstflag[i]; for (jj = 0; jj < jnum; jj++) - if (touch[jj]) npartner[i]++; + if (allflags[jj]) npartner[i]++; } - // get page chunks to store atom IDs and shear history for owned atoms + // get page chunks to store partner IDs and values for owned atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; n = npartner[i]; - partner[i] = ipage->get(n); - shearpartner[i] = dpage->get(dnum*n); - if (partner[i] == NULL || shearpartner[i] == NULL) - error->one(FLERR,"Shear history overflow, boost neigh_modify one"); + partner[i] = ipage_atom->get(n); + valuepartner[i] = dpage_atom->get(dnum*n); + if (partner[i] == NULL || valuepartner[i] == NULL) + error->one(FLERR,"Neighbor history overflow, boost neigh_modify one"); } // 2nd loop over neighbor list, I = sphere, J = tri - // store atom IDs and shear history for owned spheres + // store partner IDs and values for owned+ghost atoms // re-zero npartner to use as counter for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0; @@ -239,28 +277,28 @@ void FixShearHistory::pre_exchange_onesided() for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; - allshear = firstshear[i]; jnum = numneigh[i]; - touch = firsttouch[i]; + allflags = firstflag[i]; + allvalues = firstvalue[i]; for (jj = 0; jj < jnum; jj++) { - if (touch[jj]) { - shear = &allshear[dnum*jj]; + if (allflags[jj]) { + onevalues = &allvalues[dnum*jj]; j = jlist[jj]; j &= NEIGHMASK; m = npartner[i]++; partner[i][m] = tag[j]; - memcpy(&shearpartner[i][dnum*m],shear,dnumbytes); + memcpy(&valuepartner[i][dnum*m],onevalues,dnumbytes); } } } - // set maxtouch = max # of partners of any owned atom + // set maxpartner = max # of partners of any owned atom // bump up comm->maxexchange_fix if necessary - maxtouch = 0; - for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]); - comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum+1)*maxtouch+1); + maxpartner = 0; + for (i = 0; i < nlocal_neigh; i++) maxpartner = MAX(maxpartner,npartner[i]); + comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum+1)*maxpartner+1); // zero npartner values from previous nlocal_neigh to current nlocal @@ -269,50 +307,47 @@ void FixShearHistory::pre_exchange_onesided() } /* ---------------------------------------------------------------------- - newton on version, for sphere/sphere contacts - performs reverse comm to acquire shear partner info from ghost atoms + newton ON version + performs reverse comm to acquire partner values from ghost atoms ------------------------------------------------------------------------- */ -void FixShearHistory::pre_exchange_newton() +void FixNeighHistory::pre_exchange_newton() { int i,j,ii,jj,m,n,inum,jnum; int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *shear,*shearj,*allshear,**firstshear; + int *allflags; + double *allvalues,*onevalues,*jvalues; // NOTE: all operations until very end are with // nlocal_neigh <= current nlocal and nall_neigh - // b/c previous neigh list was built with nlocal_neigh,nghost_neigh + // b/c previous neigh list was built with nlocal_neigh & nghost_neigh // nlocal can be larger if other fixes added atoms at this pre_exchange() - // zero npartner for owned+ghost atoms - // clear 2 page data structures + // clear two paged data structures - for (i = 0; i < nall_neigh; i++) npartner[i] = 0; - - ipage->reset(); - dpage->reset(); + ipage_atom->reset(); + dpage_atom->reset(); // 1st loop over neighbor list // calculate npartner for owned+ghost atoms + for (i = 0; i < nall_neigh; i++) npartner[i] = 0; + tagint *tag = atom->tag; NeighList *list = pair->list; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - firsttouch = list->listhistory->firstneigh; - firstshear = list->listhistory->firstdouble; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; jnum = numneigh[i]; - touch = firsttouch[i]; + allflags = firstflag[i]; for (jj = 0; jj < jnum; jj++) { - if (touch[jj]) { + if (allflags[jj]) { npartner[i]++; j = jlist[jj]; j &= NEIGHMASK; @@ -326,29 +361,29 @@ void FixShearHistory::pre_exchange_newton() commflag = NPARTNER; comm->reverse_comm_fix(this,0); - // get page chunks to store atom IDs and shear history for owned+ghost atoms + // get page chunks to store partner IDs and values for owned+ghost atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; n = npartner[i]; - partner[i] = ipage->get(n); - shearpartner[i] = dpage->get(dnum*n); - if (partner[i] == NULL || shearpartner[i] == NULL) { - error->one(FLERR,"Shear history overflow, boost neigh_modify one"); + partner[i] = ipage_atom->get(n); + valuepartner[i] = dpage_atom->get(dnum*n); + if (partner[i] == NULL || valuepartner[i] == NULL) { + error->one(FLERR,"Neighbor history overflow, boost neigh_modify one"); } } for (i = nlocal_neigh; i < nall_neigh; i++) { n = npartner[i]; - partner[i] = ipage->get(n); - shearpartner[i] = dpage->get(dnum*n); - if (partner[i] == NULL || shearpartner[i] == NULL) { - error->one(FLERR,"Shear history overflow, boost neigh_modify one"); + partner[i] = ipage_atom->get(n); + valuepartner[i] = dpage_atom->get(dnum*n); + if (partner[i] == NULL || valuepartner[i] == NULL) { + error->one(FLERR,"Neighbor history overflow, boost neigh_modify one"); } } // 2nd loop over neighbor list - // store atom IDs and shear history for owned+ghost atoms + // store partner IDs and values for owned+ghost atoms // re-zero npartner to use as counter for (i = 0; i < nall_neigh; i++) npartner[i] = 0; @@ -356,40 +391,40 @@ void FixShearHistory::pre_exchange_newton() for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; - allshear = firstshear[i]; jnum = numneigh[i]; - touch = firsttouch[i]; + allflags = firstflag[i]; + allvalues = firstvalue[i]; for (jj = 0; jj < jnum; jj++) { - if (touch[jj]) { - shear = &allshear[dnum*jj]; + if (allflags[jj]) { + onevalues = &allvalues[dnum*jj]; j = jlist[jj]; j &= NEIGHMASK; m = npartner[i]++; partner[i][m] = tag[j]; - memcpy(&shearpartner[i][dnum*m],shear,dnumbytes); + memcpy(&valuepartner[i][dnum*m],onevalues,dnumbytes); m = npartner[j]++; partner[j][m] = tag[i]; - shearj = &shearpartner[j][dnum*m]; - for (n = 0; n < dnum; n++) shearj[n] = -shear[n]; + jvalues = &valuepartner[j][dnum*m]; + for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n]; } } } // perform reverse comm to augment - // owned atom partner/shearpartner with ghost info + // owned atom partner/valuepartner with ghost info // use variable variant b/c size of packed data can be arbitrarily large // if many touching neighbors for large particle commflag = PERPARTNER; comm->reverse_comm_fix_variable(this); - // set maxtouch = max # of partners of any owned atom + // set maxpartner = max # of partners of any owned atom // bump up comm->maxexchange_fix if necessary - maxtouch = 0; - for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]); - comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum+1)*maxtouch+1); + maxpartner = 0; + for (i = 0; i < nlocal_neigh; i++) maxpartner = MAX(maxpartner,npartner[i]); + comm->maxexchange_fix = MAX(comm->maxexchange_fix,4*maxpartner+1); // zero npartner values from previous nlocal_neigh to current nlocal @@ -398,49 +433,47 @@ void FixShearHistory::pre_exchange_newton() } /* ---------------------------------------------------------------------- - newton off version, for sphere/sphere contacts - newton OFF works with smaller vectors that don't include ghost info + newton OFF version + do not need partner values from ghost atoms + assume J values are negative of I values ------------------------------------------------------------------------- */ -void FixShearHistory::pre_exchange_no_newton() +void FixNeighHistory::pre_exchange_no_newton() { int i,j,ii,jj,m,n,inum,jnum; int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *shear,*shearj,*allshear,**firstshear; + int *allflags; + double *allvalues,*onevalues,*jvalues; // NOTE: all operations until very end are with nlocal_neigh <= current nlocal // b/c previous neigh list was built with nlocal_neigh // nlocal can be larger if other fixes added atoms at this pre_exchange() - // zero npartner for owned atoms - // clear 2 page data structures + // clear two paged data structures - for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0; - - ipage->reset(); - dpage->reset(); + ipage_atom->reset(); + dpage_atom->reset(); // 1st loop over neighbor list // calculate npartner for owned atoms + for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0; + tagint *tag = atom->tag; NeighList *list = pair->list; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - firsttouch = list->listhistory->firstneigh; - firstshear = list->listhistory->firstdouble; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; jnum = numneigh[i]; - touch = firsttouch[i]; + allflags = firstflag[i]; for (jj = 0; jj < jnum; jj++) { - if (touch[jj]) { + if (allflags[jj]) { npartner[i]++; j = jlist[jj]; j &= NEIGHMASK; @@ -449,19 +482,19 @@ void FixShearHistory::pre_exchange_no_newton() } } - // get page chunks to store atom IDs and shear history for owned atoms + // get page chunks to store partner IDs and values for owned atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; n = npartner[i]; - partner[i] = ipage->get(n); - shearpartner[i] = dpage->get(dnum*n); - if (partner[i] == NULL || shearpartner[i] == NULL) - error->one(FLERR,"Shear history overflow, boost neigh_modify one"); + partner[i] = ipage_atom->get(n); + valuepartner[i] = dpage_atom->get(dnum*n); + if (partner[i] == NULL || valuepartner[i] == NULL) + error->one(FLERR,"Neighbor history overflow, boost neigh_modify one"); } // 2nd loop over neighbor list - // store atom IDs and shear history for owned atoms + // store partner IDs and values for owned+ghost atoms // re-zero npartner to use as counter for (i = 0; i < nlocal_neigh; i++) npartner[i] = 0; @@ -469,34 +502,34 @@ void FixShearHistory::pre_exchange_no_newton() for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; - allshear = firstshear[i]; jnum = numneigh[i]; - touch = firsttouch[i]; + allflags = firstflag[i]; + allvalues = firstvalue[i]; for (jj = 0; jj < jnum; jj++) { - if (touch[jj]) { - shear = &allshear[dnum*jj]; + if (allflags[jj]) { + onevalues = &allvalues[dnum*jj]; j = jlist[jj]; j &= NEIGHMASK; m = npartner[i]++; partner[i][m] = tag[j]; - memcpy(&shearpartner[i][dnum*m],shear,dnumbytes); + memcpy(&valuepartner[i][dnum*m],onevalues,dnumbytes); if (j < nlocal_neigh) { m = npartner[j]++; partner[j][m] = tag[i]; - shearj = &shearpartner[j][dnum*m]; - for (n = 0; n < dnum; n++) shearj[n] = -shear[n]; + jvalues = &valuepartner[j][dnum*m]; + for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n]; } } } } - // set maxtouch = max # of partners of any owned atom + // set maxpartner = max # of partners of any owned atom // bump up comm->maxexchange_fix if necessary - maxtouch = 0; - for (i = 0; i < nlocal_neigh; i++) maxtouch = MAX(maxtouch,npartner[i]); - comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum+1)*maxtouch+1); + maxpartner = 0; + for (i = 0; i < nlocal_neigh; i++) maxpartner = MAX(maxpartner,npartner[i]); + comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum+1)*maxpartner+1); // zero npartner values from previous nlocal_neigh to current nlocal @@ -506,14 +539,109 @@ void FixShearHistory::pre_exchange_no_newton() /* ---------------------------------------------------------------------- */ -void FixShearHistory::min_pre_exchange() +void FixNeighHistory::min_pre_exchange() { pre_exchange(); } +/* ---------------------------------------------------------------------- + called after neighbor list is build + recover history info stored temporarily in per-atom partner lists + and store afresh in per-neighbor firstflag and firstvalue lists +------------------------------------------------------------------------- */ + +void FixNeighHistory::post_neighbor() +{ + int i,j,m,ii,jj,nn,np,inum,jnum,rflag; + tagint jtag; + double xtmp,ytmp,ztmp,delx,dely,delz; + double radi,rsq,radsum; + int *ilist,*jlist,*numneigh,**firstneigh; + int *allflags; + double *allvalues; + + // store atom counts used for new neighbor list which was just built + + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + nlocal_neigh = nlocal; + nall_neigh = nall; + + // realloc firstflag and firstvalue if needed + + if (maxatom < nlocal) { + memory->sfree(firstflag); + memory->sfree(firstvalue); + maxatom = nall; + firstflag = (int **) + memory->smalloc(maxatom*sizeof(int *),"neighbor_history:firstflag"); + firstvalue = (double **) + memory->smalloc(maxatom*sizeof(double *),"neighbor_history:firstvalue"); + } + + // loop over newly built neighbor list + // repopulate entire per-neighbor data structs + // whether with old-neigh partner info or zeroes + + ipage_neigh->reset(); + dpage_neigh->reset(); + + tagint *tag = atom->tag; + NeighList *list = pair->list; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + firstflag[i] = allflags = ipage_neigh->get(jnum); + firstvalue[i] = allvalues = dpage_neigh->get(jnum*dnum); + np = npartner[i]; + nn = 0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + rflag = sbmask(j); + j &= NEIGHMASK; + jlist[jj] = j; + + // rflag = 1 if r < radsum in npair_size() method + // preserve neigh history info if tag[j] is in old-neigh partner list + // this test could be more geometrically precise for two sphere/line/tri + + if (rflag) { + jtag = tag[j]; + for (m = 0; m < np; m++) + if (partner[i][m] == jtag) break; + if (m < np) { + allflags[jj] = 1; + memcpy(&allvalues[nn],&valuepartner[i][dnum*m],dnumbytes); + } else { + allflags[jj] = 0; + memcpy(&allvalues[nn],zeroes,dnumbytes); + } + } else { + allflags[jj] = 0; + memcpy(&allvalues[nn],zeroes,dnumbytes); + } + nn += dnum; + } + } +} + /* ---------------------------------------------------------------------- */ -void FixShearHistory::post_run() +void FixNeighHistory::min_post_neighbor() +{ + post_neighbor(); +} + +/* ---------------------------------------------------------------------- */ + +void FixNeighHistory::post_run() { pre_exchange(); } @@ -522,17 +650,21 @@ void FixShearHistory::post_run() memory usage of local atom-based arrays ------------------------------------------------------------------------- */ -double FixShearHistory::memory_usage() +double FixNeighHistory::memory_usage() { int nmax = atom->nmax; - double bytes = nmax * sizeof(int); - bytes += nmax * sizeof(tagint *); - bytes += nmax * sizeof(double *); + double bytes = nmax * sizeof(int); // npartner + bytes += nmax * sizeof(tagint *); // partner + bytes += nmax * sizeof(double *); // valuepartner + bytes += maxatom * sizeof(int *); // firstflag + bytes += maxatom * sizeof(double *); // firstvalue int nmypage = comm->nthreads; for (int i = 0; i < nmypage; i++) { - bytes += ipage[i].size(); - bytes += dpage[i].size(); + bytes += ipage_atom[i].size(); + bytes += dpage_atom[i].size(); + bytes += ipage_neigh[i].size(); + bytes += dpage_neigh[i].size(); } return bytes; @@ -542,38 +674,38 @@ double FixShearHistory::memory_usage() allocate local atom-based arrays ------------------------------------------------------------------------- */ -void FixShearHistory::grow_arrays(int nmax) +void FixNeighHistory::grow_arrays(int nmax) { - memory->grow(npartner,nmax,"shear_history:npartner"); + memory->grow(npartner,nmax,"neighbor_history:npartner"); partner = (tagint **) memory->srealloc(partner,nmax*sizeof(tagint *), - "shear_history:partner"); - shearpartner = (double **) memory->srealloc(shearpartner, + "neighbor_history:partner"); + valuepartner = (double **) memory->srealloc(valuepartner, nmax*sizeof(double *), - "shear_history:shearpartner"); + "neighbor_history:valuepartner"); } /* ---------------------------------------------------------------------- copy values within local atom-based arrays ------------------------------------------------------------------------- */ -void FixShearHistory::copy_arrays(int i, int j, int delflag) +void FixNeighHistory::copy_arrays(int i, int j, int delflag) { - // just copy pointers for partner and shearpartner - // b/c can't overwrite chunk allocation inside ipage,dpage + // just copy pointers for partner and valuepartner + // b/c can't overwrite chunk allocation inside ipage_atom,dpage_atom // incoming atoms in unpack_exchange just grab new chunks // so are orphaning chunks for migrating atoms - // OK, b/c will reset ipage,dpage on next reneighboring + // OK, b/c will reset ipage_atom,dpage_atom on next reneighboring npartner[j] = npartner[i]; partner[j] = partner[i]; - shearpartner[j] = shearpartner[i]; + valuepartner[j] = valuepartner[i]; } /* ---------------------------------------------------------------------- initialize one atom's array values, called when atom is created ------------------------------------------------------------------------- */ -void FixShearHistory::set_arrays(int i) +void FixNeighHistory::set_arrays(int i) { npartner[i] = 0; } @@ -582,7 +714,7 @@ void FixShearHistory::set_arrays(int i) only called by Comm::reverse_comm_fix_variable for PERPARTNER mode ------------------------------------------------------------------------- */ -int FixShearHistory::pack_reverse_comm_size(int n, int first) +int FixNeighHistory::pack_reverse_comm_size(int n, int first) { int i,last; @@ -590,7 +722,7 @@ int FixShearHistory::pack_reverse_comm_size(int n, int first) last = first + n; for (i = first; i < last; i++) - m += 1 + (dnum+1)*npartner[i]; + m += 1 + 4*npartner[i]; return m; } @@ -599,7 +731,7 @@ int FixShearHistory::pack_reverse_comm_size(int n, int first) two modes: NPARTNER and PERPARTNER ------------------------------------------------------------------------- */ -int FixShearHistory::pack_reverse_comm(int n, int first, double *buf) +int FixNeighHistory::pack_reverse_comm(int n, int first, double *buf) { int i,k,last; @@ -615,11 +747,11 @@ int FixShearHistory::pack_reverse_comm(int n, int first, double *buf) buf[m++] = npartner[i]; for (k = 0; k < npartner[i]; k++) { buf[m++] = partner[i][k]; - memcpy(&buf[m],&shearpartner[i][dnum*k],dnumbytes); + memcpy(&buf[m],&valuepartner[i][dnum*k],dnumbytes); m += dnum; } } - } else error->all(FLERR,"Unsupported comm mode in shear history"); + } else error->all(FLERR,"Unsupported comm mode in neighbor history"); return m; } @@ -628,7 +760,7 @@ int FixShearHistory::pack_reverse_comm(int n, int first, double *buf) two modes: NPARTNER and PERPARTNER ------------------------------------------------------------------------- */ -void FixShearHistory::unpack_reverse_comm(int n, int *list, double *buf) +void FixNeighHistory::unpack_reverse_comm(int n, int *list, double *buf) { int i,j,k,kk,ncount; @@ -646,18 +778,18 @@ void FixShearHistory::unpack_reverse_comm(int n, int *list, double *buf) for (k = 0; k < ncount; k++) { kk = npartner[j]++; partner[j][kk] = static_cast (buf[m++]); - memcpy(&shearpartner[j][dnum*kk],&buf[m],dnumbytes); + memcpy(&valuepartner[j][dnum*kk],&buf[m],dnumbytes); m += dnum; } } - } else error->all(FLERR,"Unsupported comm mode in shear history"); + } else error->all(FLERR,"Unsupported comm mode in neighbor history"); } /* ---------------------------------------------------------------------- pack values in local atom-based arrays for exchange with another proc ------------------------------------------------------------------------- */ -int FixShearHistory::pack_exchange(int i, double *buf) +int FixNeighHistory::pack_exchange(int i, double *buf) { // NOTE: how do I know comm buf is big enough if extreme # of touching neighs // Comm::BUFEXTRA may need to be increased @@ -666,7 +798,7 @@ int FixShearHistory::pack_exchange(int i, double *buf) buf[m++] = npartner[i]; for (int n = 0; n < npartner[i]; n++) { buf[m++] = partner[i][n]; - memcpy(&buf[m],&shearpartner[i][dnum*n],dnumbytes); + memcpy(&buf[m],&valuepartner[i][dnum*n],dnumbytes); m += dnum; } return m; @@ -676,18 +808,18 @@ int FixShearHistory::pack_exchange(int i, double *buf) unpack values in local atom-based arrays from exchange with another proc ------------------------------------------------------------------------- */ -int FixShearHistory::unpack_exchange(int nlocal, double *buf) +int FixNeighHistory::unpack_exchange(int nlocal, double *buf) { - // allocate new chunks from ipage,dpage for incoming values + // allocate new chunks from ipage_atom,dpage_atom for incoming values int m = 0; npartner[nlocal] = static_cast (buf[m++]); - maxtouch = MAX(maxtouch,npartner[nlocal]); - partner[nlocal] = ipage->get(npartner[nlocal]); - shearpartner[nlocal] = dpage->get(dnum*npartner[nlocal]); + maxpartner = MAX(maxpartner,npartner[nlocal]); + partner[nlocal] = ipage_atom->get(npartner[nlocal]); + valuepartner[nlocal] = dpage_atom->get(dnum*npartner[nlocal]); for (int n = 0; n < npartner[nlocal]; n++) { partner[nlocal][n] = static_cast (buf[m++]); - memcpy(&shearpartner[nlocal][dnum*n],&buf[m],dnumbytes); + memcpy(&valuepartner[nlocal][dnum*n],&buf[m],dnumbytes); m += dnum; } return m; @@ -697,13 +829,13 @@ int FixShearHistory::unpack_exchange(int nlocal, double *buf) pack values in local atom-based arrays for restart file ------------------------------------------------------------------------- */ -int FixShearHistory::pack_restart(int i, double *buf) +int FixNeighHistory::pack_restart(int i, double *buf) { int m = 1; buf[m++] = npartner[i]; for (int n = 0; n < npartner[i]; n++) { buf[m++] = partner[i][n]; - memcpy(&buf[m],&shearpartner[i][dnum*n],dnumbytes); + memcpy(&buf[m],&valuepartner[i][dnum*n],dnumbytes); m += dnum; } buf[0] = m; @@ -714,11 +846,11 @@ int FixShearHistory::pack_restart(int i, double *buf) unpack values from atom->extra array to restart the fix ------------------------------------------------------------------------- */ -void FixShearHistory::unpack_restart(int nlocal, int nth) +void FixNeighHistory::unpack_restart(int nlocal, int nth) { - // ipage = NULL if being called from granular pair style init() + // ipage_atom = NULL if being called from granular pair style init() - if (ipage == NULL) allocate_pages(); + if (ipage_atom == NULL) allocate_pages(); // skip to Nth set of extra values @@ -728,15 +860,15 @@ void FixShearHistory::unpack_restart(int nlocal, int nth) for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); m++; - // allocate new chunks from ipage,dpage for incoming values + // allocate new chunks from ipage_atom,dpage_atom for incoming values npartner[nlocal] = static_cast (extra[nlocal][m++]); - maxtouch = MAX(maxtouch,npartner[nlocal]); - partner[nlocal] = ipage->get(npartner[nlocal]); - shearpartner[nlocal] = dpage->get(dnum*npartner[nlocal]); + maxpartner = MAX(maxpartner,npartner[nlocal]); + partner[nlocal] = ipage_atom->get(npartner[nlocal]); + valuepartner[nlocal] = dpage_atom->get(dnum*npartner[nlocal]); for (int n = 0; n < npartner[nlocal]; n++) { partner[nlocal][n] = static_cast (extra[nlocal][m++]); - memcpy(&shearpartner[nlocal][dnum*n],&extra[nlocal][m],dnumbytes); + memcpy(&valuepartner[nlocal][dnum*n],&extra[nlocal][m],dnumbytes); m += dnum; } } @@ -745,20 +877,20 @@ void FixShearHistory::unpack_restart(int nlocal, int nth) maxsize of any atom's restart data ------------------------------------------------------------------------- */ -int FixShearHistory::maxsize_restart() +int FixNeighHistory::maxsize_restart() { - // maxtouch_all = max # of touching partners across all procs + // maxpartner_all = max # of touching partners across all procs - int maxtouch_all; - MPI_Allreduce(&maxtouch,&maxtouch_all,1,MPI_INT,MPI_MAX,world); - return (dnum+1)*maxtouch_all + 2; + int maxpartner_all; + MPI_Allreduce(&maxpartner,&maxpartner_all,1,MPI_INT,MPI_MAX,world); + return (dnum+1)*maxpartner_all + 2; } /* ---------------------------------------------------------------------- size of atom nlocal's restart data ------------------------------------------------------------------------- */ -int FixShearHistory::size_restart(int nlocal) +int FixNeighHistory::size_restart(int nlocal) { return (dnum+1)*npartner[nlocal] + 2; } diff --git a/src/fix_shear_history.h b/src/fix_neigh_history.h similarity index 62% rename from src/fix_shear_history.h rename to src/fix_neigh_history.h index 00f219f034..64cc11efec 100644 --- a/src/fix_shear_history.h +++ b/src/fix_neigh_history.h @@ -13,38 +13,35 @@ #ifdef FIX_CLASS -FixStyle(SHEAR_HISTORY,FixShearHistory) +FixStyle(NEIGH_HISTORY,FixNeighHistory) #else -#ifndef LMP_FIX_SHEAR_HISTORY_H -#define LMP_FIX_SHEAR_HISTORY_H +#ifndef LMP_FIX_NEIGH_HISTORY_H +#define LMP_FIX_NEIGH_HISTORY_H #include "fix.h" #include "my_page.h" namespace LAMMPS_NS { -class FixShearHistory : public Fix { - //friend class Neighbor; - //friend class PairGranHookeHistory; - friend class PairLineGranHookeHistory; - friend class PairTriGranHookeHistory; - +class FixNeighHistory : public Fix { public: int nlocal_neigh; // nlocal at last time neigh list was built int nall_neigh; // ditto for nlocal+nghost - int *npartner; // # of touching partners of each atom - tagint **partner; // global atom IDs for the partners - double **shearpartner; // shear values with the partner - class Pair *pair; // ptr to pair style that uses shear history + int **firstflag; // ptr to each atom's neighbor flsg + double **firstvalue; // ptr to each atom's values + class Pair *pair; // ptr to pair style that uses neighbor history - FixShearHistory(class LAMMPS *, int, char **); - ~FixShearHistory(); + FixNeighHistory(class LAMMPS *, int, char **); + ~FixNeighHistory(); int setmask(); void init(); + void setup_post_neighbor(); virtual void pre_exchange(); void min_pre_exchange(); + void post_neighbor(); + void min_post_neighbor(); void post_run(); double memory_usage(); @@ -64,20 +61,40 @@ class FixShearHistory : public Fix { protected: int newton_pair; // same as force setting - int dnum,dnumbytes; // dnum = # of shear history values + int dnum,dnumbytes; // dnum = # of values per neighbor int onesided; // 1 for line/tri history, else 0 - int maxtouch; // max # of touching partners for my atoms + int maxatom; // max size of firstflag and firstvalue int commflag; // mode of reverse comm to get ghost info + double *zeroes; + + // per-atom data structures + // partners = flagged neighbors of an atom + + int *npartner; // # of partners of each atom + tagint **partner; // global atom IDs for the partners + double **valuepartner; // values for the partners + int maxpartner; // max # of partners for any of my atoms + + // per-atom data structs pointed to by partner & valuepartner int pgsize,oneatom; // copy of settings in Neighbor - MyPage *ipage; // pages of partner atom IDs - MyPage *dpage; // pages of shear history with partners + MyPage *ipage_atom; // pages of partner atom IDs + MyPage *dpage_atom; // pages of partner values + + // per-neighbor data structs pointed to by firstflag & firstvalue + + MyPage *ipage_neigh; // pages of local atom indices + MyPage *dpage_neigh; // pages of partner values void pre_exchange_onesided(); void pre_exchange_newton(); void pre_exchange_no_newton(); void allocate_pages(); + + inline int sbmask(int j) const { + return j >> SBBITS & 3; + } }; } diff --git a/src/modify.cpp b/src/modify.cpp index 04137f110a..64970f2cf9 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -42,7 +42,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) { nfix = maxfix = 0; n_initial_integrate = n_post_integrate = 0; - n_pre_exchange = n_pre_neighbor = 0; + n_pre_exchange = n_pre_neighbor = n_post_neighbor = 0; n_pre_force = n_pre_reverse = n_post_force = 0; n_final_integrate = n_end_of_step = n_thermo_energy = 0; n_thermo_energy_atom = 0; @@ -54,14 +54,14 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) fix = NULL; fmask = NULL; list_initial_integrate = list_post_integrate = NULL; - list_pre_exchange = list_pre_neighbor = NULL; + list_pre_exchange = list_pre_neighbor = list_post_neighbor = NULL; list_pre_force = list_pre_reverse = list_post_force = NULL; list_final_integrate = list_end_of_step = NULL; list_thermo_energy = list_thermo_energy_atom = NULL; list_initial_integrate_respa = list_post_integrate_respa = NULL; list_pre_force_respa = list_post_force_respa = NULL; list_final_integrate_respa = NULL; - list_min_pre_exchange = list_min_pre_neighbor = NULL; + list_min_pre_exchange = list_min_pre_neighbor = list_min_post_neighbor = NULL; list_min_pre_force = list_min_pre_reverse = list_min_post_force = NULL; list_min_energy = NULL; @@ -123,6 +123,7 @@ Modify::~Modify() delete [] list_post_integrate; delete [] list_pre_exchange; delete [] list_pre_neighbor; + delete [] list_post_neighbor; delete [] list_pre_force; delete [] list_pre_reverse; delete [] list_post_force; @@ -137,6 +138,7 @@ Modify::~Modify() delete [] list_final_integrate_respa; delete [] list_min_pre_exchange; delete [] list_min_pre_neighbor; + delete [] list_min_post_neighbor; delete [] list_min_pre_force; delete [] list_min_pre_reverse; delete [] list_min_post_force; @@ -169,6 +171,7 @@ void Modify::init() list_init(POST_INTEGRATE,n_post_integrate,list_post_integrate); list_init(PRE_EXCHANGE,n_pre_exchange,list_pre_exchange); list_init(PRE_NEIGHBOR,n_pre_neighbor,list_pre_neighbor); + list_init(POST_NEIGHBOR,n_post_neighbor,list_post_neighbor); list_init(PRE_FORCE,n_pre_force,list_pre_force); list_init(PRE_REVERSE,n_pre_reverse,list_pre_reverse); list_init(POST_FORCE,n_post_force,list_post_force); @@ -190,6 +193,7 @@ void Modify::init() list_init(MIN_PRE_EXCHANGE,n_min_pre_exchange,list_min_pre_exchange); list_init(MIN_PRE_NEIGHBOR,n_min_pre_neighbor,list_min_pre_neighbor); + list_init(MIN_POST_NEIGHBOR,n_min_post_neighbor,list_min_post_neighbor); list_init(MIN_PRE_FORCE,n_min_pre_force,list_min_pre_force); list_init(MIN_PRE_REVERSE,n_min_pre_reverse,list_min_pre_reverse); list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force); @@ -329,6 +333,21 @@ void Modify::setup_pre_neighbor() fix[list_min_pre_neighbor[i]]->setup_pre_neighbor(); } +/* ---------------------------------------------------------------------- + setup post_neighbor call, only for fixes that define post_neighbor + called from Verlet, RESPA +------------------------------------------------------------------------- */ + +void Modify::setup_post_neighbor() +{ + if (update->whichflag == 1) + for (int i = 0; i < n_post_neighbor; i++) + fix[list_post_neighbor[i]]->setup_post_neighbor(); + else if (update->whichflag == 2) + for (int i = 0; i < n_min_post_neighbor; i++) + fix[list_min_post_neighbor[i]]->setup_post_neighbor(); +} + /* ---------------------------------------------------------------------- setup pre_force call, only for fixes that define pre_force called from Verlet, RESPA, Min @@ -399,6 +418,16 @@ void Modify::pre_neighbor() fix[list_pre_neighbor[i]]->pre_neighbor(); } +/* ---------------------------------------------------------------------- + post_neighbor call, only for relevant fixes +------------------------------------------------------------------------- */ + +void Modify::post_neighbor() +{ + for (int i = 0; i < n_post_neighbor; i++) + fix[list_post_neighbor[i]]->post_neighbor(); +} + /* ---------------------------------------------------------------------- pre_force call, only for relevant fixes ------------------------------------------------------------------------- */ @@ -589,6 +618,16 @@ void Modify::min_pre_neighbor() fix[list_min_pre_neighbor[i]]->min_pre_neighbor(); } +/* ---------------------------------------------------------------------- + minimizer post-neighbor call, only for relevant fixes +------------------------------------------------------------------------- */ + +void Modify::min_post_neighbor() +{ + for (int i = 0; i < n_min_post_neighbor; i++) + fix[list_min_post_neighbor[i]]->min_post_neighbor(); +} + /* ---------------------------------------------------------------------- minimizer pre-force call, only for relevant fixes ------------------------------------------------------------------------- */ diff --git a/src/modify.h b/src/modify.h index 4ec61f6d57..3e20df5aac 100644 --- a/src/modify.h +++ b/src/modify.h @@ -29,12 +29,13 @@ class Modify : protected Pointers { public: int nfix,maxfix; - int n_initial_integrate,n_post_integrate,n_pre_exchange,n_pre_neighbor; + int n_initial_integrate,n_post_integrate,n_pre_exchange; + int n_pre_neighbor,n_post_neighbor; int n_pre_force,n_pre_reverse,n_post_force; int n_final_integrate,n_end_of_step,n_thermo_energy,n_thermo_energy_atom; int n_initial_integrate_respa,n_post_integrate_respa; int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa; - int n_min_pre_exchange,n_min_pre_neighbor; + int n_min_pre_exchange,n_min_pre_neighbor,n_min_post_neighbor; int n_min_pre_force,n_min_pre_reverse,n_min_post_force,n_min_energy; int restart_pbc_any; // 1 if any fix sets restart_pbc @@ -53,12 +54,14 @@ class Modify : protected Pointers { virtual void setup(int); virtual void setup_pre_exchange(); virtual void setup_pre_neighbor(); + virtual void setup_post_neighbor(); virtual void setup_pre_force(int); virtual void setup_pre_reverse(int, int); virtual void initial_integrate(int); virtual void post_integrate(); virtual void pre_exchange(); virtual void pre_neighbor(); + virtual void post_neighbor(); virtual void pre_force(int); virtual void pre_reverse(int,int); virtual void post_force(int); @@ -78,6 +81,7 @@ class Modify : protected Pointers { virtual void min_pre_exchange(); virtual void min_pre_neighbor(); + virtual void min_post_neighbor(); virtual void min_pre_force(int); virtual void min_pre_reverse(int,int); virtual void min_post_force(int); @@ -123,14 +127,14 @@ class Modify : protected Pointers { // lists of fixes to apply at different stages of timestep int *list_initial_integrate,*list_post_integrate; - int *list_pre_exchange,*list_pre_neighbor; + int *list_pre_exchange,*list_pre_neighbor,*list_post_neighbor; int *list_pre_force,*list_pre_reverse,*list_post_force; int *list_final_integrate,*list_end_of_step,*list_thermo_energy; int *list_thermo_energy_atom; int *list_initial_integrate_respa,*list_post_integrate_respa; int *list_pre_force_respa,*list_post_force_respa; int *list_final_integrate_respa; - int *list_min_pre_exchange,*list_min_pre_neighbor; + int *list_min_pre_exchange,*list_min_pre_neighbor,*list_min_post_neighbor; int *list_min_pre_force,*list_min_pre_reverse,*list_min_post_force; int *list_min_energy; diff --git a/src/neigh_list.cpp b/src/neigh_list.cpp index dde544a69f..3e347360c3 100644 --- a/src/neigh_list.cpp +++ b/src/neigh_list.cpp @@ -40,16 +40,15 @@ NeighList::NeighList(LAMMPS *lmp) : Pointers(lmp) ilist = NULL; numneigh = NULL; firstneigh = NULL; - firstdouble = NULL; // defaults, but may be reset by post_constructor() occasional = 0; ghost = 0; ssa = 0; + history = 0; copy = 0; copymode = 0; - dnum = 0; // ptrs @@ -60,9 +59,6 @@ NeighList::NeighList(LAMMPS *lmp) : Pointers(lmp) listskip = NULL; listfull = NULL; - listhistory = NULL; - fix_history = NULL; - respamiddle = 0; listinner = NULL; listmiddle = NULL; @@ -70,7 +66,6 @@ NeighList::NeighList(LAMMPS *lmp) : Pointers(lmp) fix_bond = NULL; ipage = NULL; - dpage = NULL; // Kokkos package @@ -92,10 +87,7 @@ NeighList::~NeighList() memory->destroy(ilist); memory->destroy(numneigh); memory->sfree(firstneigh); - memory->sfree(firstdouble); - delete [] ipage; - delete [] dpage; } delete [] iskip; @@ -108,7 +100,6 @@ NeighList::~NeighList() copy -> set listcopy for list to copy from skip -> set listskip for list to skip from, create copy of itype,ijtype halffull -> set listfull for full list to derive from - history -> set LH and FH ptrs in partner list that uses the history info respaouter -> set listinner/listmiddle for other rRESPA lists bond -> set fix_bond to Fix that made the request ------------------------------------------------------------------------- */ @@ -120,8 +111,8 @@ void NeighList::post_constructor(NeighRequest *nq) occasional = nq->occasional; ghost = nq->ghost; ssa = nq->ssa; + history = nq->history; copy = nq->copy; - dnum = nq->dnum; if (nq->copy) listcopy = neighbor->lists[nq->copylist]; @@ -141,13 +132,6 @@ void NeighList::post_constructor(NeighRequest *nq) if (nq->halffull) listfull = neighbor->lists[nq->halffulllist]; - if (nq->history) { - neighbor->lists[nq->historylist]->listhistory = this; - int tmp; - neighbor->lists[nq->historylist]->fix_history = - (Fix *) ((Pair *) nq->requestor)->extract("history",tmp); - } - if (nq->respaouter) { if (nq->respamiddlelist < 0) { respamiddle = 0; @@ -173,12 +157,6 @@ void NeighList::setup_pages(int pgsize_caller, int oneatom_caller) ipage = new MyPage[nmypage]; for (int i = 0; i < nmypage; i++) ipage[i].init(oneatom,pgsize,PGDELTA); - - if (dnum) { - dpage = new MyPage[nmypage]; - for (int i = 0; i < nmypage; i++) - dpage[i].init(dnum*oneatom,dnum*pgsize,PGDELTA); - } else dpage = NULL; } /* ---------------------------------------------------------------------- @@ -195,11 +173,11 @@ void NeighList::grow(int nlocal, int nall) { // trigger grow() in children before possible return - if (listhistory) listhistory->grow(nlocal,nall); if (listinner) listinner->grow(nlocal,nall); if (listmiddle) listmiddle->grow(nlocal,nall); // skip if data structs are already big enough + if (ssa) { if ((nlocal * 3) + nall <= maxatom) return; } else if (ghost) { @@ -218,11 +196,6 @@ void NeighList::grow(int nlocal, int nall) memory->create(numneigh,maxatom,"neighlist:numneigh"); firstneigh = (int **) memory->smalloc(maxatom*sizeof(int *), "neighlist:firstneigh"); - if (dnum) { - memory->sfree(firstdouble); - firstdouble = (double **) memory->smalloc(maxatom*sizeof(double *), - "neighlist:firstdouble"); - } } /* ---------------------------------------------------------------------- @@ -262,13 +235,11 @@ void NeighList::print_attributes() printf(" %d = kokkos host\n",rq->kokkos_host); printf(" %d = kokkos device\n",rq->kokkos_device); printf(" %d = ssa flag\n",ssa); - printf(" %d = dnum\n",dnum); printf("\n"); printf(" %d = skip flag\n",rq->skip); printf(" %d = off2on\n",rq->off2on); printf(" %d = copy flag\n",rq->copy); printf(" %d = half/full\n",rq->halffull); - printf(" %d = history/partner\n",rq->history_partner); printf("\n"); } @@ -292,12 +263,5 @@ bigint NeighList::memory_usage() bytes += ipage[i].size(); } - if (dnum && dpage) { - for (int i = 0; i < nmypage; i++) { - bytes += maxatom * sizeof(double *); - bytes += dpage[i].size(); - } - } - return bytes; } diff --git a/src/neigh_list.h b/src/neigh_list.h index 4010a68857..c5b685ca39 100644 --- a/src/neigh_list.h +++ b/src/neigh_list.h @@ -34,9 +34,9 @@ class NeighList : protected Pointers { int occasional; // 0 if build every reneighbor, 1 if not int ghost; // 1 if list stores neighbors of ghosts int ssa; // 1 if list stores Shardlow data - int copy; // 1 if this list is (host) copied from another list + int history; // 1 if there is neigh history (FixNeighHist) + int copy; // 1 if this list is copied from another list int copymode; // 1 if this is a Kokkos on-device copy - int dnum; // # of doubles per neighbor, 0 if none // data structs to store neighbor pairs I,J and associated values @@ -45,13 +45,11 @@ class NeighList : protected Pointers { int *ilist; // local indices of I atoms int *numneigh; // # of J neighbors for each I atom int **firstneigh; // ptr to 1st J int value of each I atom - double **firstdouble; // ptr to 1st J double value of each I atom int maxatom; // size of allocated per-atom arrays int pgsize; // size of each page int oneatom; // max size for one atom MyPage *ipage; // pages of neighbor indices - MyPage *dpage; // pages of neighbor doubles, if dnum > 0 // atom types to skip when building list // copied info from corresponding request into realloced vec/array @@ -65,9 +63,6 @@ class NeighList : protected Pointers { NeighList *listskip; // me = skip list, point to list I skip from NeighList *listfull; // me = half list, point to full I derive from - NeighList *listhistory; // list storing neigh history - class Fix *fix_history; // fix that stores history info - int respamiddle; // 1 if this respaouter has middle list NeighList *listinner; // me = respaouter, point to respainner NeighList *listmiddle; // me = respaouter, point to respamiddle diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp index 8d720e766c..d5397dcf7d 100644 --- a/src/neigh_request.cpp +++ b/src/neigh_request.cpp @@ -42,7 +42,7 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp) // default is use newton_pair setting in force // default is no neighbors of ghosts // default is use cutoffs, not size of particles - // default is no additional neighbor history info + // default is no associated neighbor history info in FixNeighHistory // default is no one-sided sphere/surface interactions (when size = 1) // default is neighbors of atoms, not bonds // default is no multilevel rRESPA neighbors @@ -68,8 +68,6 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp) cut = 0; cutoff = 0.0; - dnum = 0; - // skip info, default is no skipping skip = 0; @@ -88,8 +86,6 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp) copylist = -1; halffull = 0; halffulllist = -1; - history_partner = 0; - historylist = -1; respaouterlist = -1; respamiddlelist = -1; respainnerlist = -1; @@ -158,8 +154,6 @@ int NeighRequest::identical(NeighRequest *other) if (copy != other->copy) same = 0; if (cutoff != other->cutoff) same = 0; - if (dnum != other->dnum) same = 0; - if (skip != other->skip) same = 0; if (skip) same = same_skip(other); @@ -226,8 +220,6 @@ void NeighRequest::copy_request(NeighRequest *other, int skipflag) cut = other->cut; cutoff = other->cutoff; - dnum = other->dnum; - iskip = NULL; ijskip = NULL; diff --git a/src/neigh_request.h b/src/neigh_request.h index 70f7783a70..615b60dfe7 100644 --- a/src/neigh_request.h +++ b/src/neigh_request.h @@ -59,7 +59,7 @@ class NeighRequest : protected Pointers { int ghost; // 1 if includes ghost atom neighbors int size; // 1 if pair cutoff set by particle radius - int history; // 1 if stores neighbor history info + int history; // 1 if there is also neigh history info (FixNeighHist) int granonesided; // 1 if one-sided granular list for // sphere/surf interactions int respainner; // 1 if a rRESPA inner list @@ -74,8 +74,6 @@ class NeighRequest : protected Pointers { int cut; // 1 if use a non-standard cutoff length double cutoff; // special cutoff distance for this list - int dnum; // # of extra floating point values stored in list - // flags set by pair hybrid int skip; // 1 if this list skips atom types from another list @@ -100,10 +98,6 @@ class NeighRequest : protected Pointers { int halffull; // 1 if half list computed from another full list int halffulllist; // index of full list to derive half from - int history_partner; // 1 if this list partners with a history list - int historylist; // index of the associated history list - // for history = 1, index of the non-history partner - int respaouterlist; // index of respaouter/middle/inner lists int respamiddlelist; // which this rREPSA list is associated with int respainnerlist; // each rRESPA style list points at the others diff --git a/src/neighbor.cpp b/src/neighbor.cpp index a460be0065..2d3ed0990e 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -133,8 +133,6 @@ pairclass(NULL), pairnames(NULL), pairmasks(NULL) old_pgsize = pgsize; old_oneatom = oneatom; - zeroes = NULL; - binclass = NULL; binnames = NULL; binmasks = NULL; @@ -208,8 +206,6 @@ Neighbor::~Neighbor() if (old_requests[i]) delete old_requests[i]; memory->sfree(old_requests); - delete [] zeroes; - delete [] binclass; delete [] binnames; delete [] binmasks; @@ -666,14 +662,12 @@ int Neighbor::init_pair() // purpose is to avoid duplicate or inefficient builds // may add new requests if a needed request to derive from does not exist // methods: - // (1) other = point history and rRESPA lists at their partner lists + // (1) other = point rRESPA lists at their partner lists // (2) skip = create any new non-skip lists needed by pair hybrid skip lists // (3) granular = adjust parent and skip lists for granular onesided usage // (4) h/f = pair up any matching half/full lists // (5) copy = convert as many lists as possible to copy lists // order of morph methods matters: - // (1) before (2), b/c (2) needs to know history partner pairings - // (2) after (1), b/c (2) may also need to create new history lists // (3) after (2), b/c it adjusts lists created by (2) // (4) after (2) and (3), // b/c (2) may create new full lists, (3) may change them @@ -827,23 +821,14 @@ int Neighbor::init_pair() } // allocate initial pages for each list, except if copy flag set - // allocate dnum vector of zeroes if set - int dnummax = 0; for (i = 0; i < nlist; i++) { if (lists[i]->copy) continue; lists[i]->setup_pages(pgsize,oneatom); - dnummax = MAX(dnummax,lists[i]->dnum); - } - - if (dnummax) { - delete [] zeroes; - zeroes = new double[dnummax]; - for (i = 0; i < dnummax; i++) zeroes[i] = 0.0; } // first-time allocation of per-atom data for lists that are built and store - // lists that are not built: granhistory, respa inner/middle (no neigh_pair) + // lists that are not built: respa inner/middle (no neigh_pair) // lists that do not store: copy // use atom->nmax for both grow() args // i.e. grow first time to expanded size to avoid future reallocs @@ -923,7 +908,7 @@ int Neighbor::init_pair() /* ---------------------------------------------------------------------- scan NeighRequests to set additional flags - only for history, respaouter, custom cutoff lists + only for respaouter, custom cutoff lists ------------------------------------------------------------------------- */ void Neighbor::morph_other() @@ -933,14 +918,6 @@ void Neighbor::morph_other() for (int i = 0; i < nrequest; i++) { irq = requests[i]; - // if history, point this list and partner list at each other - - if (irq->history) { - irq->historylist = i-1; - requests[i-1]->history_partner = 1; - requests[i-1]->historylist = i; - } - // if respaouter, point all associated rRESPA lists at each other if (irq->respaouter) { @@ -987,7 +964,6 @@ void Neighbor::morph_skip() // halffull list and its full parent may both skip, // but are checked to insure matching skip info - if (irq->history) continue; if (irq->respainner || irq->respamiddle) continue; if (irq->halffull) continue; if (irq->copy) continue; @@ -1022,11 +998,11 @@ void Neighbor::morph_skip() // or their data structures are different // this includes custom cutoff set by requestor // no need to check respaouter b/c it stores same pairs - // no need to check dnum b/c only set for history // NOTE: need check for 2 Kokkos flags? if (irq->ghost != jrq->ghost) continue; if (irq->size != jrq->size) continue; + if (irq->history != jrq->history) continue; if (irq->bond != jrq->bond) continue; if (irq->omp != jrq->omp) continue; if (irq->intel != jrq->intel) continue; @@ -1045,8 +1021,8 @@ void Neighbor::morph_skip() // else create a new identical list except non-skip // for new list, set neigh = 1, skip = 0, no skip vec/array, // copy unique flag (since copy_request() will not do it) - // note: parents of skip lists do not have associated history list - // b/c child skip lists store their own history info + // note: parents of skip lists do not have associated history + // b/c child skip lists have the associated history if (j < nrequest) irq->skiplist = j; else { @@ -1107,7 +1083,6 @@ void Neighbor::morph_granular() if (onesided == 2) break; } - // if onesided = 2, parent has children with both granonesided = 0/1 // force parent newton off (newton = 2) to enable onesided skip by child // set parent granonesided = 0, so it stores all neighs in usual manner @@ -1159,7 +1134,6 @@ void Neighbor::morph_halffull() // these lists are created other ways, no need for halffull // do want to process skip lists - if (irq->history) continue; if (irq->respainner || irq->respamiddle) continue; if (irq->copy) continue; @@ -1180,10 +1154,10 @@ void Neighbor::morph_halffull() // or their data structures are different // this includes custom cutoff set by requestor // no need to check respaouter b/c it stores same pairs - // no need to check dnum b/c only set for history if (irq->ghost != jrq->ghost) continue; if (irq->size != jrq->size) continue; + if (irq->history != jrq->history) continue; if (irq->bond != jrq->bond) continue; if (irq->omp != jrq->omp) continue; if (irq->intel != jrq->intel) continue; @@ -1233,7 +1207,6 @@ void Neighbor::morph_copy() // these lists are created other ways, no need to copy // skip lists are eligible to become a copy list - if (irq->history) continue; if (irq->respainner || irq->respamiddle) continue; // check all other lists @@ -1272,9 +1245,8 @@ void Neighbor::morph_copy() if (irq->ghost && !jrq->ghost) continue; - // do not copy from a history list or a respa middle/inner list + // do not copy from a respa middle/inner list - if (jrq->history) continue; if (jrq->respamiddle) continue; if (jrq->respainner) continue; @@ -1284,10 +1256,10 @@ void Neighbor::morph_copy() // this includes custom cutoff set by requestor // no need to check respaouter b/c it stores same pairs // no need to check omp b/c it stores same pairs - // no need to check dnum b/c only set for history // NOTE: need check for 2 Kokkos flags? if (irq->size != jrq->size) continue; + if (irq->history != jrq->history) continue; if (irq->bond != jrq->bond) continue; if (irq->intel != jrq->intel) continue; if (irq->kokkos_host != jrq->kokkos_host) continue; @@ -1535,9 +1507,7 @@ void Neighbor::print_pairwise_info() // order these to get single output of most relevant - if (rq->history) - fprintf(out,", history for (%d)",rq->historylist+1); - else if (rq->copy) + if (rq->copy) fprintf(out,", copy from (%d)",rq->copylist+1); else if (rq->halffull) fprintf(out,", half/full from (%d)",rq->halffulllist+1); @@ -1659,7 +1629,6 @@ int Neighbor::choose_bin(NeighRequest *rq) if (style == NSQ) return 0; if (rq->skip || rq->copy || rq->halffull) return 0; - if (rq->history) return 0; if (rq->respainner || rq->respamiddle) return 0; // use request settings to match exactly one NBin class mask @@ -1701,7 +1670,6 @@ int Neighbor::choose_stencil(NeighRequest *rq) if (style == NSQ) return 0; if (rq->skip || rq->copy || rq->halffull) return 0; - if (rq->history) return 0; if (rq->respainner || rq->respamiddle) return 0; // convert newton request to newtflag = on or off @@ -1795,7 +1763,6 @@ int Neighbor::choose_pair(NeighRequest *rq) { // no neighbor list build performed - if (rq->history) return 0; if (rq->respainner || rq->respamiddle) return 0; // error check for includegroup with ghost neighbor request diff --git a/src/neighbor.h b/src/neighbor.h index 64bced2293..7520e191ec 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -54,7 +54,6 @@ class Neighbor : protected Pointers { double *bboxlo,*bboxhi; // ptrs to full domain bounding box // different for orthog vs triclinic - double *zeroes; // vector of zeroes for shear history init // exclusion info, used by NeighPair diff --git a/src/npair.cpp b/src/npair.cpp index 9fbb4d219d..dd3a73926e 100644 --- a/src/npair.cpp +++ b/src/npair.cpp @@ -66,7 +66,6 @@ void NPair::copy_neighbor_info() cut_inner_sq = neighbor->cut_inner_sq; cut_middle_sq = neighbor->cut_middle_sq; cut_middle_inside_sq = neighbor->cut_middle_inside_sq; - zeroes = neighbor->zeroes; bboxlo = neighbor->bboxlo; bboxhi = neighbor->bboxhi; diff --git a/src/npair.h b/src/npair.h index 6941b86164..289a1348d9 100644 --- a/src/npair.h +++ b/src/npair.h @@ -47,7 +47,6 @@ class NPair : protected Pointers { double cut_inner_sq; double cut_middle_sq; double cut_middle_inside_sq; - double *zeroes; double *bboxlo,*bboxhi; // exclusion data from Neighbor class diff --git a/src/npair_copy.cpp b/src/npair_copy.cpp index 1799d48fed..9426d22ed3 100644 --- a/src/npair_copy.cpp +++ b/src/npair_copy.cpp @@ -40,7 +40,5 @@ void NPairCopy::build(NeighList *list) list->ilist = listcopy->ilist; list->numneigh = listcopy->numneigh; list->firstneigh = listcopy->firstneigh; - list->firstdouble = listcopy->firstdouble; list->ipage = listcopy->ipage; - list->dpage = listcopy->dpage; } diff --git a/src/npair_half_size_bin_newtoff.cpp b/src/npair_half_size_bin_newtoff.cpp index 571b2484ea..cf608b5d59 100644 --- a/src/npair_half_size_bin_newtoff.cpp +++ b/src/npair_half_size_bin_newtoff.cpp @@ -17,9 +17,6 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "domain.h" -#include "fix_shear_history.h" #include "my_page.h" #include "error.h" @@ -32,7 +29,6 @@ NPairHalfSizeBinNewtoff::NPairHalfSizeBinNewtoff(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- size particles binned neighbor list construction with partial Newton's 3rd law - shear history must be accounted for when a neighbor pair is added each owned atom i checks own bin and surrounding bins in non-Newton stencil 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) @@ -40,20 +36,10 @@ NPairHalfSizeBinNewtoff::NPairHalfSizeBinNewtoff(LAMMPS *lmp) : NPair(lmp) {} void NPairHalfSizeBinNewtoff::build(NeighList *list) { - int i,j,k,m,n,nn,ibin,dnum,dnumbytes; + int i,j,k,m,n,nn,ibin; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; - int *neighptr,*touchptr; - double *shearptr; - - int *npartner; - tagint **partner; - double **shearpartner; - int **firsttouch; - double **firstshear; - MyPage *ipage_touch; - MyPage *dpage_shear; - NeighList *listhistory; + int *neighptr; double **x = atom->x; double *radius = atom->radius; @@ -64,42 +50,20 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list) 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; - FixShearHistory *fix_history = (FixShearHistory *) list->fix_history; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nlocal + atom->nghost; - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - listhistory = list->listhistory; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage; - dpage_shear = listhistory->dpage; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - } + int mask_history = 3 << SBBITS; int inum = 0; ipage->reset(); - if (fix_history) { - ipage_touch->reset(); - dpage_shear->reset(); - } for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - if (fix_history) { - nn = 0; - touchptr = ipage_touch->vget(); - shearptr = dpage_shear->vget(); - } xtmp = x[i][0]; ytmp = x[i][1]; @@ -116,38 +80,19 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list) for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { if (j <= i) continue; if (exclude && exclusion(i,j,type[i],type[j],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]; cutsq = (radsum+skin) * (radsum+skin); - + if (rsq <= cutsq) { - neighptr[n] = j; - - if (fix_history) { - if (rsq < radsum*radsum) { - for (m = 0; m < npartner[i]; m++) - if (partner[i][m] == tag[j]) break; - if (m < npartner[i]) { - touchptr[n] = 1; - memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes); - nn += dnum; - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } - - n++; + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; } } } @@ -158,13 +103,6 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list) ipage->vgot(n); if (ipage->status()) error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); - - if (fix_history) { - firsttouch[i] = touchptr; - firstshear[i] = shearptr; - ipage_touch->vgot(n); - dpage_shear->vgot(nn); - } } list->inum = inum; diff --git a/src/npair_half_size_bin_newton.cpp b/src/npair_half_size_bin_newton.cpp index 4f4ecccb16..662bf91d6e 100644 --- a/src/npair_half_size_bin_newton.cpp +++ b/src/npair_half_size_bin_newton.cpp @@ -17,9 +17,6 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "domain.h" -#include "fix_shear_history.h" #include "my_page.h" #include "error.h" @@ -32,27 +29,16 @@ NPairHalfSizeBinNewton::NPairHalfSizeBinNewton(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- size particles binned neighbor list construction with full Newton's 3rd law - shear history must be accounted for when a neighbor pair is added each owned atom i checks its own bin and other bins in Newton stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void NPairHalfSizeBinNewton::build(NeighList *list) { - int i,j,k,m,n,nn,ibin,dnum,dnumbytes; + int i,j,k,m,n,nn,ibin; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; - int *neighptr,*touchptr; - double *shearptr; - - int *npartner; - tagint **partner; - double **shearpartner; - int **firsttouch; - double **firstshear; - MyPage *ipage_touch; - MyPage *dpage_shear; - NeighList *listhistory; + int *neighptr; double **x = atom->x; double *radius = atom->radius; @@ -63,42 +49,20 @@ void NPairHalfSizeBinNewton::build(NeighList *list) 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; - FixShearHistory *fix_history = (FixShearHistory *) list->fix_history; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nlocal + atom->nghost; - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - listhistory = list->listhistory; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage; - dpage_shear = listhistory->dpage; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - } + int mask_history = 3 << SBBITS; int inum = 0; ipage->reset(); - if (fix_history) { - ipage_touch->reset(); - dpage_shear->reset(); - } for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - if (fix_history) { - nn = 0; - touchptr = ipage_touch->vget(); - shearptr = dpage_shear->vget(); - } xtmp = x[i][0]; ytmp = x[i][1]; @@ -128,29 +92,10 @@ void NPairHalfSizeBinNewton::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { - neighptr[n] = j; - - if (fix_history) { - if (rsq < radsum*radsum) { - for (m = 0; m < npartner[i]; m++) - if (partner[i][m] == tag[j]) break; - if (m < npartner[i]) { - touchptr[n] = 1; - memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes); - nn += dnum; - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } - - n++; + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; } } @@ -169,29 +114,10 @@ void NPairHalfSizeBinNewton::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { - neighptr[n] = j; - - if (fix_history) { - if (rsq < radsum*radsum) { - for (m = 0; m < npartner[i]; m++) - if (partner[i][m] == tag[j]) break; - if (m < npartner[i]) { - touchptr[n] = 1; - memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes); - nn += dnum; - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } - - n++; + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; } } } @@ -202,13 +128,6 @@ void NPairHalfSizeBinNewton::build(NeighList *list) ipage->vgot(n); if (ipage->status()) error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); - - if (fix_history) { - firsttouch[i] = touchptr; - firstshear[i] = shearptr; - ipage_touch->vgot(n); - dpage_shear->vgot(nn); - } } list->inum = inum; diff --git a/src/npair_half_size_bin_newton_tri.cpp b/src/npair_half_size_bin_newton_tri.cpp index 559eb09a7a..e70c072280 100644 --- a/src/npair_half_size_bin_newton_tri.cpp +++ b/src/npair_half_size_bin_newton_tri.cpp @@ -17,9 +17,6 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "domain.h" -#include "fix_shear_history.h" #include "my_page.h" #include "error.h" @@ -33,27 +30,16 @@ NPairHalfSizeBinNewtonTri::NPairHalfSizeBinNewtonTri(LAMMPS *lmp) : /* ---------------------------------------------------------------------- size particles binned neighbor list construction with Newton's 3rd law for triclinic - shear history must be accounted for when a neighbor pair is added each owned atom i checks its own bin and other bins in triclinic stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void NPairHalfSizeBinNewtonTri::build(NeighList *list) { - int i,j,k,m,n,nn,ibin,dnum,dnumbytes; + int i,j,k,m,n,nn,ibin; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; - int *neighptr,*touchptr; - double *shearptr; - - int *npartner; - tagint **partner; - double **shearpartner; - int **firsttouch; - double **firstshear; - MyPage *ipage_touch; - MyPage *dpage_shear; - NeighList *listhistory; + int *neighptr; double **x = atom->x; double *radius = atom->radius; @@ -64,42 +50,20 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list) 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; - FixShearHistory *fix_history = (FixShearHistory *) list->fix_history; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nlocal + atom->nghost; - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - listhistory = list->listhistory; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage; - dpage_shear = listhistory->dpage; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - } + int mask_history = 3 << SBBITS; int inum = 0; ipage->reset(); - if (fix_history) { - ipage_touch->reset(); - dpage_shear->reset(); - } for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - if (fix_history) { - nn = 0; - touchptr = ipage_touch->vget(); - shearptr = dpage_shear->vget(); - } xtmp = x[i][0]; ytmp = x[i][1]; @@ -134,29 +98,10 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { - neighptr[n++] = j; - - if (fix_history) { - if (rsq < radsum*radsum) { - for (m = 0; m < npartner[i]; m++) - if (partner[i][m] == tag[j]) break; - if (m < npartner[i]) { - touchptr[n] = 1; - memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes); - nn += dnum; - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } - - n++; + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; } } } @@ -167,13 +112,6 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list) ipage->vgot(n); if (ipage->status()) error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); - - if (fix_history) { - firsttouch[i] = touchptr; - firstshear[i] = shearptr; - ipage_touch->vgot(n); - dpage_shear->vgot(nn); - } } list->inum = inum; diff --git a/src/npair_half_size_nsq_newtoff.cpp b/src/npair_half_size_nsq_newtoff.cpp index 56630a9dc8..e6f5cba657 100644 --- a/src/npair_half_size_nsq_newtoff.cpp +++ b/src/npair_half_size_nsq_newtoff.cpp @@ -18,9 +18,6 @@ #include "atom.h" #include "atom_vec.h" #include "group.h" -#include "molecule.h" -#include "domain.h" -#include "fix_shear_history.h" #include "my_page.h" #include "error.h" @@ -33,27 +30,16 @@ NPairHalfSizeNsqNewtoff::NPairHalfSizeNsqNewtoff(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- size particles N^2 / 2 search for neighbor pairs with partial Newton's 3rd law - shear history must be accounted for when a neighbor pair is added pair added to list if atoms i and j are both owned and i < j pair added if j is ghost (also stored by proc owning j) ------------------------------------------------------------------------- */ void NPairHalfSizeNsqNewtoff::build(NeighList *list) { - int i,j,m,n,nn,bitmask,dnum,dnumbytes; + int i,j,m,n,nn,bitmask; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; - int *neighptr,*touchptr; - double *shearptr; - - int *npartner; - tagint **partner; - double **shearpartner; - int **firsttouch; - double **firstshear; - MyPage *ipage_touch; - MyPage *dpage_shear; - NeighList *listhistory; + int *neighptr; double **x = atom->x; double *radius = atom->radius; @@ -68,42 +54,20 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list) bitmask = group->bitmask[includegroup]; } + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - FixShearHistory *fix_history = (FixShearHistory *) list->fix_history; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nall; - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - listhistory = list->listhistory; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage; - dpage_shear = listhistory->dpage; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - } + int mask_history = 3 << SBBITS; int inum = 0; ipage->reset(); - if (fix_history) { - ipage_touch->reset(); - dpage_shear->reset(); - } for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - if (fix_history) { - nn = 0; - touchptr = ipage_touch->vget(); - shearptr = dpage_shear->vget(); - } xtmp = x[i][0]; ytmp = x[i][1]; @@ -124,29 +88,10 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { - neighptr[n] = j; - - if (fix_history) { - if (rsq < radsum*radsum) { - for (m = 0; m < npartner[i]; m++) - if (partner[i][m] == tag[j]) break; - if (m < npartner[i]) { - touchptr[n] = 1; - memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes); - nn += dnum; - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } - - n++; + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; } } @@ -156,13 +101,6 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list) ipage->vgot(n); if (ipage->status()) error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); - - if (fix_history) { - firsttouch[i] = touchptr; - firstshear[i] = shearptr; - ipage_touch->vgot(n); - dpage_shear->vgot(nn); - } } list->inum = inum; diff --git a/src/npair_half_size_nsq_newton.cpp b/src/npair_half_size_nsq_newton.cpp index 177685b9fc..78811170cb 100644 --- a/src/npair_half_size_nsq_newton.cpp +++ b/src/npair_half_size_nsq_newton.cpp @@ -18,9 +18,6 @@ #include "atom.h" #include "atom_vec.h" #include "group.h" -#include "molecule.h" -#include "domain.h" -#include "fix_shear_history.h" #include "my_page.h" #include "error.h" @@ -33,7 +30,6 @@ NPairHalfSizeNsqNewton::NPairHalfSizeNsqNewton(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- size particles N^2 / 2 search for neighbor pairs with full Newton's 3rd law - shear history must be accounted for when a neighbor pair is added pair added to list if atoms i and j are both owned and i < j if j is ghost only me or other proc adds pair decision based on itag,jtag tests @@ -41,20 +37,10 @@ NPairHalfSizeNsqNewton::NPairHalfSizeNsqNewton(LAMMPS *lmp) : NPair(lmp) {} void NPairHalfSizeNsqNewton::build(NeighList *list) { - int i,j,m,n,nn,itag,jtag,bitmask,dnum,dnumbytes; + int i,j,m,n,nn,itag,jtag,bitmask; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; - int *neighptr,*touchptr; - double *shearptr; - - int *npartner; - tagint **partner; - double **shearpartner; - int **firsttouch; - double **firstshear; - MyPage *ipage_touch; - MyPage *dpage_shear; - NeighList *listhistory; + int *neighptr; double **x = atom->x; double *radius = atom->radius; @@ -69,42 +55,20 @@ void NPairHalfSizeNsqNewton::build(NeighList *list) bitmask = group->bitmask[includegroup]; } + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - FixShearHistory *fix_history = (FixShearHistory *) list->fix_history; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nall; - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - listhistory = list->listhistory; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage; - dpage_shear = listhistory->dpage; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - } + int mask_history = 3 << SBBITS; int inum = 0; ipage->reset(); - if (fix_history) { - ipage_touch->reset(); - dpage_shear->reset(); - } for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - if (fix_history) { - nn = 0; - touchptr = ipage_touch->vget(); - shearptr = dpage_shear->vget(); - } itag = tag[i]; xtmp = x[i][0]; @@ -142,29 +106,10 @@ void NPairHalfSizeNsqNewton::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { - neighptr[n] = j; - - if (fix_history) { - if (rsq < radsum*radsum) { - for (m = 0; m < npartner[i]; m++) - if (partner[i][m] == tag[j]) break; - if (m < npartner[i]) { - touchptr[n] = 1; - memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes); - nn += dnum; - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } - - n++; + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; } } @@ -174,13 +119,6 @@ void NPairHalfSizeNsqNewton::build(NeighList *list) ipage->vgot(n); if (ipage->status()) error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); - - if (fix_history) { - firsttouch[i] = touchptr; - firstshear[i] = shearptr; - ipage_touch->vgot(n); - dpage_shear->vgot(nn); - } } list->inum = inum; diff --git a/src/npair_skip_size.cpp b/src/npair_skip_size.cpp index e8d19dedca..075387f5b0 100644 --- a/src/npair_skip_size.cpp +++ b/src/npair_skip_size.cpp @@ -17,9 +17,6 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "domain.h" -#include "fix_shear_history.h" #include "my_page.h" #include "error.h" @@ -32,24 +29,13 @@ NPairSkipSize::NPairSkipSize(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- build skip list for subset of types from parent list iskip and ijskip flag which atom types and type pairs to skip - if list requests it, preserve shear history via fix shear/history ------------------------------------------------------------------------- */ void NPairSkipSize::build(NeighList *list) { int i,j,ii,jj,m,n,nn,itype,jnum,joriginal,dnum,dnumbytes; tagint jtag; - int *neighptr,*jlist,*touchptr; - double *shearptr; - - int *npartner; - tagint **partner; - double **shearpartner; - int **firsttouch; - double **firstshear; - MyPage *ipage_touch; - MyPage *dpage_shear; - NeighList *listhistory; + int *neighptr,*jlist; tagint *tag = atom->tag; int *type = atom->type; @@ -68,28 +54,8 @@ void NPairSkipSize::build(NeighList *list) int *iskip = list->iskip; int **ijskip = list->ijskip; - FixShearHistory *fix_history = (FixShearHistory *) list->fix_history; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nlocal + atom->nghost; - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - listhistory = list->listhistory; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage; - dpage_shear = listhistory->dpage; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - } - int inum = 0; ipage->reset(); - if (fix_history) { - ipage_touch->reset(); - dpage_shear->reset(); - } // loop over atoms in other list // skip I atom entirely if iskip is set for type[I] @@ -102,13 +68,8 @@ void NPairSkipSize::build(NeighList *list) n = 0; neighptr = ipage->vget(); - if (fix_history) { - nn = 0; - touchptr = ipage_touch->vget(); - shearptr = dpage_shear->vget(); - } - // loop over parent non-skip size list and optionally its history info + // loop over parent non-skip size list jlist = firstneigh_skip[i]; jnum = numneigh_skip[i]; @@ -117,29 +78,7 @@ void NPairSkipSize::build(NeighList *list) joriginal = jlist[jj]; j = joriginal & NEIGHMASK; if (ijskip[itype][type[j]]) continue; - neighptr[n] = joriginal; - - // no numeric test for current touch - // just use FSH partner list to infer it - // would require distance calculation for spheres - // more complex calculation for surfs - - if (fix_history) { - jtag = tag[j]; - for (m = 0; m < npartner[i]; m++) - if (partner[i][m] == jtag) break; - if (m < npartner[i]) { - touchptr[n] = 1; - memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes); - nn += dnum; - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } - - n++; + neighptr[n++] = joriginal; } ilist[inum++] = i; @@ -148,13 +87,6 @@ void NPairSkipSize::build(NeighList *list) ipage->vgot(n); if (ipage->status()) error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); - - if (fix_history) { - firsttouch[i] = touchptr; - firstshear[i] = shearptr; - ipage_touch->vgot(n); - dpage_shear->vgot(nn); - } } list->inum = inum; diff --git a/src/npair_skip_size_off2on.cpp b/src/npair_skip_size_off2on.cpp index da9dd57047..92eae285d0 100644 --- a/src/npair_skip_size_off2on.cpp +++ b/src/npair_skip_size_off2on.cpp @@ -17,9 +17,6 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "domain.h" -#include "fix_shear_history.h" #include "my_page.h" #include "error.h" @@ -33,24 +30,13 @@ NPairSkipSizeOff2on::NPairSkipSizeOff2on(LAMMPS *lmp) : NPair(lmp) {} build skip list for subset of types from parent list iskip and ijskip flag which atom types and type pairs to skip parent non-skip list used newton off, this skip list is newton on - if list requests it, preserve shear history via fix shear/history ------------------------------------------------------------------------- */ void NPairSkipSizeOff2on::build(NeighList *list) { int i,j,ii,jj,m,n,nn,itype,jnum,joriginal,dnum,dnumbytes; tagint itag,jtag; - int *neighptr,*jlist,*touchptr; - double *shearptr; - - int *npartner; - tagint **partner; - double **shearpartner; - int **firsttouch; - double **firstshear; - MyPage *ipage_touch; - MyPage *dpage_shear; - NeighList *listhistory; + int *neighptr,*jlist; tagint *tag = atom->tag; int *type = atom->type; @@ -69,28 +55,8 @@ void NPairSkipSizeOff2on::build(NeighList *list) int *iskip = list->iskip; int **ijskip = list->ijskip; - FixShearHistory *fix_history = (FixShearHistory *) list->fix_history; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nlocal + atom->nghost; - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - listhistory = list->listhistory; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage; - dpage_shear = listhistory->dpage; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - } - int inum = 0; ipage->reset(); - if (fix_history) { - ipage_touch->reset(); - dpage_shear->reset(); - } // loop over atoms in other list // skip I atom entirely if iskip is set for type[I] @@ -104,11 +70,6 @@ void NPairSkipSizeOff2on::build(NeighList *list) n = 0; neighptr = ipage->vget(); - if (fix_history) { - nn = 0; - touchptr = ipage_touch->vget(); - shearptr = dpage_shear->vget(); - } // loop over parent non-skip size list and optionally its history info @@ -125,28 +86,7 @@ void NPairSkipSizeOff2on::build(NeighList *list) jtag = tag[j]; if (j >= nlocal && jtag < itag) continue; - neighptr[n] = joriginal; - - // no numeric test for current touch - // just use FSH partner list to infer it - // would require distance calculation for spheres - // more complex calculation for surfs - - if (fix_history) { - for (m = 0; m < npartner[i]; m++) - if (partner[i][m] == jtag) break; - if (m < npartner[i]) { - touchptr[n] = 1; - memcpy(&shearptr[nn],&shearpartner[i][dnum*m],dnumbytes); - nn += dnum; - } else { - touchptr[n] = 0; - memcpy(&shearptr[nn],zeroes,dnumbytes); - nn += dnum; - } - } - - n++; + neighptr[n++] = joriginal; } ilist[inum++] = i; @@ -155,13 +95,6 @@ void NPairSkipSizeOff2on::build(NeighList *list) ipage->vgot(n); if (ipage->status()) error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); - - if (fix_history) { - firsttouch[i] = touchptr; - firstshear[i] = shearptr; - ipage_touch->vgot(n); - dpage_shear->vgot(nn); - } } list->inum = inum; diff --git a/src/npair_skip_size_off2on_oneside.cpp b/src/npair_skip_size_off2on_oneside.cpp index 7377feec5b..f2fca7b128 100644 --- a/src/npair_skip_size_off2on_oneside.cpp +++ b/src/npair_skip_size_off2on_oneside.cpp @@ -17,9 +17,7 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" -#include "molecule.h" #include "domain.h" -#include "fix_shear_history.h" #include "my_page.h" #include "error.h" @@ -35,7 +33,6 @@ NPairSkipSizeOff2onOneside::NPairSkipSizeOff2onOneside(LAMMPS *lmp) : iskip and ijskip flag which atom types and type pairs to skip parent non-skip list used newton off and was not onesided, this skip list is newton on and onesided - if list requests it, preserve shear history via fix shear/history ------------------------------------------------------------------------- */ void NPairSkipSizeOff2onOneside::build(NeighList *list) @@ -44,15 +41,6 @@ void NPairSkipSizeOff2onOneside::build(NeighList *list) tagint jtag; int *surf,*jlist; - int *npartner; - tagint **partner; - double **shearpartner; - int **firsttouch; - double **firstshear; - MyPage *ipage_touch; - MyPage *dpage_shear; - NeighList *listhistory; - tagint *tag = atom->tag; int *type = atom->type; int nlocal = atom->nlocal; @@ -73,28 +61,8 @@ void NPairSkipSizeOff2onOneside::build(NeighList *list) if (domain->dimension == 2) surf = atom->line; else surf = atom->tri; - FixShearHistory *fix_history = (FixShearHistory *) list->fix_history; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nlocal + atom->nghost; - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - listhistory = list->listhistory; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage; - dpage_shear = listhistory->dpage; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - } - int inum = 0; ipage->reset(); - if (fix_history) { - ipage_touch->reset(); - dpage_shear->reset(); - } // two loops over parent list required, one to count, one to store // because onesided constraint means pair I,J may be stored with I or J @@ -139,7 +107,7 @@ void NPairSkipSizeOff2onOneside::build(NeighList *list) } } - // allocate all per-atom neigh list chunks, including history + // allocate all per-atom neigh list chunks for (i = 0; i < nlocal; i++) { if (numneigh[i] == 0) continue; @@ -147,10 +115,6 @@ void NPairSkipSizeOff2onOneside::build(NeighList *list) firstneigh[i] = ipage->get(n); if (ipage->status()) error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); - if (fix_history) { - firsttouch[i] = ipage_touch->get(n); - firstshear[i] = dpage_shear->get(dnum*n); - } } // second loop over atoms in other list to store neighbors @@ -189,32 +153,11 @@ void NPairSkipSizeOff2onOneside::build(NeighList *list) // OK, b/c there is no special list flagging for surfs firstneigh[i][numneigh[i]] = j; - - // no numeric test for current touch - // just use FSH partner list to infer it - // would require complex calculation for surfs - - if (fix_history) { - jtag = tag[j]; - n = numneigh[i]; - nn = dnum*n; - for (m = 0; m < npartner[i]; m++) - if (partner[i][m] == jtag) break; - if (m < npartner[i]) { - firsttouch[i][n] = 1; - memcpy(&firstshear[i][nn],&shearpartner[i][dnum*m],dnumbytes); - } else { - firsttouch[i][n] = 0; - memcpy(&firstshear[i][nn],zeroes,dnumbytes); - } - } - numneigh[i]++; if (flip) i = j; } // only add atom I to ilist if it has neighbors - // fix shear/history allows for this in pre_exchange_onesided() if (numneigh[i]) ilist[inum++] = i; } diff --git a/src/verlet.cpp b/src/verlet.cpp index b242b00722..d74906556b 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -90,10 +90,9 @@ void Verlet::setup(int flag) if (comm->me == 0 && screen) { fprintf(screen,"Setting up Verlet run ...\n"); if (flag) { - fprintf(screen," Unit style : %s\n", update->unit_style); - fprintf(screen," Current step : " BIGINT_FORMAT "\n", - update->ntimestep); - fprintf(screen," Time step : %g\n", update->dt); + fprintf(screen," Unit style : %s\n",update->unit_style); + fprintf(screen," Current step : " BIGINT_FORMAT "\n",update->ntimestep); + fprintf(screen," Time step : %g\n",update->dt); timer->print_timeout(screen); } } @@ -122,6 +121,7 @@ void Verlet::setup(int flag) domain->box_too_small_check(); modify->setup_pre_neighbor(); neighbor->build(); + modify->setup_post_neighbor(); neighbor->ncalls = 0; // compute all forces @@ -183,6 +183,7 @@ void Verlet::setup_minimal(int flag) domain->box_too_small_check(); modify->setup_pre_neighbor(); neighbor->build(); + modify->setup_post_neighbor(); neighbor->ncalls = 0; } @@ -227,6 +228,7 @@ void Verlet::run(int n) int n_post_integrate = modify->n_post_integrate; int n_pre_exchange = modify->n_pre_exchange; int n_pre_neighbor = modify->n_pre_neighbor; + int n_post_neighbor = modify->n_post_neighbor; int n_pre_force = modify->n_pre_force; int n_pre_reverse = modify->n_pre_reverse; int n_post_force = modify->n_post_force; @@ -284,6 +286,10 @@ void Verlet::run(int n) } neighbor->build(); timer->stamp(Timer::NEIGH); + if (n_post_neighbor) { + modify->post_neighbor(); + timer->stamp(Timer::MODIFY); + } } // force computations