jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #688 from lammps/history 

refactoring of neighbor history
This commit is contained in:
Steve Plimpton 2017-10-18 13:24:36 -06:00 committed by GitHub
parent f8f13d929f 17c17ac409
commit 7dfc6b7eab
87 changed files with 1136 additions and 2081 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
src/CORESHELL/pair_lj_cut_coul_long_cs.cpp
View File

@ -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];

16
src/DIPOLE/pair_lj_long_dipole_long.cpp
View File

@ -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
------------------------------------------------------------------------- */

1
src/DIPOLE/pair_lj_long_dipole_long.h
View File

@ -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 *);

7
src/GRANULAR/pair_gran_hertz_history.cpp
View File

@ -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++;

57
src/GRANULAR/pair_gran_hooke_history.cpp
View File

@ -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;
}

4
src/GRANULAR/pair_gran_hooke_history.h
View File

@ -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

9
src/KOKKOS/neigh_list_kokkos.cpp
View File

@ -49,15 +49,6 @@ void NeighListKokkos<Device>::grow(int nmax)
d_neighbors =
typename ArrayTypes<Device>::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");
}
/* ---------------------------------------------------------------------- */

3
src/KOKKOS/npair_copy_kokkos.cpp
View File

@ -41,10 +41,7 @@ void NPairCopyKokkos<DeviceType>::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<DeviceType>* list_kk = (NeighListKokkos<DeviceType>*) list;
NeighListKokkos<DeviceType>* listcopy_kk = (NeighListKokkos<DeviceType>*) list->listcopy;

57
src/KSPACE/pair_buck_long_coul_long.cpp
View File

@ -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;
}
} 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 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 (; ineigh<ineighn; ++ineigh) { // loop over my atoms
i = *ineigh; fi = f0+3*i;
if (order1) qri = qqrd2e*q[i];
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 (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;
@ -740,7 +715,7 @@ void PairBuckLongCoulLong::compute_middle()
double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi;
vector xi, d;
ineighn = (ineigh = listmiddle->ilist)+listmiddle->inum;
ineighn = (ineigh = list->ilist_middle)+list->inum_middle;
for (; ineigh<ineighn; ++ineigh) { // loop over my atoms
i = *ineigh; fi = f0+3*i;
@ -748,7 +723,7 @@ void PairBuckLongCoulLong::compute_middle()
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 (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;
@ -839,7 +814,7 @@ void PairBuckLongCoulLong::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;
ineighn = (ineigh = listouter->ilist)+listouter->inum;
ineighn = (ineigh = list->ilist)+list->inum;
for (; ineigh<ineighn; ++ineigh) { // loop over my atoms
i = *ineigh; fi = f0+3*i;
@ -849,7 +824,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
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<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;

1
src/KSPACE/pair_buck_long_coul_long.h
View File

@ -35,7 +35,6 @@ class PairBuckLongCoulLong : 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 *);

68
src/KSPACE/pair_lj_charmm_coul_long.cpp
View File

@ -242,10 +242,10 @@ void PairLJCharmmCoulLong::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;
// 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;
}
} 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;
// 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
------------------------------------------------------------------------- */

1
src/KSPACE/pair_lj_charmm_coul_long.h
View File

@ -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 *);

8
src/KSPACE/pair_lj_charmm_coul_msm.cpp
View File

@ -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];

37
src/KSPACE/pair_lj_charmmfsw_coul_long.cpp
View File

@ -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
------------------------------------------------------------------------- */

1
src/KSPACE/pair_lj_charmmfsw_coul_long.h
View File

@ -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 *);

68
src/KSPACE/pair_lj_cut_coul_long.cpp
View File

@ -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;
}
} 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;
@ -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
------------------------------------------------------------------------- */

1
src/KSPACE/pair_lj_cut_coul_long.h
View File

@ -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 *);

8
src/KSPACE/pair_lj_cut_coul_msm.cpp
View File

@ -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];

54
src/KSPACE/pair_lj_long_coul_long.cpp
View File

