diff --git a/src/CORESHELL/pair_lj_cut_coul_long_cs.cpp b/src/CORESHELL/pair_lj_cut_coul_long_cs.cpp index e2ffda148f..d418cf20af 100644 --- a/src/CORESHELL/pair_lj_cut_coul_long_cs.cpp +++ b/src/CORESHELL/pair_lj_cut_coul_long_cs.cpp @@ -225,10 +225,10 @@ void PairLJCutCoulLongCS::compute_inner() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -311,10 +311,10 @@ void PairLJCutCoulLongCS::compute_middle() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; @@ -412,10 +412,10 @@ void PairLJCutCoulLongCS::compute_outer(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; diff --git a/src/DIPOLE/pair_lj_long_dipole_long.cpp b/src/DIPOLE/pair_lj_long_dipole_long.cpp index b833b250d4..c9b2b3f4af 100644 --- a/src/DIPOLE/pair_lj_long_dipole_long.cpp +++ b/src/DIPOLE/pair_lj_long_dipole_long.cpp @@ -263,22 +263,6 @@ void PairLJLongDipoleLong::init_style() if (force->kspace) g_ewald = force->kspace->g_ewald; } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJLongDipoleLong::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; - - if (id) - error->all(FLERR,"Pair style lj/long/dipole/long does not currently support respa"); -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/DIPOLE/pair_lj_long_dipole_long.h b/src/DIPOLE/pair_lj_long_dipole_long.h index f9fa10af11..2ace9ca301 100644 --- a/src/DIPOLE/pair_lj_long_dipole_long.h +++ b/src/DIPOLE/pair_lj_long_dipole_long.h @@ -34,7 +34,6 @@ class PairLJLongDipoleLong : 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 *); 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/KOKKOS/neigh_list_kokkos.cpp b/src/KOKKOS/neigh_list_kokkos.cpp index caf2dfee56..04454e53cb 100644 --- a/src/KOKKOS/neigh_list_kokkos.cpp +++ b/src/KOKKOS/neigh_list_kokkos.cpp @@ -49,15 +49,6 @@ void NeighListKokkos::grow(int nmax) d_neighbors = typename ArrayTypes::t_neighbors_2d("neighlist:neighbors", maxatoms,maxneighs); - - memory->sfree(firstneigh); - memory->sfree(firstdouble); - - firstneigh = (int **) memory->smalloc(maxatoms*sizeof(int *), - "neighlist:firstneigh"); - if (dnum) - firstdouble = (double **) memory->smalloc(maxatoms*sizeof(double *), - "neighlist:firstdouble"); } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/npair_copy_kokkos.cpp b/src/KOKKOS/npair_copy_kokkos.cpp index 6835d8c1b5..8702816033 100644 --- a/src/KOKKOS/npair_copy_kokkos.cpp +++ b/src/KOKKOS/npair_copy_kokkos.cpp @@ -41,10 +41,7 @@ void NPairCopyKokkos::build(NeighList *list) list->gnum = listcopy->gnum; list->ilist = listcopy->ilist; list->numneigh = listcopy->numneigh; - list->firstneigh = listcopy->firstneigh; - list->firstdouble = listcopy->firstdouble; list->ipage = listcopy->ipage; - list->dpage = listcopy->dpage; NeighListKokkos* list_kk = (NeighListKokkos*) list; NeighListKokkos* listcopy_kk = (NeighListKokkos*) list->listcopy; diff --git a/src/KSPACE/pair_buck_long_coul_long.cpp b/src/KSPACE/pair_buck_long_coul_long.cpp index 4cfb9b7267..7df8ebac68 100644 --- a/src/KSPACE/pair_buck_long_coul_long.cpp +++ b/src/KSPACE/pair_buck_long_coul_long.cpp @@ -233,7 +233,8 @@ void PairBuckLongCoulLong::init_style() if (!atom->q_flag && (ewald_order&(1<<1))) error->all(FLERR, - "Invoking coulombic in pair style buck/long/coul/long requires atom attribute q"); + "Invoking coulombic in pair style buck/long/coul/long " + "requires atom attribute q"); // ensure use of KSpace long-range solver, set two g_ewalds @@ -258,51 +259,25 @@ void PairBuckLongCoulLong::init_style() if (force->kspace->neighrequest_flag) { int irequest; + int respa = 0; if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { - int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + } - if (respa == 0) irequest = neighbor->request(this,instance_me); - else if (respa == 1) { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } else { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 2; - neighbor->requests[irequest]->respamiddle = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } + irequest = neighbor->request(this,instance_me); - } else irequest = neighbor->request(this,instance_me); + if (respa >= 1) { + neighbor->requests[irequest]->respaouter = 1; + neighbor->requests[irequest]->respainner = 1; + } + if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; } cut_coulsq = cut_coul * cut_coul; } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairBuckLongCoulLong::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ @@ -651,14 +626,14 @@ void PairBuckLongCoulLong::compute_inner() double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi; vector xi, d; - ineighn = (ineigh = listinner->ilist) + listinner->inum; + ineighn = (ineigh = list->ilist_inner) + list->inum_inner; for (; ineighfirstneigh[i])+listinner->numneigh[i]; + jneighn = (jneigh = list->firstneigh_inner[i])+list->numneigh_inner[i]; for (; jneighilist)+listmiddle->inum; + ineighn = (ineigh = list->ilist_middle)+list->inum_middle; for (; ineighfirstneigh[i])+listmiddle->numneigh[i]; + jneighn = (jneigh = list->firstneigh_middle[i])+list->numneigh_middle[i]; for (; jneighilist)+listouter->inum; + ineighn = (ineigh = list->ilist)+list->inum; for (; ineighfirstneigh[i])+listouter->numneigh[i]; + jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i]; for (; jneighnewton_pair; double qqrd2e = force->qqrd2e; - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; // loop over neighbors of my atoms @@ -320,10 +320,10 @@ void PairLJCharmmCoulLong::compute_middle() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; // loop over neighbors of my atoms @@ -417,10 +417,10 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; // loop over neighbors of my atoms @@ -687,36 +687,23 @@ void PairLJCharmmCoulLong::init_style() error->all(FLERR, "Pair style lj/charmm/coul/long requires atom attribute q"); - // request regular or rRESPA neighbor lists + // request regular or rRESPA neighbor list int irequest; + int respa = 0; if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { - int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + } - if (respa == 0) irequest = neighbor->request(this,instance_me); - else if (respa == 1) { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } else { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 2; - neighbor->requests[irequest]->respamiddle = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } + irequest = neighbor->request(this,instance_me); - } else irequest = neighbor->request(this,instance_me); + if (respa >= 1) { + neighbor->requests[irequest]->respaouter = 1; + neighbor->requests[irequest]->respainner = 1; + } + if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; // require cut_lj_inner < cut_lj @@ -767,19 +754,6 @@ void PairLJCharmmCoulLong::init_style() if (ncoultablebits) init_tables(cut_coul,cut_respa); } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJCharmmCoulLong::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_lj_charmm_coul_long.h b/src/KSPACE/pair_lj_charmm_coul_long.h index 1544f3bc14..95c6d0d1c7 100644 --- a/src/KSPACE/pair_lj_charmm_coul_long.h +++ b/src/KSPACE/pair_lj_charmm_coul_long.h @@ -33,7 +33,6 @@ class PairLJCharmmCoulLong : public Pair { virtual void settings(int, char **); void coeff(int, char **); virtual void init_style(); - void init_list(int, class NeighList *); virtual double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); diff --git a/src/KSPACE/pair_lj_charmm_coul_msm.cpp b/src/KSPACE/pair_lj_charmm_coul_msm.cpp index 76c9ef0cc7..00617c0bf2 100644 --- a/src/KSPACE/pair_lj_charmm_coul_msm.cpp +++ b/src/KSPACE/pair_lj_charmm_coul_msm.cpp @@ -278,10 +278,10 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; diff --git a/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp b/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp index 30d8ab64b6..859f421763 100644 --- a/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp +++ b/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp @@ -274,10 +274,10 @@ void PairLJCharmmfswCoulLong::compute_inner() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -359,10 +359,10 @@ void PairLJCharmmfswCoulLong::compute_middle() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; @@ -465,10 +465,10 @@ void PairLJCharmmfswCoulLong::compute_outer(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; @@ -824,19 +824,6 @@ void PairLJCharmmfswCoulLong::init_style() if (ncoultablebits) init_tables(cut_coul,cut_respa); } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJCharmmfswCoulLong::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_lj_charmmfsw_coul_long.h b/src/KSPACE/pair_lj_charmmfsw_coul_long.h index 650a908e48..135b82ea72 100644 --- a/src/KSPACE/pair_lj_charmmfsw_coul_long.h +++ b/src/KSPACE/pair_lj_charmmfsw_coul_long.h @@ -33,7 +33,6 @@ class PairLJCharmmfswCoulLong : public Pair { virtual void settings(int, char **); void coeff(int, char **); virtual void init_style(); - void init_list(int, class NeighList *); virtual double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); diff --git a/src/KSPACE/pair_lj_cut_coul_long.cpp b/src/KSPACE/pair_lj_cut_coul_long.cpp index f8be9fdb79..3096df2b01 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long.cpp @@ -224,10 +224,10 @@ void PairLJCutCoulLong::compute_inner() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -309,10 +309,10 @@ void PairLJCutCoulLong::compute_middle() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; @@ -410,10 +410,10 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; @@ -656,36 +656,23 @@ void PairLJCutCoulLong::init_style() if (!atom->q_flag) error->all(FLERR,"Pair style lj/cut/coul/long requires atom attribute q"); - // request regular or rRESPA neighbor lists + // request regular or rRESPA neighbor list int irequest; + int respa = 0; if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { - int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + } - if (respa == 0) irequest = neighbor->request(this,instance_me); - else if (respa == 1) { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } else { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 2; - neighbor->requests[irequest]->respamiddle = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } + irequest = neighbor->request(this,instance_me); - } else irequest = neighbor->request(this,instance_me); + if (respa >= 1) { + neighbor->requests[irequest]->respaouter = 1; + neighbor->requests[irequest]->respainner = 1; + } + if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; cut_coulsq = cut_coul * cut_coul; @@ -707,19 +694,6 @@ void PairLJCutCoulLong::init_style() if (ncoultablebits) init_tables(cut_coul,cut_respa); } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_lj_cut_coul_long.h b/src/KSPACE/pair_lj_cut_coul_long.h index 886542d075..e6f97c088d 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.h +++ b/src/KSPACE/pair_lj_cut_coul_long.h @@ -33,7 +33,6 @@ class PairLJCutCoulLong : public Pair { virtual void settings(int, char **); void coeff(int, char **); virtual void init_style(); - void init_list(int, class NeighList *); virtual double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); diff --git a/src/KSPACE/pair_lj_cut_coul_msm.cpp b/src/KSPACE/pair_lj_cut_coul_msm.cpp index e3b3f58fcb..9f901db9fc 100644 --- a/src/KSPACE/pair_lj_cut_coul_msm.cpp +++ b/src/KSPACE/pair_lj_cut_coul_msm.cpp @@ -265,10 +265,10 @@ void PairLJCutCoulMSM::compute_outer(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; diff --git a/src/KSPACE/pair_lj_long_coul_long.cpp b/src/KSPACE/pair_lj_long_coul_long.cpp index 7c6adfcb41..61b69011f1 100644 --- a/src/KSPACE/pair_lj_long_coul_long.cpp +++ b/src/KSPACE/pair_lj_long_coul_long.cpp @@ -253,51 +253,25 @@ void PairLJLongCoulLong::init_style() if (force->kspace->neighrequest_flag) { int irequest; - + int respa = 0; + if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { - int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; - - if (respa == 0) irequest = neighbor->request(this,instance_me); - else if (respa == 1) { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } else { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 2; - neighbor->requests[irequest]->respamiddle = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } - - } else irequest = neighbor->request(this,instance_me); + } + + irequest = neighbor->request(this,instance_me); + + if (respa >= 1) { + neighbor->requests[irequest]->respaouter = 1; + neighbor->requests[irequest]->respainner = 1; + } + if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; } cut_coulsq = cut_coul * cut_coul; } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJLongCoulLong::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ @@ -649,13 +623,13 @@ void PairLJLongCoulLong::compute_inner() double qri, *cut_ljsqi, *lj1i, *lj2i; vector xi, d; - ineighn = (ineigh = listinner->ilist)+listinner->inum; + ineighn = (ineigh = list->ilist_inner)+list->inum_inner; for (; ineighfirstneigh[i])+listinner->numneigh[i]; + jneighn = (jneigh = list->firstneigh_inner[i])+list->numneigh_inner[i]; for (; jneighilist)+listmiddle->inum; + ineighn = (ineigh = list->ilist_middle)+list->inum_middle; for (; ineighfirstneigh[i])+listmiddle->numneigh[i]; + jneighn = (jneigh = list->firstneigh_middle[i])+list->numneigh_middle[i]; for (; jneighilist)+listouter->inum; + ineighn = (ineigh = list->ilist)+list->inum; for (; ineighfirstneigh[i])+listouter->numneigh[i]; + jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i]; for (; jneighinum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; // loop over neighbors of my atoms @@ -769,10 +769,10 @@ void PairLJLongTIP4PLong::compute_middle() int ni; double *lj1i, *lj2i; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; // loop over neighbors of my atoms @@ -1049,10 +1049,10 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag) double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; // loop over neighbors of my atoms diff --git a/src/OPT/pair_lj_long_coul_long_opt.cpp b/src/OPT/pair_lj_long_coul_long_opt.cpp index 9004e5c93c..678d2d8bc4 100644 --- a/src/OPT/pair_lj_long_coul_long_opt.cpp +++ b/src/OPT/pair_lj_long_coul_long_opt.cpp @@ -726,7 +726,7 @@ void PairLJLongCoulLongOpt::eval_outer() double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; - ineighn = (ineigh = listouter->ilist)+listouter->inum; + ineighn = (ineigh = list->ilist)+list->inum; for (; ineighfirstneigh[i])+listouter->numneigh[i]; + jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i]; for (; jneighnewton_pair; double qqrd2e = force->qqrd2e; - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -315,10 +315,10 @@ void PairLJCharmmCoulLongSoft::compute_middle() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; @@ -428,10 +428,10 @@ void PairLJCharmmCoulLongSoft::compute_outer(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; @@ -758,19 +758,6 @@ void PairLJCharmmCoulLongSoft::init_style() g_ewald = force->kspace->g_ewald; } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJCharmmCoulLongSoft::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/USER-FEP/pair_lj_charmm_coul_long_soft.h b/src/USER-FEP/pair_lj_charmm_coul_long_soft.h index 7e52ec54b5..252c9f66f5 100644 --- a/src/USER-FEP/pair_lj_charmm_coul_long_soft.h +++ b/src/USER-FEP/pair_lj_charmm_coul_long_soft.h @@ -33,7 +33,6 @@ class PairLJCharmmCoulLongSoft : public Pair { 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 *); diff --git a/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp b/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp index f7c4084fe2..7be2ebabea 100644 --- a/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp +++ b/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp @@ -209,10 +209,10 @@ void PairLJCutCoulLongSoft::compute_inner() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -299,10 +299,10 @@ void PairLJCutCoulLongSoft::compute_middle() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; @@ -403,10 +403,10 @@ void PairLJCutCoulLongSoft::compute_outer(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; @@ -686,19 +686,6 @@ void PairLJCutCoulLongSoft::init_style() g_ewald = force->kspace->g_ewald; } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJCutCoulLongSoft::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/USER-FEP/pair_lj_cut_coul_long_soft.h b/src/USER-FEP/pair_lj_cut_coul_long_soft.h index a03be3814a..d49d1c8641 100644 --- a/src/USER-FEP/pair_lj_cut_coul_long_soft.h +++ b/src/USER-FEP/pair_lj_cut_coul_long_soft.h @@ -32,7 +32,6 @@ class PairLJCutCoulLongSoft : public Pair { virtual void settings(int, char **); void coeff(int, char **); virtual void init_style(); - void init_list(int, class NeighList *); virtual double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); diff --git a/src/USER-FEP/pair_lj_cut_soft.cpp b/src/USER-FEP/pair_lj_cut_soft.cpp index 8b6280a61a..9ae108fa33 100644 --- a/src/USER-FEP/pair_lj_cut_soft.cpp +++ b/src/USER-FEP/pair_lj_cut_soft.cpp @@ -164,10 +164,10 @@ void PairLJCutSoft::compute_inner() double *special_lj = force->special_lj; int newton_pair = force->newton_pair; - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -242,10 +242,10 @@ void PairLJCutSoft::compute_middle() double *special_lj = force->special_lj; int newton_pair = force->newton_pair; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; @@ -333,10 +333,10 @@ void PairLJCutSoft::compute_outer(int eflag, int vflag) double *special_lj = force->special_lj; int newton_pair = force->newton_pair; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; @@ -556,19 +556,6 @@ void PairLJCutSoft::init_style() } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJCutSoft::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/USER-FEP/pair_lj_cut_soft.h b/src/USER-FEP/pair_lj_cut_soft.h index 50ce685e5c..46202d78a8 100644 --- a/src/USER-FEP/pair_lj_cut_soft.h +++ b/src/USER-FEP/pair_lj_cut_soft.h @@ -32,7 +32,6 @@ class PairLJCutSoft : public Pair { virtual void settings(int, char **); void coeff(int, char **); virtual void init_style(); - void init_list(int, class NeighList *); virtual double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); diff --git a/src/USER-OMP/fix_shear_history_omp.cpp b/src/USER-OMP/fix_neigh_history_omp.cpp similarity index 61% rename from src/USER-OMP/fix_shear_history_omp.cpp rename to src/USER-OMP/fix_neigh_history_omp.cpp index 4180e0af41..05b8fa41b9 100644 --- a/src/USER-OMP/fix_shear_history_omp.cpp +++ b/src/USER-OMP/fix_neigh_history_omp.cpp @@ -13,7 +13,7 @@ #include #include -#include "fix_shear_history_omp.h" +#include "fix_neigh_history_omp.h" #include "atom.h" #include "comm.h" #include "neighbor.h" @@ -31,15 +31,38 @@ using namespace LAMMPS_NS; using namespace FixConst; + +FixNeighHistoryOMP::FixNeighHistoryOMP(class LAMMPS *lmp,int narg,char **argv) + : FixNeighHistory(lmp,narg,argv) { + + if (onesided) + error->all(FLERR,"tri/lj and line/lj are not supported by USER-OMP"); + if (!newton_pair) + error->all(FLERR,"Newton off for granular is not supported by USER-OMP"); +} + + /* ---------------------------------------------------------------------- - copy shear partner info from neighbor lists to per-atom arrays - so it can be exchanged with those atoms + 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 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 + + the USER-OMP version only supports newton on ------------------------------------------------------------------------- */ -void FixShearHistoryOMP::pre_exchange() +void FixNeighHistoryOMP::pre_exchange() { const int nthreads = comm->nthreads; - maxtouch = 0; + maxpartner = 0; #if defined(_OPENMP) #pragma omp parallel default(none) @@ -54,16 +77,16 @@ void FixShearHistoryOMP::pre_exchange() int i,j,ii,jj,m,n,inum,jnum; int *ilist,*jlist,*numneigh,**firstneigh; - int *touch,**firsttouch; - double *shear,*shearj,*allshear,**firstshear; + int *allflags,**firstflag; + double *allvalues,*onevalues,*jvalues; - MyPage &ipg = ipage[tid]; - MyPage &dpg = dpage[tid]; + MyPage &ipg = ipage_atom[tid]; + MyPage &dpg = dpage_atom[tid]; ipg.reset(); dpg.reset(); // 1st loop over neighbor list - // calculate nparter for each owned atom + // calculate npartner for each owned atom tagint *tag = atom->tag; @@ -72,8 +95,6 @@ void FixShearHistoryOMP::pre_exchange() ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - firsttouch = list->listhistory->firstneigh; - firstshear = list->listhistory->firstdouble; int nlocal_neigh = 0; if (inum) nlocal_neigh = ilist[inum-1] + 1; @@ -94,10 +115,10 @@ void FixShearHistoryOMP::pre_exchange() 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]) { if ((i >= lfrom) && (i < lto)) npartner[i]++; @@ -116,55 +137,55 @@ void FixShearHistoryOMP::pre_exchange() if ((i >= lfrom) && (i < lto)) { n = npartner[i]; partner[i] = ipg.get(n); - shearpartner[i] = dpg.get(dnum*n); - if (partner[i] == NULL || shearpartner[i] == NULL) - error->one(FLERR,"Shear history overflow, boost neigh_modify one"); + valuepartner[i] = dpg.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 my atoms - // re-zero npartner to use as counter for all my atoms + // store partner IDs and values for owned+ghost atoms + // re-zero npartner to use as counter for (i = lfrom; i < lto; i++) npartner[i] = 0; 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[3*jj]; + if (allflags[jj]) { + onevalues = &allvalues[3*jj]; j = jlist[jj]; j &= NEIGHMASK; if ((i >= lfrom) && (i < lto)) { m = npartner[i]++; partner[i][m] = tag[j]; - memcpy(&shearpartner[i][dnum*m],shear,dnumbytes); + memcpy(&valuepartner[i][dnum*m],onevalues,dnumbytes); } if ((j >= lfrom) && (j < lto)) { 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 - maxtouch = m = 0; + // set maxpartner = max # of partners of any owned atom + maxpartner = m = 0; for (i = lfrom; i < lto; i++) m = MAX(m,npartner[i]); #if defined(_OPENMP) #pragma omp critical #endif - maxtouch = MAX(m,maxtouch); + maxpartner = MAX(m,maxpartner); } } diff --git a/src/USER-OMP/fix_shear_history_omp.h b/src/USER-OMP/fix_neigh_history_omp.h similarity index 70% rename from src/USER-OMP/fix_shear_history_omp.h rename to src/USER-OMP/fix_neigh_history_omp.h index 95281b2afc..83472391a0 100644 --- a/src/USER-OMP/fix_shear_history_omp.h +++ b/src/USER-OMP/fix_neigh_history_omp.h @@ -13,22 +13,21 @@ #ifdef FIX_CLASS -FixStyle(SHEAR_HISTORY/omp,FixShearHistoryOMP) +FixStyle(NEIGH_HISTORY/omp,FixNeighHistoryOMP) #else -#ifndef LMP_FIX_SHEAR_HISTORY_OMP_H -#define LMP_FIX_SHEAR_HISTORY_OMP_H +#ifndef LMP_FIX_NEIGH_HISTORY_OMP_H +#define LMP_FIX_NEIGH_HISTORY_OMP_H -#include "fix_shear_history.h" +#include "fix_neigh_history.h" namespace LAMMPS_NS { -class FixShearHistoryOMP : public FixShearHistory { +class FixNeighHistoryOMP : public FixNeighHistory { public: - FixShearHistoryOMP(class LAMMPS *lmp, int narg, char **argv) - : FixShearHistory(lmp,narg,argv) {}; + FixNeighHistoryOMP(class LAMMPS *lmp, int narg, char **argv); virtual void pre_exchange(); }; diff --git a/src/USER-OMP/npair_half_respa_bin_newtoff_omp.cpp b/src/USER-OMP/npair_half_respa_bin_newtoff_omp.cpp index 45add87092..f094691b71 100644 --- a/src/USER-OMP/npair_half_respa_bin_newtoff_omp.cpp +++ b/src/USER-OMP/npair_half_respa_bin_newtoff_omp.cpp @@ -45,12 +45,10 @@ void NPairHalfRespaBinNewtoffOmp::build(NeighList *list) NPAIR_OMP_INIT; - NeighList *listinner = list->listinner; - NeighList *listmiddle = list->listmiddle; const int respamiddle = list->respamiddle; #if defined(_OPENMP) -#pragma omp parallel default(none) shared(list,listinner,listmiddle) +#pragma omp parallel default(none) shared(list) #endif NPAIR_OMP_SETUP(nlocal); @@ -77,26 +75,26 @@ void NPairHalfRespaBinNewtoffOmp::build(NeighList *list) int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; int *ilist_middle,*numneigh_middle,**firstneigh_middle; if (respamiddle) { - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; } // each thread has its own page allocator MyPage &ipage = list->ipage[tid]; - MyPage &ipage_inner = listinner->ipage[tid]; + MyPage &ipage_inner = list->ipage_inner[tid]; ipage.reset(); ipage_inner.reset(); MyPage *ipage_middle; if (respamiddle) { - ipage_middle = listmiddle->ipage + tid; + ipage_middle = list->ipage_middle + tid; ipage_middle->reset(); } @@ -199,6 +197,6 @@ void NPairHalfRespaBinNewtoffOmp::build(NeighList *list) } NPAIR_OMP_CLOSE; list->inum = nlocal; - listinner->inum = nlocal; - if (respamiddle) listmiddle->inum = nlocal; + list->inum_inner = nlocal; + if (respamiddle) list->inum_middle = nlocal; } diff --git a/src/USER-OMP/npair_half_respa_bin_newton_omp.cpp b/src/USER-OMP/npair_half_respa_bin_newton_omp.cpp index ee6b9b7501..de7ef5f7d5 100644 --- a/src/USER-OMP/npair_half_respa_bin_newton_omp.cpp +++ b/src/USER-OMP/npair_half_respa_bin_newton_omp.cpp @@ -44,12 +44,10 @@ void NPairHalfRespaBinNewtonOmp::build(NeighList *list) NPAIR_OMP_INIT; - NeighList *listinner = list->listinner; - NeighList *listmiddle = list->listmiddle; const int respamiddle = list->respamiddle; #if defined(_OPENMP) -#pragma omp parallel default(none) shared(list,listinner,listmiddle) +#pragma omp parallel default(none) shared(list) #endif NPAIR_OMP_SETUP(nlocal); @@ -76,26 +74,26 @@ void NPairHalfRespaBinNewtonOmp::build(NeighList *list) int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; int *ilist_middle,*numneigh_middle,**firstneigh_middle; if (respamiddle) { - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; } // each thread has its own page allocator MyPage &ipage = list->ipage[tid]; - MyPage &ipage_inner = listinner->ipage[tid]; + MyPage &ipage_inner = list->ipage_inner[tid]; ipage.reset(); ipage_inner.reset(); MyPage *ipage_middle; if (respamiddle) { - ipage_middle = listmiddle->ipage + tid; + ipage_middle = list->ipage_middle + tid; ipage_middle->reset(); } @@ -245,6 +243,6 @@ void NPairHalfRespaBinNewtonOmp::build(NeighList *list) } NPAIR_OMP_CLOSE; list->inum = nlocal; - listinner->inum = nlocal; - if (respamiddle) listmiddle->inum = nlocal; + list->inum_inner = nlocal; + if (respamiddle) list->inum_middle = nlocal; } diff --git a/src/USER-OMP/npair_half_respa_bin_newton_tri_omp.cpp b/src/USER-OMP/npair_half_respa_bin_newton_tri_omp.cpp index fbb512ba64..f20d101bc9 100644 --- a/src/USER-OMP/npair_half_respa_bin_newton_tri_omp.cpp +++ b/src/USER-OMP/npair_half_respa_bin_newton_tri_omp.cpp @@ -44,12 +44,10 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list) NPAIR_OMP_INIT; - NeighList *listinner = list->listinner; - NeighList *listmiddle = list->listmiddle; const int respamiddle = list->respamiddle; #if defined(_OPENMP) -#pragma omp parallel default(none) shared(list,listinner,listmiddle) +#pragma omp parallel default(none) shared(list) #endif NPAIR_OMP_SETUP(nlocal); @@ -76,26 +74,26 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list) int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; int *ilist_middle,*numneigh_middle,**firstneigh_middle; if (respamiddle) { - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; } // each thread has its own page allocator MyPage &ipage = list->ipage[tid]; - MyPage &ipage_inner = listinner->ipage[tid]; + MyPage &ipage_inner = list->ipage_inner[tid]; ipage.reset(); ipage_inner.reset(); MyPage *ipage_middle; if (respamiddle) { - ipage_middle = listmiddle->ipage + tid; + ipage_middle = list->ipage_middle + tid; ipage_middle->reset(); } @@ -206,6 +204,6 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list) } NPAIR_OMP_CLOSE; list->inum = nlocal; - listinner->inum = nlocal; - if (respamiddle) listmiddle->inum = nlocal; + list->inum_inner = nlocal; + if (respamiddle) list->inum_middle = nlocal; } diff --git a/src/USER-OMP/npair_half_respa_nsq_newtoff_omp.cpp b/src/USER-OMP/npair_half_respa_nsq_newtoff_omp.cpp index 5ee71bebad..0f726cdd7f 100644 --- a/src/USER-OMP/npair_half_respa_nsq_newtoff_omp.cpp +++ b/src/USER-OMP/npair_half_respa_nsq_newtoff_omp.cpp @@ -46,12 +46,10 @@ void NPairHalfRespaNsqNewtoffOmp::build(NeighList *list) NPAIR_OMP_INIT; - NeighList *listinner = list->listinner; - NeighList *listmiddle = list->listmiddle; const int respamiddle = list->respamiddle; #if defined(_OPENMP) -#pragma omp parallel default(none) shared(list,listinner,listmiddle) +#pragma omp parallel default(none) shared(list) #endif NPAIR_OMP_SETUP(nlocal); @@ -80,26 +78,26 @@ void NPairHalfRespaNsqNewtoffOmp::build(NeighList *list) int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; int *ilist_middle,*numneigh_middle,**firstneigh_middle; if (respamiddle) { - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; } // each thread has its own page allocator MyPage &ipage = list->ipage[tid]; - MyPage &ipage_inner = listinner->ipage[tid]; + MyPage &ipage_inner = list->ipage_inner[tid]; ipage.reset(); ipage_inner.reset(); MyPage *ipage_middle; if (respamiddle) { - ipage_middle = listmiddle->ipage + tid; + ipage_middle = list->ipage_middle + tid; ipage_middle->reset(); } @@ -193,6 +191,6 @@ void NPairHalfRespaNsqNewtoffOmp::build(NeighList *list) } NPAIR_OMP_CLOSE; list->inum = nlocal; - listinner->inum = nlocal; - if (respamiddle) listmiddle->inum = nlocal; + list->inum_inner = nlocal; + if (respamiddle) list->inum_middle = nlocal; } diff --git a/src/USER-OMP/npair_half_respa_nsq_newton_omp.cpp b/src/USER-OMP/npair_half_respa_nsq_newton_omp.cpp index 89cff732c9..2783e1255e 100644 --- a/src/USER-OMP/npair_half_respa_nsq_newton_omp.cpp +++ b/src/USER-OMP/npair_half_respa_nsq_newton_omp.cpp @@ -47,12 +47,10 @@ void NPairHalfRespaNsqNewtonOmp::build(NeighList *list) NPAIR_OMP_INIT; - NeighList *listinner = list->listinner; - NeighList *listmiddle = list->listmiddle; const int respamiddle = list->respamiddle; #if defined(_OPENMP) -#pragma omp parallel default(none) shared(list,listinner,listmiddle) +#pragma omp parallel default(none) shared(list) #endif NPAIR_OMP_SETUP(nlocal); @@ -81,26 +79,26 @@ void NPairHalfRespaNsqNewtonOmp::build(NeighList *list) int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; int *ilist_middle,*numneigh_middle,**firstneigh_middle; if (respamiddle) { - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; } // each thread has its own page allocator MyPage &ipage = list->ipage[tid]; - MyPage &ipage_inner = listinner->ipage[tid]; + MyPage &ipage_inner = list->ipage_inner[tid]; ipage.reset(); ipage_inner.reset(); MyPage *ipage_middle; if (respamiddle) { - ipage_middle = listmiddle->ipage + tid; + ipage_middle = list->ipage_middle + tid; ipage_middle->reset(); } @@ -212,6 +210,6 @@ void NPairHalfRespaNsqNewtonOmp::build(NeighList *list) } NPAIR_OMP_CLOSE; list->inum = nlocal; - listinner->inum = nlocal; - if (respamiddle) listmiddle->inum = nlocal; + list->inum_inner = nlocal; + if (respamiddle) list->inum_middle = nlocal; } diff --git a/src/USER-OMP/npair_half_size_bin_newtoff_omp.cpp b/src/USER-OMP/npair_half_size_bin_newtoff_omp.cpp index 120658b714..6a1cb46ea6 100644 --- a/src/USER-OMP/npair_half_size_bin_newtoff_omp.cpp +++ b/src/USER-OMP/npair_half_size_bin_newtoff_omp.cpp @@ -18,9 +18,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" @@ -34,7 +31,6 @@ NPairHalfSizeBinNewtoffOmp::NPairHalfSizeBinNewtoffOmp(LAMMPS *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) @@ -43,30 +39,20 @@ NPairHalfSizeBinNewtoffOmp::NPairHalfSizeBinNewtoffOmp(LAMMPS *lmp) : void NPairHalfSizeBinNewtoffOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; - - FixShearHistory * const fix_history = (FixShearHistory *) list->fix_history; - NeighList * listhistory = list->listhistory; + const int history = list->history; + const int mask_history = 3 << SBBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) -#pragma omp parallel default(none) shared(list,listhistory) +#pragma omp parallel default(none) shared(list) #endif NPAIR_OMP_SETUP(nlocal); - 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; - MyPage *ipage_touch; - MyPage *dpage_shear; - - int *npartner; - tagint **partner; - double **shearpartner; - int **firsttouch; - double **firstshear; + int *neighptr; // loop over each atom, storing neighbors @@ -85,29 +71,10 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list) MyPage &ipage = list->ipage[tid]; ipage.reset(); - if (fix_history) { - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage+tid; - dpage_shear = listhistory->dpage+tid; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - ipage_touch->reset(); - dpage_shear->reset(); - } - for (i = ifrom; i < ito; 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]; @@ -133,29 +100,10 @@ void NPairHalfSizeBinNewtoffOmp::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; } } } @@ -166,13 +114,6 @@ void NPairHalfSizeBinNewtoffOmp::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); - } } NPAIR_OMP_CLOSE; list->inum = nlocal; diff --git a/src/USER-OMP/npair_half_size_bin_newton_omp.cpp b/src/USER-OMP/npair_half_size_bin_newton_omp.cpp index cf0c6d20fe..d8e1e6da44 100644 --- a/src/USER-OMP/npair_half_size_bin_newton_omp.cpp +++ b/src/USER-OMP/npair_half_size_bin_newton_omp.cpp @@ -18,9 +18,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" @@ -34,7 +31,6 @@ NPairHalfSizeBinNewtonOmp::NPairHalfSizeBinNewtonOmp(LAMMPS *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 ------------------------------------------------------------------------- */ @@ -42,36 +38,20 @@ NPairHalfSizeBinNewtonOmp::NPairHalfSizeBinNewtonOmp(LAMMPS *lmp) : void NPairHalfSizeBinNewtonOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; - - FixShearHistory * const fix_history = (FixShearHistory *) list->fix_history; - NeighList * listhistory = list->listhistory; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nlocal + atom->nghost; - } + const int history = list->history; + const int mask_history = 3 << SBBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) -#pragma omp parallel default(none) shared(list,listhistory) +#pragma omp parallel default(none) shared(list) #endif NPAIR_OMP_SETUP(nlocal); - 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; - MyPage *ipage_touch; - MyPage *dpage_shear; - - int *npartner; - tagint **partner; - double **shearpartner; - int **firsttouch; - double **firstshear; - - // loop over each atom, storing neighbors + int *neighptr; double **x = atom->x; double *radius = atom->radius; @@ -88,29 +68,10 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list) MyPage &ipage = list->ipage[tid]; ipage.reset(); - if (fix_history) { - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage+tid; - dpage_shear = listhistory->dpage+tid; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - ipage_touch->reset(); - dpage_shear->reset(); - } - for (i = ifrom; i < ito; 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]; @@ -140,29 +101,10 @@ void NPairHalfSizeBinNewtonOmp::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; } } @@ -181,29 +123,10 @@ void NPairHalfSizeBinNewtonOmp::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; } } } @@ -214,13 +137,6 @@ void NPairHalfSizeBinNewtonOmp::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); - } } NPAIR_OMP_CLOSE; list->inum = nlocal; diff --git a/src/USER-OMP/npair_half_size_bin_newton_tri_omp.cpp b/src/USER-OMP/npair_half_size_bin_newton_tri_omp.cpp index da04eebd1e..b02bfa345e 100644 --- a/src/USER-OMP/npair_half_size_bin_newton_tri_omp.cpp +++ b/src/USER-OMP/npair_half_size_bin_newton_tri_omp.cpp @@ -17,8 +17,6 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "domain.h" #include "my_page.h" #include "error.h" @@ -32,7 +30,6 @@ NPairHalfSizeBinNewtonTriOmp::NPairHalfSizeBinNewtonTriOmp(LAMMPS *lmp) : /* ---------------------------------------------------------------------- size particles binned neighbor list construction with Newton's 3rd law for triclinic - no shear history is allowed for this option each owned atom i checks its own bin and other bins in triclinic stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ @@ -40,6 +37,8 @@ NPairHalfSizeBinNewtonTriOmp::NPairHalfSizeBinNewtonTriOmp(LAMMPS *lmp) : void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int history = list->history; + const int mask_history = 3 << SBBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -105,7 +104,12 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) radsum = radi + radius[j]; cutsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutsq) neighptr[n++] = j; + if (rsq <= cutsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } } } diff --git a/src/USER-OMP/npair_half_size_nsq_newtoff_omp.cpp b/src/USER-OMP/npair_half_size_nsq_newtoff_omp.cpp index f898ec3828..3c7b6b118f 100644 --- a/src/USER-OMP/npair_half_size_nsq_newtoff_omp.cpp +++ b/src/USER-OMP/npair_half_size_nsq_newtoff_omp.cpp @@ -19,9 +19,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" @@ -44,34 +41,20 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; - - FixShearHistory * const fix_history = (FixShearHistory *) list->fix_history; - NeighList * listhistory = list->listhistory; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nlocal + atom->nghost; - } + const int history = list->history; + const int mask_history = 3 << SBBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) -#pragma omp parallel default(none) shared(list,listhistory) +#pragma omp parallel default(none) shared(list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,m,n,nn,dnum,dnumbytes; + int i,j,m,n,nn; 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; + int *neighptr; double **x = atom->x; double *radius = atom->radius; @@ -89,29 +72,10 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) MyPage &ipage = list->ipage[tid]; ipage.reset(); - if (fix_history) { - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage+tid; - dpage_shear = listhistory->dpage+tid; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - ipage_touch->reset(); - dpage_shear->reset(); - } - for (i = ifrom; i < ito; 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]; @@ -132,29 +96,10 @@ void NPairHalfSizeNsqNewtoffOmp::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; } } @@ -164,13 +109,6 @@ void NPairHalfSizeNsqNewtoffOmp::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); - } } NPAIR_OMP_CLOSE; list->inum = nlocal; diff --git a/src/USER-OMP/npair_half_size_nsq_newton_omp.cpp b/src/USER-OMP/npair_half_size_nsq_newton_omp.cpp index a7caac372a..37a4181af7 100644 --- a/src/USER-OMP/npair_half_size_nsq_newton_omp.cpp +++ b/src/USER-OMP/npair_half_size_nsq_newton_omp.cpp @@ -19,9 +19,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" @@ -45,34 +42,20 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;; - - FixShearHistory * const fix_history = (FixShearHistory *) list->fix_history; - NeighList * listhistory = list->listhistory; - if (fix_history) { - fix_history->nlocal_neigh = nlocal; - fix_history->nall_neigh = nlocal+atom->nghost; - } + const int history = list->history; + const int mask_history = 3 << SBBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) -#pragma omp parallel default(none) shared(list,listhistory) +#pragma omp parallel default(none) shared(list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,m,n,nn,itag,jtag,dnum,dnumbytes; + int i,j,m,n,nn,itag,jtag; 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; + int *neighptr; double **x = atom->x; double *radius = atom->radius; @@ -90,29 +73,10 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) MyPage &ipage = list->ipage[tid]; ipage.reset(); - if (fix_history) { - npartner = fix_history->npartner; - partner = fix_history->partner; - shearpartner = fix_history->shearpartner; - firsttouch = listhistory->firstneigh; - firstshear = listhistory->firstdouble; - ipage_touch = listhistory->ipage+tid; - dpage_shear = listhistory->dpage+tid; - dnum = listhistory->dnum; - dnumbytes = dnum * sizeof(double); - ipage_touch->reset(); - dpage_shear->reset(); - } - for (i = ifrom; i < ito; 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]; @@ -150,29 +114,10 @@ void NPairHalfSizeNsqNewtonOmp::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; } } @@ -183,12 +128,6 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) 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); - } } NPAIR_OMP_CLOSE; list->inum = nlocal; diff --git a/src/USER-OMP/pair_buck_long_coul_long_omp.cpp b/src/USER-OMP/pair_buck_long_coul_long_omp.cpp index 87f9e2e321..f996372409 100644 --- a/src/USER-OMP/pair_buck_long_coul_long_omp.cpp +++ b/src/USER-OMP/pair_buck_long_coul_long_omp.cpp @@ -319,7 +319,7 @@ void PairBuckLongCoulLongOMP::compute_inner() const int nall = atom->nlocal + atom->nghost; const int nthreads = comm->nthreads; - const int inum = listinner->inum; + const int inum = list->inum_inner; #if defined(_OPENMP) #pragma omp parallel default(none) #endif @@ -343,7 +343,7 @@ void PairBuckLongCoulLongOMP::compute_middle() const int nall = atom->nlocal + atom->nghost; const int nthreads = comm->nthreads; - const int inum = listmiddle->inum; + const int inum = list->inum_middle; #if defined(_OPENMP) #pragma omp parallel default(none) @@ -373,7 +373,7 @@ void PairBuckLongCoulLongOMP::compute_outer(int eflag, int vflag) const int nall = atom->nlocal + atom->nghost; const int nthreads = comm->nthreads; - const int inum = listouter->inum; + const int inum = list->inum; #if defined(_OPENMP) #pragma omp parallel default(none) shared(eflag,vflag) @@ -811,7 +811,7 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t const double *x0 = x[0]; double *f0 = f[0], *fi = 0; - int *ilist = listinner->ilist; + int *ilist = list->ilist_inner; const int newton_pair = force->newton_pair; @@ -835,7 +835,7 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t memcpy(xi, x0+(i+(i<<1)), sizeof(vector)); cut_bucksqi = cut_bucksq[typei = type[i]]; buck1i = buck1[typei]; buck2i = buck2[typei]; rhoinvi = rhoinv[typei]; - jneighn = (jneigh = listinner->firstneigh[i])+listinner->numneigh[i]; + jneighn = (jneigh = list->firstneigh_inner[i])+list->numneigh_inner[i]; for (; jneighilist; + int *ilist = list->ilist_middle; const int newton_pair = force->newton_pair; @@ -932,7 +932,7 @@ void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const memcpy(xi, x0+(i+(i<<1)), sizeof(vector)); cut_bucksqi = cut_bucksq[typei = type[i]]; buck1i = buck1[typei]; buck2i = buck2[typei]; rhoinvi = rhoinv[typei]; - jneighn = (jneigh = listmiddle->firstneigh[i])+listmiddle->numneigh[i]; + jneighn = (jneigh = list->firstneigh_middle[i])+list->numneigh_middle[i]; for (; jneighilist; + int *ilist = list->ilist; int i, j, ii; int *jneigh, *jneighn, typei, typej, ni, respa_flag; @@ -1035,7 +1035,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t buckai = buck_a[typei]; buckci = buck_c[typei]; rhoinvi = rhoinv[typei]; cutsqi = cutsq[typei]; cut_bucksqi = cut_bucksq[typei]; memcpy(xi, x0+(i+(i<<1)), sizeof(vector)); - jneighn = (jneigh = listouter->firstneigh[i])+listouter->numneigh[i]; + jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i]; for (; jneigh #include "pair_gran_hertz_history_omp.h" +#include "fix_neigh_history.h" #include "atom.h" #include "comm.h" #include "fix.h" @@ -134,8 +135,8 @@ void PairGranHertzHistoryOMP::eval(int iifrom, int iito, ThrData * const thr) 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 diff --git a/src/USER-OMP/pair_gran_hooke_history_omp.cpp b/src/USER-OMP/pair_gran_hooke_history_omp.cpp index e507a63f7c..2e7d55aff0 100644 --- a/src/USER-OMP/pair_gran_hooke_history_omp.cpp +++ b/src/USER-OMP/pair_gran_hooke_history_omp.cpp @@ -14,6 +14,7 @@ #include #include "pair_gran_hooke_history_omp.h" +#include "fix_neigh_history.h" #include "atom.h" #include "comm.h" #include "fix.h" @@ -137,8 +138,8 @@ void PairGranHookeHistoryOMP::eval(int iifrom, int iito, ThrData * const thr) 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 diff --git a/src/USER-OMP/pair_lj_long_coul_long_omp.cpp b/src/USER-OMP/pair_lj_long_coul_long_omp.cpp index 28d4f229c8..c0c87e7481 100644 --- a/src/USER-OMP/pair_lj_long_coul_long_omp.cpp +++ b/src/USER-OMP/pair_lj_long_coul_long_omp.cpp @@ -317,7 +317,7 @@ void PairLJLongCoulLongOMP::compute_inner() const int nall = atom->nlocal + atom->nghost; const int nthreads = comm->nthreads; - const int inum = listinner->inum; + const int inum = list->inum_inner; #if defined(_OPENMP) #pragma omp parallel default(none) #endif @@ -341,7 +341,7 @@ void PairLJLongCoulLongOMP::compute_middle() const int nall = atom->nlocal + atom->nghost; const int nthreads = comm->nthreads; - const int inum = listmiddle->inum; + const int inum = list->inum_middle; #if defined(_OPENMP) #pragma omp parallel default(none) @@ -371,7 +371,7 @@ void PairLJLongCoulLongOMP::compute_outer(int eflag, int vflag) const int nall = atom->nlocal + atom->nghost; const int nthreads = comm->nthreads; - const int inum = listouter->inum; + const int inum = list->inum; #if defined(_OPENMP) #pragma omp parallel default(none) shared(eflag,vflag) @@ -805,7 +805,7 @@ void PairLJLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const thr const double *x0 = x[0]; double *f0 = f[0], *fi = 0; - int *ilist = listinner->ilist; + int *ilist = list->ilist_inner; const int newton_pair = force->newton_pair; @@ -828,7 +828,7 @@ void PairLJLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const thr memcpy(xi, x0+(i+(i<<1)), sizeof(vector)); cut_ljsqi = cut_ljsq[typei = type[i]]; lj1i = lj1[typei]; lj2i = lj2[typei]; - jneighn = (jneigh = listinner->firstneigh[i])+listinner->numneigh[i]; + jneighn = (jneigh = list->firstneigh_inner[i])+list->numneigh_inner[i]; for (; jneighilist; + int *ilist = list->ilist_middle; const int newton_pair = force->newton_pair; @@ -925,7 +925,7 @@ void PairLJLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const th memcpy(xi, x0+(i+(i<<1)), sizeof(vector)); cut_ljsqi = cut_ljsq[typei = type[i]]; lj1i = lj1[typei]; lj2i = lj2[typei]; - jneighn = (jneigh = listmiddle->firstneigh[i])+listmiddle->numneigh[i]; + jneighn = (jneigh = list->firstneigh_middle[i])+list->numneigh_middle[i]; for (; jneighilist; + int *ilist = list->ilist; int i, j, ii; int *jneigh, *jneighn, typei, typej, ni, respa_flag; @@ -1027,7 +1027,7 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr lj1i = lj1[typei]; lj2i = lj2[typei]; lj3i = lj3[typei]; lj4i = lj4[typei]; cutsqi = cutsq[typei]; cut_ljsqi = cut_ljsq[typei]; memcpy(xi, x0+(i+(i<<1)), sizeof(vector)); - jneighn = (jneigh = listouter->firstneigh[i])+listouter->numneigh[i]; + jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i]; for (; jneighnthreads; - const int inum = listinner->inum; + const int inum = list->inum_inner; #if defined(_OPENMP) #pragma omp parallel default(none) #endif @@ -403,7 +403,7 @@ void PairLJLongTIP4PLongOMP::compute_middle() const int nall = atom->nlocal + atom->nghost; const int nthreads = comm->nthreads; - const int inum = listmiddle->inum; + const int inum = list->inum_middle; #if defined(_OPENMP) #pragma omp parallel default(none) @@ -457,7 +457,7 @@ void PairLJLongTIP4PLongOMP::compute_outer(int eflag, int vflag) } const int nthreads = comm->nthreads; - const int inum = listouter->inum; + const int inum = list->inum; #if defined(_OPENMP) #pragma omp parallel default(none) shared(eflag,vflag) @@ -1126,9 +1126,9 @@ void PairLJLongTIP4PLongOMP::eval_inner(int iifrom, int iito, ThrData * const th double *lj1i, *lj2i; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; // loop over neighbors of my atoms @@ -1388,9 +1388,9 @@ void PairLJLongTIP4PLongOMP::eval_middle(int iifrom, int iito, ThrData * const t int ni; double *lj1i, *lj2i; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; // loop over neighbors of my atoms @@ -1656,9 +1656,9 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th double fxtmp,fytmp,fztmp; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; // loop over neighbors of my atoms diff --git a/src/USER-OMP/respa_omp.cpp b/src/USER-OMP/respa_omp.cpp index 738538a209..aa4aa65a4f 100644 --- a/src/USER-OMP/respa_omp.cpp +++ b/src/USER-OMP/respa_omp.cpp @@ -108,6 +108,7 @@ void RespaOMP::setup() domain->box_too_small_check(); modify->setup_pre_neighbor(); neighbor->build(); + modify->setup_post_neighbor(); neighbor->ncalls = 0; // compute all forces @@ -200,6 +201,7 @@ void RespaOMP::setup_minimal(int flag) domain->box_too_small_check(); modify->setup_pre_neighbor(); neighbor->build(); + modify->setup_post_neighbor(); neighbor->ncalls = 0; } @@ -311,6 +313,10 @@ void RespaOMP::recurse(int ilevel) } neighbor->build(); timer->stamp(Timer::NEIGH); + if (modify->n_post_neighbor) { + modify->post_neighbor(); + timer->stamp(Timer::MODIFY); + } } else if (ilevel == 0) { timer->stamp(); comm->forward_comm(); 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/min.cpp b/src/min.cpp index af23629cad..653cac71e6 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -246,6 +246,7 @@ void Min::setup(int flag) domain->box_too_small_check(); modify->setup_pre_neighbor(); neighbor->build(); + modify->setup_post_neighbor(); neighbor->ncalls = 0; // remove these restriction eventually @@ -345,6 +346,7 @@ void Min::setup_minimal(int flag) domain->box_too_small_check(); modify->setup_pre_neighbor(); neighbor->build(); + modify->setup_post_neighbor(); neighbor->ncalls = 0; } @@ -503,12 +505,15 @@ double Min::energy_force(int resetflag) if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); timer->stamp(Timer::COMM); if (modify->n_min_pre_neighbor) { - timer->stamp(); modify->min_pre_neighbor(); timer->stamp(Timer::MODIFY); } neighbor->build(); timer->stamp(Timer::NEIGH); + if (modify->n_min_post_neighbor) { + modify->min_post_neighbor(); + timer->stamp(Timer::MODIFY); + } } ev_set(update->ntimestep); 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..934b9f7d9b 100644 --- a/src/neigh_list.cpp +++ b/src/neigh_list.cpp @@ -40,16 +40,18 @@ 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; + respaouter = 0; + respamiddle = 0; + respainner = 0; copy = 0; copymode = 0; - dnum = 0; // ptrs @@ -60,17 +62,24 @@ NeighList::NeighList(LAMMPS *lmp) : Pointers(lmp) listskip = NULL; listfull = NULL; - listhistory = NULL; - fix_history = NULL; - - respamiddle = 0; - listinner = NULL; - listmiddle = NULL; - fix_bond = NULL; ipage = NULL; - dpage = NULL; + + // extra rRESPA lists + + inum_inner = gnum_inner = 0; + ilist_inner = NULL; + numneigh_inner = NULL; + firstneigh_inner = NULL; + + inum_middle = gnum_middle = 0; + ilist_middle = NULL; + numneigh_middle = NULL; + firstneigh_middle = NULL; + + ipage_inner = NULL; + ipage_middle = NULL; // Kokkos package @@ -92,10 +101,21 @@ NeighList::~NeighList() memory->destroy(ilist); memory->destroy(numneigh); memory->sfree(firstneigh); - memory->sfree(firstdouble); - delete [] ipage; - delete [] dpage; + } + + if (respainner) { + memory->destroy(ilist_inner); + memory->destroy(numneigh_inner); + memory->sfree(firstneigh_inner); + delete [] ipage_inner; + } + + if (respamiddle) { + memory->destroy(ilist_middle); + memory->destroy(numneigh_middle); + memory->sfree(firstneigh_middle); + delete [] ipage_middle; } delete [] iskip; @@ -108,8 +128,7 @@ 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 + respaouter -> set all 3 outer/middle/inner flags bond -> set fix_bond to Fix that made the request ------------------------------------------------------------------------- */ @@ -120,8 +139,11 @@ void NeighList::post_constructor(NeighRequest *nq) occasional = nq->occasional; ghost = nq->ghost; ssa = nq->ssa; + history = nq->history; + respaouter = nq->respaouter; + respamiddle = nq->respamiddle; + respainner = nq->respainner; copy = nq->copy; - dnum = nq->dnum; if (nq->copy) listcopy = neighbor->lists[nq->copylist]; @@ -141,24 +163,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; - listinner = neighbor->lists[nq->respainnerlist]; - } else { - respamiddle = 1; - listmiddle = neighbor->lists[nq->respamiddlelist]; - listinner = neighbor->lists[nq->respainnerlist]; - } - } - if (nq->bond) fix_bond = (Fix *) nq->requestor; } @@ -174,32 +178,29 @@ void NeighList::setup_pages(int pgsize_caller, int oneatom_caller) for (int i = 0; i < nmypage; i++) ipage[i].init(oneatom,pgsize,PGDELTA); - if (dnum) { - dpage = new MyPage[nmypage]; + if (respainner) { + ipage_inner = new MyPage[nmypage]; for (int i = 0; i < nmypage; i++) - dpage[i].init(dnum*oneatom,dnum*pgsize,PGDELTA); - } else dpage = NULL; + ipage_inner[i].init(oneatom,pgsize,PGDELTA); + } + + if (respamiddle) { + ipage_middle = new MyPage[nmypage]; + for (int i = 0; i < nmypage; i++) + ipage_middle[i].init(oneatom,pgsize,PGDELTA); + } } /* ---------------------------------------------------------------------- grow per-atom data to allow for nlocal/nall atoms - for parent lists: - also trigger grow in child list(s) which are not built themselves - history calls grow() in listhistory - respaouter calls grow() in respainner, respamiddle triggered by neighbor list build not called if a copy list ------------------------------------------------------------------------- */ 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,10 +219,25 @@ 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"); + + if (respainner) { + memory->destroy(ilist_inner); + memory->destroy(numneigh_inner); + memory->sfree(firstneigh_inner); + memory->create(ilist_inner,maxatom,"neighlist:ilist_inner"); + memory->create(numneigh_inner,maxatom,"neighlist:numneigh_inner"); + firstneigh_inner = (int **) memory->smalloc(maxatom*sizeof(int *), + "neighlist:firstneigh_inner"); + } + + if (respamiddle) { + memory->destroy(ilist_middle); + memory->destroy(numneigh_middle); + memory->sfree(firstneigh_middle); + memory->create(ilist_middle,maxatom,"neighlist:ilist_middle"); + memory->create(numneigh_middle,maxatom,"neighlist:numneigh_middle"); + firstneigh_middle = (int **) memory->smalloc(maxatom*sizeof(int *), + "neighlist:firstneigh_middle"); } } @@ -253,22 +269,20 @@ void NeighList::print_attributes() printf(" %d = size\n",rq->size); printf(" %d = history\n",rq->history); printf(" %d = granonesided\n",rq->granonesided); - printf(" %d = respainner\n",rq->respainner); - printf(" %d = respamiddle\n",rq->respamiddle); printf(" %d = respaouter\n",rq->respaouter); + printf(" %d = respamiddle\n",rq->respamiddle); + printf(" %d = respainner\n",rq->respainner); printf(" %d = bond\n",rq->bond); printf(" %d = omp\n",rq->omp); printf(" %d = intel\n",rq->intel); 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,10 +306,23 @@ 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(); + if (respainner) { + bytes += memory->usage(ilist_inner,maxatom); + bytes += memory->usage(numneigh_inner,maxatom); + bytes += maxatom * sizeof(int *); + if (ipage_inner) { + for (int i = 0; i < nmypage; i++) + bytes += ipage_inner[i].size(); + } + } + + if (respamiddle) { + bytes += memory->usage(ilist_middle,maxatom); + bytes += memory->usage(numneigh_middle,maxatom); + bytes += maxatom * sizeof(int *); + if (ipage_middle) { + for (int i = 0; i < nmypage; i++) + bytes += ipage_middle[i].size(); } } diff --git a/src/neigh_list.h b/src/neigh_list.h index 4010a68857..d633ba839e 100644 --- a/src/neigh_list.h +++ b/src/neigh_list.h @@ -34,9 +34,12 @@ 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 respaouter; // 1 if list is a rRespa outer list + int respamiddle; // 1 if there is also a rRespa middle list + int respainner; // 1 if there is also a rRespa inner list + 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 +48,28 @@ 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 + + // data structs to store rRESPA neighbor pairs I,J and associated values + + int inum_inner; // # of I atoms neighbors are stored for + int gnum_inner; // # of ghost atoms neighbors are stored for + int *ilist_inner; // local indices of I atoms + int *numneigh_inner; // # of J neighbors for each I atom + int **firstneigh_inner; // ptr to 1st J int value of each I atom + + int inum_middle; // # of I atoms neighbors are stored for + int gnum_middle; // # of ghost atoms neighbors are stored for + int *ilist_middle; // local indices of I atoms + int *numneigh_middle; // # of J neighbors for each I atom + int **firstneigh_middle; // ptr to 1st J int value of each I atom + + MyPage *ipage_inner; // pages of neighbor indices for inner + MyPage *ipage_middle; // pages of neighbor indices for middle // atom types to skip when building list // copied info from corresponding request into realloced vec/array @@ -65,13 +83,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 - class Fix *fix_bond; // fix that stores bond info // Kokkos package diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp index 8d720e766c..6325eec566 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,11 +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; unique = 0; // internal settings @@ -158,8 +151,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 +217,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..16e6f1a8c0 100644 --- a/src/neigh_request.h +++ b/src/neigh_request.h @@ -59,12 +59,12 @@ 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 - int respamiddle; // 1 if a rRESPA middle list - int respaouter; // 1 if a rRESPA outer list + int respainner; // 1 if need a rRESPA inner list + int respamiddle; // 1 if need a rRESPA middle list + int respaouter; // 1 if need a rRESPA outer list int bond; // 1 if store bond neighbors instead of atom neighs int omp; // set by USER-OMP package int intel; // set by USER-INTEL package @@ -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,21 +98,9 @@ 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 - int unique; // 1 if this list requires its own // NStencil, Nbin class - because of requestor cutoff - // pointer to FSH class, set by requestor class (not by Neighbor) - - class FixShearHistory *fix_history; // fix that stores per-atom history info - // ----------------------------- // internal settings made by Neighbor class // ----------------------------- diff --git a/src/neighbor.cpp b/src/neighbor.cpp index a460be0065..cc2e5d6d11 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) unique = create unique lists if cutoff is explicitly set // (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 @@ -681,7 +675,7 @@ int Neighbor::init_pair() int nrequest_original = nrequest; - morph_other(); + morph_unique(); morph_skip(); morph_granular(); // this method can change flags set by requestor morph_halffull(); @@ -827,23 +821,13 @@ 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 do not store: copy // use atom->nmax for both grow() args // i.e. grow first time to expanded size to avoid future reallocs @@ -923,40 +907,16 @@ int Neighbor::init_pair() /* ---------------------------------------------------------------------- scan NeighRequests to set additional flags - only for history, respaouter, custom cutoff lists + only for custom cutoff lists ------------------------------------------------------------------------- */ -void Neighbor::morph_other() +void Neighbor::morph_unique() { NeighRequest *irq; 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) { - if (requests[i-1]->respainner) { - irq->respainnerlist = i-1; - requests[i-1]->respaouterlist = i; - } else { - irq->respamiddlelist = i-1; - requests[i-1]->respaouterlist = i; - requests[i-1]->respainnerlist = i-1; - irq->respainnerlist = i-2; - requests[i-2]->respaouterlist = i; - requests[i-2]->respamiddlelist = i-1; - } - } - // if cut flag set by requestor, set unique flag // this forces Pair,Stencil,Bin styles to be instantiated separately @@ -987,8 +947,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; @@ -1021,12 +979,11 @@ void Neighbor::morph_skip() // else 2 lists do not store same pairs // 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 +1002,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 +1064,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,8 +1115,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; // check all other lists @@ -1179,11 +1133,10 @@ void Neighbor::morph_halffull() // else 2 lists do not store same pairs // 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; @@ -1230,12 +1183,6 @@ void Neighbor::morph_copy() if (irq->copy) continue; - // 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 for (j = 0; j < nrequest; j++) { @@ -1272,9 +1219,9 @@ 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 list with respa middle/inner + // b/c its outer list will not be complete - if (jrq->history) continue; if (jrq->respamiddle) continue; if (jrq->respainner) continue; @@ -1282,12 +1229,11 @@ void Neighbor::morph_copy() // else 2 lists do not store same pairs // 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 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 +1481,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); @@ -1562,9 +1506,8 @@ void Neighbor::print_pairwise_info() if (rq->size) fprintf(out,", size"); if (rq->history) fprintf(out,", history"); if (rq->granonesided) fprintf(out,", onesided"); - if (rq->respainner) fprintf(out,", respa outer"); - if (rq->respamiddle) fprintf(out,", respa middle"); - if (rq->respaouter) fprintf(out,", respa inner"); + if (rq->respamiddle) fprintf(out,", respa outer/middle/inner"); + else if (rq->respainner) fprintf(out,", respa outer/inner"); if (rq->bond) fprintf(out,", bond"); if (rq->omp) fprintf(out,", omp"); if (rq->intel) fprintf(out,", intel"); @@ -1659,8 +1602,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 // checks are bitwise using NeighConst bit masks @@ -1701,8 +1642,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 @@ -1793,11 +1732,6 @@ int Neighbor::choose_stencil(NeighRequest *rq) 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 if (includegroup && rq->ghost) diff --git a/src/neighbor.h b/src/neighbor.h index 64bced2293..9244bc575d 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 @@ -205,7 +204,7 @@ class Neighbor : protected Pointers { int init_pair(); virtual void init_topology(); - void morph_other(); + void morph_unique(); void morph_skip(); void morph_granular(); void morph_halffull(); 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_respa_bin_newtoff.cpp b/src/npair_half_respa_bin_newtoff.cpp index 11246b4af8..0145771f4a 100644 --- a/src/npair_half_respa_bin_newtoff.cpp +++ b/src/npair_half_respa_bin_newtoff.cpp @@ -63,22 +63,19 @@ void NPairHalfRespaBinNewtoff::build(NeighList *list) int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - NeighList *listinner = list->listinner; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; - MyPage *ipage_inner = listinner->ipage; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; + MyPage *ipage_inner = list->ipage_inner; - NeighList *listmiddle; int *ilist_middle,*numneigh_middle,**firstneigh_middle; MyPage *ipage_middle; int respamiddle = list->respamiddle; if (respamiddle) { - listmiddle = list->listmiddle; - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; - ipage_middle = listmiddle->ipage; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; + ipage_middle = list->ipage_middle; } int inum = 0; @@ -185,6 +182,6 @@ void NPairHalfRespaBinNewtoff::build(NeighList *list) } list->inum = inum; - listinner->inum = inum; - if (respamiddle) listmiddle->inum = inum; + list->inum_inner = inum; + if (respamiddle) list->inum_middle = inum; } diff --git a/src/npair_half_respa_bin_newton.cpp b/src/npair_half_respa_bin_newton.cpp index db76678036..72a613204d 100644 --- a/src/npair_half_respa_bin_newton.cpp +++ b/src/npair_half_respa_bin_newton.cpp @@ -62,22 +62,19 @@ void NPairHalfRespaBinNewton::build(NeighList *list) int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - NeighList *listinner = list->listinner; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; - MyPage *ipage_inner = listinner->ipage; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; + MyPage *ipage_inner = list->ipage_inner; - NeighList *listmiddle; int *ilist_middle,*numneigh_middle,**firstneigh_middle; MyPage *ipage_middle; int respamiddle = list->respamiddle; if (respamiddle) { - listmiddle = list->listmiddle; - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; - ipage_middle = listmiddle->ipage; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; + ipage_middle = list->ipage_middle; } int inum = 0; @@ -231,6 +228,6 @@ void NPairHalfRespaBinNewton::build(NeighList *list) } list->inum = inum; - listinner->inum = inum; - if (respamiddle) listmiddle->inum = inum; + list->inum_inner = inum; + if (respamiddle) list->inum_middle = inum; } diff --git a/src/npair_half_respa_bin_newton_tri.cpp b/src/npair_half_respa_bin_newton_tri.cpp index 4ec6685e1d..add1cf6e5c 100644 --- a/src/npair_half_respa_bin_newton_tri.cpp +++ b/src/npair_half_respa_bin_newton_tri.cpp @@ -63,22 +63,19 @@ void NPairHalfRespaBinNewtonTri::build(NeighList *list) int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - NeighList *listinner = list->listinner; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; - MyPage *ipage_inner = listinner->ipage; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; + MyPage *ipage_inner = list->ipage_inner; - NeighList *listmiddle; int *ilist_middle,*numneigh_middle,**firstneigh_middle; MyPage *ipage_middle; int respamiddle = list->respamiddle; if (respamiddle) { - listmiddle = list->listmiddle; - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; - ipage_middle = listmiddle->ipage; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; + ipage_middle = list->ipage_middle; } int inum = 0; @@ -193,6 +190,6 @@ void NPairHalfRespaBinNewtonTri::build(NeighList *list) } list->inum = inum; - listinner->inum = inum; - if (respamiddle) listmiddle->inum = inum; + list->inum_inner = inum; + if (respamiddle) list->inum_middle = inum; } diff --git a/src/npair_half_respa_nsq_newtoff.cpp b/src/npair_half_respa_nsq_newtoff.cpp index 1bb2034384..c0e932f0ae 100644 --- a/src/npair_half_respa_nsq_newtoff.cpp +++ b/src/npair_half_respa_nsq_newtoff.cpp @@ -67,22 +67,19 @@ void NPairHalfRespaNsqNewtoff::build(NeighList *list) int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - NeighList *listinner = list->listinner; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; - MyPage *ipage_inner = listinner->ipage; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; + MyPage *ipage_inner = list->ipage_inner; - NeighList *listmiddle; int *ilist_middle,*numneigh_middle,**firstneigh_middle; MyPage *ipage_middle; int respamiddle = list->respamiddle; if (respamiddle) { - listmiddle = list->listmiddle; - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; - ipage_middle = listmiddle->ipage; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; + ipage_middle = list->ipage_middle; } int inum = 0; @@ -180,6 +177,6 @@ void NPairHalfRespaNsqNewtoff::build(NeighList *list) } list->inum = inum; - listinner->inum = inum; - if (respamiddle) listmiddle->inum = inum; + list->inum_inner = inum; + if (respamiddle) list->inum_middle = inum; } diff --git a/src/npair_half_respa_nsq_newton.cpp b/src/npair_half_respa_nsq_newton.cpp index 9aacc702cc..f7d161896d 100644 --- a/src/npair_half_respa_nsq_newton.cpp +++ b/src/npair_half_respa_nsq_newton.cpp @@ -69,22 +69,19 @@ void NPairHalfRespaNsqNewton::build(NeighList *list) int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - NeighList *listinner = list->listinner; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; - MyPage *ipage_inner = listinner->ipage; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; + MyPage *ipage_inner = list->ipage_inner; - NeighList *listmiddle; int *ilist_middle,*numneigh_middle,**firstneigh_middle; MyPage *ipage_middle; int respamiddle = list->respamiddle; if (respamiddle) { - listmiddle = list->listmiddle; - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; - ipage_middle = listmiddle->ipage; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; + ipage_middle = list->ipage_middle; } int inum = 0; @@ -200,6 +197,6 @@ void NPairHalfRespaNsqNewton::build(NeighList *list) } list->inum = inum; - listinner->inum = inum; - if (respamiddle) listmiddle->inum = inum; + list->inum_inner = inum; + if (respamiddle) list->inum_middle = inum; } 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_respa.cpp b/src/npair_skip_respa.cpp index 31420b32d1..1d4eda5354 100644 --- a/src/npair_skip_respa.cpp +++ b/src/npair_skip_respa.cpp @@ -53,28 +53,24 @@ void NPairSkipRespa::build(NeighList *list) int *iskip = list->iskip; int **ijskip = list->ijskip; - NeighList *listinner = list->listinner; - int *ilist_inner = listinner->ilist; - int *numneigh_inner = listinner->numneigh; - int **firstneigh_inner = listinner->firstneigh; - MyPage *ipage_inner = listinner->ipage; + int *ilist_inner = list->ilist_inner; + int *numneigh_inner = list->numneigh_inner; + int **firstneigh_inner = list->firstneigh_inner; + MyPage *ipage_inner = list->ipage_inner; + int *numneigh_inner_skip = list->listskip->numneigh_inner; + int **firstneigh_inner_skip = list->listskip->firstneigh_inner; - int *numneigh_inner_skip = list->listskip->listinner->numneigh; - int **firstneigh_inner_skip = list->listskip->listinner->firstneigh; - - NeighList *listmiddle; int *ilist_middle,*numneigh_middle,**firstneigh_middle; MyPage *ipage_middle; int *numneigh_middle_skip,**firstneigh_middle_skip; int respamiddle = list->respamiddle; if (respamiddle) { - listmiddle = list->listmiddle; - ilist_middle = listmiddle->ilist; - numneigh_middle = listmiddle->numneigh; - firstneigh_middle = listmiddle->firstneigh; - ipage_middle = listmiddle->ipage; - numneigh_middle_skip = list->listskip->listmiddle->numneigh; - firstneigh_middle_skip = list->listskip->listmiddle->firstneigh; + ilist_middle = list->ilist_middle; + numneigh_middle = list->numneigh_middle; + firstneigh_middle = list->firstneigh_middle; + ipage_middle = list->ipage_middle; + numneigh_middle_skip = list->listskip->numneigh_middle; + firstneigh_middle_skip = list->listskip->firstneigh_middle; } int inum = 0; @@ -164,6 +160,6 @@ void NPairSkipRespa::build(NeighList *list) } list->inum = inum; - listinner->inum = inum; - if (respamiddle) listmiddle->inum = inum; + list->inum_inner = inum; + if (respamiddle) list->inum_middle = 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/pair.h b/src/pair.h index 0f7b0f85b6..eb71e88224 100644 --- a/src/pair.h +++ b/src/pair.h @@ -92,10 +92,6 @@ class Pair : protected Pointers { class NeighList *list; // standard neighbor list used by most pairs class NeighList *listhalf; // half list used by some pairs class NeighList *listfull; // full list used by some pairs - class NeighList *listhistory; // neighbor history list used by some pairs - class NeighList *listinner; // rRESPA lists used by some pairs - class NeighList *listmiddle; - class NeighList *listouter; int allocated; // 0/1 = whether arrays are allocated // public so external driver can check diff --git a/src/pair_lj96_cut.cpp b/src/pair_lj96_cut.cpp index 83fc5bcdda..842b918fe1 100644 --- a/src/pair_lj96_cut.cpp +++ b/src/pair_lj96_cut.cpp @@ -157,10 +157,10 @@ void PairLJ96Cut::compute_inner() double *special_lj = force->special_lj; int newton_pair = force->newton_pair; - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -231,10 +231,10 @@ void PairLJ96Cut::compute_middle() double *special_lj = force->special_lj; int newton_pair = force->newton_pair; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; @@ -318,10 +318,10 @@ void PairLJ96Cut::compute_outer(int eflag, int vflag) double *special_lj = force->special_lj; int newton_pair = force->newton_pair; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; @@ -487,36 +487,23 @@ void PairLJ96Cut::coeff(int narg, char **arg) void PairLJ96Cut::init_style() { - // request regular or rRESPA neighbor lists + // request regular or rRESPA neighbor list int irequest; + int respa = 0; if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { - int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + } - if (respa == 0) irequest = neighbor->request(this,instance_me); - else if (respa == 1) { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } else { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 2; - neighbor->requests[irequest]->respamiddle = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } + irequest = neighbor->request(this,instance_me); - } else irequest = neighbor->request(this,instance_me); + if (respa >= 1) { + neighbor->requests[irequest]->respaouter = 1; + neighbor->requests[irequest]->respainner = 1; + } + if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; // set rRESPA cutoffs @@ -526,19 +513,6 @@ void PairLJ96Cut::init_style() else cut_respa = NULL; } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJ96Cut::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/pair_lj96_cut.h b/src/pair_lj96_cut.h index 6b677c6429..4d6df02127 100644 --- a/src/pair_lj96_cut.h +++ b/src/pair_lj96_cut.h @@ -33,7 +33,6 @@ class PairLJ96Cut : public Pair { 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 *); diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index 7f838061f1..215fabecbb 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -156,10 +156,10 @@ void PairLJCut::compute_inner() double *special_lj = force->special_lj; int newton_pair = force->newton_pair; - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -229,10 +229,10 @@ void PairLJCut::compute_middle() double *special_lj = force->special_lj; int newton_pair = force->newton_pair; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; @@ -315,10 +315,10 @@ void PairLJCut::compute_outer(int eflag, int vflag) double *special_lj = force->special_lj; int newton_pair = force->newton_pair; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; @@ -481,36 +481,23 @@ void PairLJCut::coeff(int narg, char **arg) void PairLJCut::init_style() { - // request regular or rRESPA neighbor lists + // request regular or rRESPA neighbor list int irequest; + int respa = 0; if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { - int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + } - if (respa == 0) irequest = neighbor->request(this,instance_me); - else if (respa == 1) { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } else { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 2; - neighbor->requests[irequest]->respamiddle = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } + irequest = neighbor->request(this,instance_me); - } else irequest = neighbor->request(this,instance_me); + if (respa >= 1) { + neighbor->requests[irequest]->respaouter = 1; + neighbor->requests[irequest]->respainner = 1; + } + if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; // set rRESPA cutoffs @@ -520,19 +507,6 @@ void PairLJCut::init_style() else cut_respa = NULL; } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJCut::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/pair_lj_cut.h b/src/pair_lj_cut.h index 43eeda09ca..3724685db6 100644 --- a/src/pair_lj_cut.h +++ b/src/pair_lj_cut.h @@ -32,7 +32,6 @@ class PairLJCut : public Pair { 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 *); diff --git a/src/pair_mie_cut.cpp b/src/pair_mie_cut.cpp index 320f21248d..04f8de8d7d 100644 --- a/src/pair_mie_cut.cpp +++ b/src/pair_mie_cut.cpp @@ -159,10 +159,10 @@ void PairMIECut::compute_inner() double *special_mie = force->special_lj; int newton_pair = force->newton_pair; - inum = listinner->inum; - ilist = listinner->ilist; - numneigh = listinner->numneigh; - firstneigh = listinner->firstneigh; + inum = list->inum_inner; + ilist = list->ilist_inner; + numneigh = list->numneigh_inner; + firstneigh = list->firstneigh_inner; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -233,10 +233,10 @@ void PairMIECut::compute_middle() double *special_mie = force->special_lj; int newton_pair = force->newton_pair; - inum = listmiddle->inum; - ilist = listmiddle->ilist; - numneigh = listmiddle->numneigh; - firstneigh = listmiddle->firstneigh; + inum = list->inum_middle; + ilist = list->ilist_middle; + numneigh = list->numneigh_middle; + firstneigh = list->firstneigh_middle; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; @@ -320,10 +320,10 @@ void PairMIECut::compute_outer(int eflag, int vflag) double *special_mie = force->special_lj; int newton_pair = force->newton_pair; - inum = listouter->inum; - ilist = listouter->ilist; - numneigh = listouter->numneigh; - firstneigh = listouter->firstneigh; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; @@ -496,36 +496,23 @@ void PairMIECut::coeff(int narg, char **arg) void PairMIECut::init_style() { - // request regular or rRESPA neighbor lists + // request regular or rRESPA neighbor list int irequest; + int respa = 0; if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { - int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + } - if (respa == 0) irequest = neighbor->request(this,instance_me); - else if (respa == 1) { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } else { - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 2; - neighbor->requests[irequest]->respamiddle = 1; - irequest = neighbor->request(this,instance_me); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->respaouter = 1; - } + irequest = neighbor->request(this,instance_me); - } else irequest = neighbor->request(this,instance_me); + if (respa >= 1) { + neighbor->requests[irequest]->respaouter = 1; + neighbor->requests[irequest]->respainner = 1; + } + if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; // set rRESPA cutoffs @@ -535,19 +522,6 @@ void PairMIECut::init_style() else cut_respa = NULL; } -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairMIECut::init_list(int id, NeighList *ptr) -{ - if (id == 0) list = ptr; - else if (id == 1) listinner = ptr; - else if (id == 2) listmiddle = ptr; - else if (id == 3) listouter = ptr; -} - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ diff --git a/src/pair_mie_cut.h b/src/pair_mie_cut.h index 2a0a29843e..9e12438d14 100644 --- a/src/pair_mie_cut.h +++ b/src/pair_mie_cut.h @@ -32,7 +32,6 @@ class PairMIECut : public Pair { 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 *); diff --git a/src/respa.cpp b/src/respa.cpp index 5d51ff64ee..23cd941834 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -442,6 +442,7 @@ void Respa::setup(int flag) domain->box_too_small_check(); modify->setup_pre_neighbor(); neighbor->build(); + modify->setup_post_neighbor(); neighbor->ncalls = 0; // compute all forces @@ -517,6 +518,7 @@ void Respa::setup_minimal(int flag) domain->box_too_small_check(); modify->setup_pre_neighbor(); neighbor->build(); + modify->setup_post_neighbor(); neighbor->ncalls = 0; } @@ -668,6 +670,11 @@ void Respa::recurse(int ilevel) } neighbor->build(); timer->stamp(Timer::NEIGH); + if (modify->n_post_neighbor) { + modify->post_neighbor(); + timer->stamp(Timer::MODIFY); + } + } else if (ilevel == 0) { timer->stamp(); comm->forward_comm(); 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