@ -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 (; ineigh<ineighn; ++ineigh) { // loop over my atoms
i = *ineigh; fi = f0+3*i;
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 (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;
ni = sbmask(j);
@ -736,7 +710,7 @@ void PairLJLongCoulLong::compute_middle()
double qri, *cut_ljsqi, *lj1i, *lj2i;
vector xi, d;
ineighn = (ineigh = listmiddle->ilist)+listmiddle->inum;
ineighn = (ineigh = list->ilist_middle)+list->inum_middle;
for (; ineigh<ineighn; ++ineigh) { // loop over my atoms
i = *ineigh; fi = f0+3*i;
@ -744,7 +718,7 @@ void PairLJLongCoulLong::compute_middle()
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 (; jneigh<jneighn; ++jneigh) {
j = *jneigh;
@ -833,7 +807,7 @@ void PairLJLongCoulLong::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;
ineighn = (ineigh = listouter->ilist)+listouter->inum;
ineighn = (ineigh = list->ilist)+list->inum;
for (; ineigh<ineighn; ++ineigh) { // loop over my atoms
i = *ineigh; fi = f0+3*i;
@ -842,7 +816,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
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 (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;

1
src/KSPACE/pair_lj_long_coul_long.h
View File

@ -34,7 +34,6 @@ class PairLJLongCoulLong : 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 *);

24
src/KSPACE/pair_lj_long_tip4p_long.cpp
View File

@ -516,10 +516,10 @@ void PairLJLongTIP4PLong::compute_inner()
int ni;
double *lj1i, *lj2i;
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
@ -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

4
src/OPT/pair_lj_long_coul_long_opt.cpp
View File

@ -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 (; ineigh<ineighn; ++ineigh) { // loop over my atoms
i = *ineigh; fi = f0+3*i;
@ -735,7 +735,7 @@ void PairLJLongCoulLongOpt::eval_outer()
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 (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;

3
src/Purge.list
View File

@ -16,6 +16,9 @@ style_region.h
style_neigh_bin.h
style_neigh_pair.h
style_neigh_stencil.h
# deleted on 11 October 2017
fix_shear_history_omp.cpp
fix_shear_history_omp.h
# deleted on 5 September 2017
npair_halffull_newton_ssa.cpp
npair_halffull_newton_ssa.h

13
src/USER-DRUDE/pair_lj_cut_thole_long.cpp
View File

@ -378,19 +378,6 @@ void PairLJCutTholeLong::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 PairLJCutTholeLong::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
------------------------------------------------------------------------- */

1
src/USER-DRUDE/pair_lj_cut_thole_long.h
View File

@ -34,7 +34,6 @@ class PairLJCutTholeLong : 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 *);

37
src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp
View File

@ -226,10 +226,10 @@ void PairLJCharmmCoulLongSoft::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];
@ -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
------------------------------------------------------------------------- */

1
src/USER-FEP/pair_lj_charmm_coul_long_soft.h
View File

@ -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 *);

37
src/USER-FEP/pair_lj_cut_coul_long_soft.cpp
View File

@ -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
------------------------------------------------------------------------- */

1
src/USER-FEP/pair_lj_cut_coul_long_soft.h
View File

@ -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 *);

37
src/USER-FEP/pair_lj_cut_soft.cpp
View File

@ -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
------------------------------------------------------------------------- */

1
src/USER-FEP/pair_lj_cut_soft.h
View File

@ -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 *);

79
src/USER-OMP/fix_shear_history_omp.cpp → src/USER-OMP/fix_neigh_history_omp.cpp
View File

@ -13,7 +13,7 @@
#include <string.h>
#include <stdio.h>
#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 <tagint> &ipg = ipage[tid];
MyPage <double> &dpg = dpage[tid];
MyPage <tagint> &ipg = ipage_atom[tid];
MyPage <double> &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);
}
}

13
src/USER-OMP/fix_shear_history_omp.h → src/USER-OMP/fix_neigh_history_omp.h
View File

@ -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();
};

24
src/USER-OMP/npair_half_respa_bin_newtoff_omp.cpp
View File

@ -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<int> &ipage = list->ipage[tid];
MyPage<int> &ipage_inner = listinner->ipage[tid];
MyPage<int> &ipage_inner = list->ipage_inner[tid];
ipage.reset();
ipage_inner.reset();
MyPage<int> *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;
}

24
src/USER-OMP/npair_half_respa_bin_newton_omp.cpp
View File

@ -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<int> &ipage = list->ipage[tid];
MyPage<int> &ipage_inner = listinner->ipage[tid];
MyPage<int> &ipage_inner = list->ipage_inner[tid];
ipage.reset();
ipage_inner.reset();
MyPage<int> *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;
}

24
src/USER-OMP/npair_half_respa_bin_newton_tri_omp.cpp
View File

@ -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<int> &ipage = list->ipage[tid];
MyPage<int> &ipage_inner = listinner->ipage[tid];
MyPage<int> &ipage_inner = list->ipage_inner[tid];
ipage.reset();
ipage_inner.reset();
MyPage<int> *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;
}

24
src/USER-OMP/npair_half_respa_nsq_newtoff_omp.cpp
View File

@ -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<int> &ipage = list->ipage[tid];
MyPage<int> &ipage_inner = listinner->ipage[tid];
MyPage<int> &ipage_inner = list->ipage_inner[tid];
ipage.reset();
ipage_inner.reset();
MyPage<int> *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;
}

24
src/USER-OMP/npair_half_respa_nsq_newton_omp.cpp
View File

@ -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<int> &ipage = list->ipage[tid];
MyPage<int> &ipage_inner = listinner->ipage[tid];
MyPage<int> &ipage_inner = list->ipage_inner[tid];
ipage.reset();
ipage_inner.reset();
MyPage<int> *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;
}

77
src/USER-OMP/npair_half_size_bin_newtoff_omp.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *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<int> &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;

110
src/USER-OMP/npair_half_size_bin_newton_omp.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *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<int> &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;

12
src/USER-OMP/npair_half_size_bin_newton_tri_omp.cpp
View File

@ -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;
}
}
}

80
src/USER-OMP/npair_half_size_nsq_newtoff_omp.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *dpage_shear;
int *neighptr;
double **x = atom->x;
double *radius = atom->radius;
@ -89,29 +72,10 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list)
MyPage<int> &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;

79
src/USER-OMP/npair_half_size_nsq_newton_omp.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *dpage_shear;
int *neighptr;
double **x = atom->x;
double *radius = atom->radius;
@ -90,29 +73,10 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list)
MyPage<int> &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;

18
src/USER-OMP/pair_buck_long_coul_long_omp.cpp
View File

@ -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 (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;
@ -904,7 +904,7 @@ void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const
const double *x0 = x[0];
double *f0 = f[0], *fi = 0;
int *ilist = listmiddle->ilist;
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 (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;
@ -1009,7 +1009,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
const double *x0 = x[0];
double *f0 = f[0], *fi = f0;
int *ilist = listouter->ilist;
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<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;

5
src/USER-OMP/pair_gran_hertz_history_omp.cpp
View File

@ -14,6 +14,7 @@
#include <math.h>
#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

5
src/USER-OMP/pair_gran_hooke_history_omp.cpp
View File

@ -14,6 +14,7 @@
#include <math.h>
#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

18
src/USER-OMP/pair_lj_long_coul_long_omp.cpp
View File

@ -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 (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;
ni = sbmask(j);
@ -896,7 +896,7 @@ void PairLJLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const th
const double *x0 = x[0];
double *f0 = f[0], *fi = 0;
int *ilist = listmiddle->ilist;
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 (; jneigh<jneighn; ++jneigh) {
j = *jneigh;
@ -1000,7 +1000,7 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
const double *x0 = x[0];
double *f0 = f[0], *fi = f0;
int *ilist = listouter->ilist;
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 (; jneigh<jneighn; ++jneigh) { // loop over neighbors
j = *jneigh;

24
src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
View File

@ -379,7 +379,7 @@ void PairLJLongTIP4PLongOMP::compute_inner()
for (i = 0; i < nall; i++) hneigh_thr[i].t = 0;
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
@ -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

6
src/USER-OMP/respa_omp.cpp
View File

@ -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();

43
src/fix.h
View File

@ -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;
}
}

498
src/fix_shear_history.cpp → src/fix_neigh_history.cpp
View File

@ -14,7 +14,7 @@
#include <mpi.h>
#include <string.h>
#include <stdio.h>
#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<tagint>[nmypage];
dpage = new MyPage<double>[nmypage];
ipage_atom = new MyPage<tagint>[nmypage];
dpage_atom = new MyPage<double>[nmypage];
ipage_neigh = new MyPage<int>[nmypage];
dpage_neigh = new MyPage<double>[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<tagint> (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<int> (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<tagint> (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<int> (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<int> (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<tagint> (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;
}

55
src/fix_shear_history.h → src/fix_neigh_history.h
View File

@ -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<tagint> *ipage; // pages of partner atom IDs
MyPage<double> *dpage; // pages of shear history with partners
MyPage<tagint> *ipage_atom; // pages of partner atom IDs
MyPage<double> *dpage_atom; // pages of partner values
// per-neighbor data structs pointed to by firstflag & firstvalue
MyPage<int> *ipage_neigh; // pages of local atom indices
MyPage<double> *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;
}
};
}

7
src/min.cpp
View File

@ -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);

45
src/modify.cpp
View File

@ -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
------------------------------------------------------------------------- */

12
src/modify.h
View File

@ -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;

147
src/neigh_list.cpp
View File

@ -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<double>[nmypage];
if (respainner) {
ipage_inner = new MyPage<int>[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<int>[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();
}
}

33
src/neigh_list.h
View File

@ -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<int> *ipage; // pages of neighbor indices
MyPage<double> *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<int> *ipage_inner; // pages of neighbor indices for inner
MyPage<int> *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

13
src/neigh_request.cpp
View File

@ -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;

22
src/neigh_request.h
View File

@ -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
// -----------------------------

94
src/neighbor.cpp
View File

@ -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)

3
src/neighbor.h
View File

@ -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();

1
src/npair.cpp
View File

@ -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;

1
src/npair.h
View File

@ -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

2
src/npair_copy.cpp
View File

@ -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;
}

23
src/npair_half_respa_bin_newtoff.cpp
View File

@ -63,22 +63,19 @@ void NPairHalfRespaBinNewtoff::build(NeighList *list)
int **firstneigh = list->firstneigh;
MyPage<int> *ipage = list->ipage;
NeighList *listinner = list->listinner;
int *ilist_inner = listinner->ilist;
int *numneigh_inner = listinner->numneigh;
int **firstneigh_inner = listinner->firstneigh;
MyPage<int> *ipage_inner = listinner->ipage;
int *ilist_inner = list->ilist_inner;
int *numneigh_inner = list->numneigh_inner;
int **firstneigh_inner = list->firstneigh_inner;
MyPage<int> *ipage_inner = list->ipage_inner;
NeighList *listmiddle;
int *ilist_middle,*numneigh_middle,**firstneigh_middle;
MyPage<int> *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;
}

23
src/npair_half_respa_bin_newton.cpp
View File

@ -62,22 +62,19 @@ void NPairHalfRespaBinNewton::build(NeighList *list)
int **firstneigh = list->firstneigh;
MyPage<int> *ipage = list->ipage;
NeighList *listinner = list->listinner;
int *ilist_inner = listinner->ilist;
int *numneigh_inner = listinner->numneigh;
int **firstneigh_inner = listinner->firstneigh;
MyPage<int> *ipage_inner = listinner->ipage;
int *ilist_inner = list->ilist_inner;
int *numneigh_inner = list->numneigh_inner;
int **firstneigh_inner = list->firstneigh_inner;
MyPage<int> *ipage_inner = list->ipage_inner;
NeighList *listmiddle;
int *ilist_middle,*numneigh_middle,**firstneigh_middle;
MyPage<int> *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;
}

23
src/npair_half_respa_bin_newton_tri.cpp
View File

@ -63,22 +63,19 @@ void NPairHalfRespaBinNewtonTri::build(NeighList *list)
int **firstneigh = list->firstneigh;
MyPage<int> *ipage = list->ipage;
NeighList *listinner = list->listinner;
int *ilist_inner = listinner->ilist;
int *numneigh_inner = listinner->numneigh;
int **firstneigh_inner = listinner->firstneigh;
MyPage<int> *ipage_inner = listinner->ipage;
int *ilist_inner = list->ilist_inner;
int *numneigh_inner = list->numneigh_inner;
int **firstneigh_inner = list->firstneigh_inner;
MyPage<int> *ipage_inner = list->ipage_inner;
NeighList *listmiddle;
int *ilist_middle,*numneigh_middle,**firstneigh_middle;
MyPage<int> *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;
}

23
src/npair_half_respa_nsq_newtoff.cpp
View File

@ -67,22 +67,19 @@ void NPairHalfRespaNsqNewtoff::build(NeighList *list)
int **firstneigh = list->firstneigh;
MyPage<int> *ipage = list->ipage;
NeighList *listinner = list->listinner;
int *ilist_inner = listinner->ilist;
int *numneigh_inner = listinner->numneigh;
int **firstneigh_inner = listinner->firstneigh;
MyPage<int> *ipage_inner = listinner->ipage;
int *ilist_inner = list->ilist_inner;
int *numneigh_inner = list->numneigh_inner;
int **firstneigh_inner = list->firstneigh_inner;
MyPage<int> *ipage_inner = list->ipage_inner;
NeighList *listmiddle;
int *ilist_middle,*numneigh_middle,**firstneigh_middle;
MyPage<int> *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;
}

23
src/npair_half_respa_nsq_newton.cpp
View File

@ -69,22 +69,19 @@ void NPairHalfRespaNsqNewton::build(NeighList *list)
int **firstneigh = list->firstneigh;
MyPage<int> *ipage = list->ipage;
NeighList *listinner = list->listinner;
int *ilist_inner = listinner->ilist;
int *numneigh_inner = listinner->numneigh;
int **firstneigh_inner = listinner->firstneigh;
MyPage<int> *ipage_inner = listinner->ipage;
int *ilist_inner = list->ilist_inner;
int *numneigh_inner = list->numneigh_inner;
int **firstneigh_inner = list->firstneigh_inner;
MyPage<int> *ipage_inner = list->ipage_inner;
NeighList *listmiddle;
int *ilist_middle,*numneigh_middle,**firstneigh_middle;
MyPage<int> *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;
}

78
src/npair_half_size_bin_newtoff.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *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<int> *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];
@ -125,29 +89,10 @@ void NPairHalfSizeBinNewtoff::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;
}
}
}
@ -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;

105
src/npair_half_size_bin_newton.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *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<int> *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;

76
src/npair_half_size_bin_newton_tri.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *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<int> *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) {
if (history && rsq < radsum*radsum)
neighptr[n++] = j ^ mask_history;
else
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++;
}
}
}
@ -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;

78
src/npair_half_size_nsq_newtoff.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *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<int> *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;

78
src/npair_half_size_nsq_newton.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *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<int> *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;

32
src/npair_skip_respa.cpp
View File

@ -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<int> *ipage_inner = listinner->ipage;
int *ilist_inner = list->ilist_inner;
int *numneigh_inner = list->numneigh_inner;
int **firstneigh_inner = list->firstneigh_inner;
MyPage<int> *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<int> *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;
}

74
src/npair_skip_size.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *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;

71
src/npair_skip_size_off2on.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *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;

59
src/npair_skip_size_off2on_oneside.cpp
View File

@ -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<int> *ipage_touch;
MyPage<double> *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;
}

4
src/pair.h
View File

@ -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

68
src/pair_lj96_cut.cpp
View File

@ -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;
}
} 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;
// 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
------------------------------------------------------------------------- */

1
src/pair_lj96_cut.h
View File

@ -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 *);

68
src/pair_lj_cut.cpp
View File

@ -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;
}
} 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;
// 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
------------------------------------------------------------------------- */

1
src/pair_lj_cut.h
View File

@ -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 *);

68
src/pair_mie_cut.cpp
View File

@ -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;
}
} 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;
// 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
------------------------------------------------------------------------- */

1
src/pair_mie_cut.h
View File

@ -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 *);

7
src/respa.cpp
View File

@ -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();

10
src/verlet.cpp
View File

@ -91,8 +91,7 @@ void Verlet::setup(int flag)
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," 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