From 159d131e37cd697a4506053f7d2cc545e8e8477a Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 18 Mar 2011 15:09:03 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5797 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/COLLOID/pair_colloid.cpp | 4 +- src/MANYBODY/pair_airebo.cpp | 119 +++-------------- src/MANYBODY/pair_airebo.h | 2 + src/PERI/pair_peri_lps.cpp | 5 - src/PERI/pair_peri_pmb.cpp | 5 - src/USER-REAXC/fix_reax_c.cpp | 2 +- src/USER-REAXC/reaxc_types.h | 1 - src/atom.cpp | 9 +- src/comm.cpp | 6 +- src/neigh_derive.cpp | 16 ++- src/neigh_full.cpp | 241 +++++++++++++++++++++++++++++++++- src/neigh_list.cpp | 44 +++++-- src/neigh_list.h | 7 +- src/neigh_request.cpp | 8 ++ src/neigh_request.h | 6 +- src/neigh_stencil.cpp | 44 +++++++ src/neighbor.cpp | 111 ++++++++++++---- src/neighbor.h | 9 ++ src/pair.cpp | 2 + src/pair.h | 4 +- src/pair_hybrid.cpp | 15 ++- 21 files changed, 493 insertions(+), 167 deletions(-) diff --git a/src/COLLOID/pair_colloid.cpp b/src/COLLOID/pair_colloid.cpp index b67463432d..7af1410ed8 100644 --- a/src/COLLOID/pair_colloid.cpp +++ b/src/COLLOID/pair_colloid.cpp @@ -348,7 +348,6 @@ double PairColloid::init_one(int i, int j) sigma6[j][i] = sigma6[i][j]; diameter[j][i] = diameter[i][j]; cut[j][i] = cut[i][j]; - cutsq[j][i] = cutsq[i][j] = cut[i][j] * cut[i][j]; double epsilon = a12[i][j]/144.0; lj1[j][i] = lj1[i][j] = 48.0 * epsilon * sigma6[i][j] * sigma6[i][j]; @@ -359,7 +358,8 @@ double PairColloid::init_one(int i, int j) offset[j][i] = offset[i][j] = 0.0; if (offset_flag) { double tmp; - offset[j][i] = offset[i][j] = single(0,0,i,j,cutsq[i][j],0.0,1.0,tmp); + offset[j][i] = offset[i][j] = + single(0,0,i,j,cut[i][j]*cut[i][j],0.0,1.0,tmp); } return cut[i][j]; diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 37081d05de..a60c19521a 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -48,14 +48,17 @@ PairAIREBO::PairAIREBO(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; one_coeff = 1; + ghostneigh = 1; maxlocal = 0; REBO_numneigh = NULL; REBO_firstneigh = NULL; - closestdistsq = NULL; maxpage = 0; pages = NULL; nC = nH = NULL; + + time1 = 0.0; + time2 = 0.0; } /* ---------------------------------------------------------------------- @@ -68,13 +71,13 @@ PairAIREBO::~PairAIREBO() memory->sfree(REBO_firstneigh); for (int i = 0; i < maxpage; i++) memory->sfree(pages[i]); memory->sfree(pages); - memory->sfree(closestdistsq); memory->sfree(nC); memory->sfree(nH); if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); + memory->destroy_2d_double_array(cutghost); memory->destroy_2d_double_array(cutljsq); memory->destroy_2d_double_array(lj1); @@ -92,12 +95,18 @@ void PairAIREBO::compute(int eflag, int vflag) if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = vflag_atom = 0; + double tstart = MPI_Wtime(); REBO_neigh(); + double tmid = MPI_Wtime(); FREBO(eflag,vflag); if (ljflag) FLJ(eflag,vflag); if (torflag) TORSION(eflag,vflag); if (vflag_fdotr) virial_compute(); + double tstop = MPI_Wtime(); + + time1 += tmid-tstart; + time2 += tstop-tmid; } /* ---------------------------------------------------------------------- @@ -115,6 +124,7 @@ void PairAIREBO::allocate() setflag[i][j] = 0; cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + cutghost = memory->create_2d_double_array(n+1,n+1,"pair:cutghost"); // only sized by C,H = 2 types @@ -210,11 +220,12 @@ void PairAIREBO::init_style() if (force->newton_pair == 0) error->all("Pair style AIREBO requires newton pair on"); - // need a full neighbor list + // need a full neighbor list, including neighbors of ghosts int irequest = neighbor->request(this); neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->full = 1; + neighbor->requests[irequest]->ghost = 1; // local REBO neighbor list memory @@ -261,6 +272,7 @@ double PairAIREBO::init_one(int i, int j) // rcLJmax + 2*rcmax, since I-J < rcLJmax and J-L,L-N = REBO distances // long interaction = I-J with I = owned and J = ghost // cutlj*sigma, since I-J < LJ cutoff + // cutghost = REBO cutoff used in REBO_neigh() for neighbors of ghosts double cutmax = cut3rebo; if (ljflag) { @@ -268,14 +280,14 @@ double PairAIREBO::init_one(int i, int j) cutmax = MAX(cutmax,cutlj*sigma[0][0]); } - cutsq[i][j] = cutmax * cutmax; + cutghost[i][j] = rcmax[ii][jj]; cutljsq[ii][jj] = cutlj*sigma[ii][jj] * cutlj*sigma[ii][jj]; lj1[ii][jj] = 48.0 * epsilon[ii][jj] * pow(sigma[ii][jj],12.0); lj2[ii][jj] = 24.0 * epsilon[ii][jj] * pow(sigma[ii][jj],6.0); lj3[ii][jj] = 4.0 * epsilon[ii][jj] * pow(sigma[ii][jj],12.0); lj4[ii][jj] = 4.0 * epsilon[ii][jj] * pow(sigma[ii][jj],6.0); - cutsq[j][i] = cutsq[i][j]; + cutghost[j][i] = cutghost[i][j]; cutljsq[jj][ii] = cutljsq[ii][jj]; lj1[jj][ii] = lj1[ii][jj]; lj2[jj][ii] = lj2[ii][jj]; @@ -292,7 +304,7 @@ double PairAIREBO::init_one(int i, int j) void PairAIREBO::REBO_neigh() { - int i,j,ii,jj,m,n,inum,jnum,itype,jtype; + int i,j,ii,jj,m,n,allnum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS; int *ilist,*jlist,*numneigh,**firstneigh; int *neighptr; @@ -306,41 +318,28 @@ void PairAIREBO::REBO_neigh() maxlocal = atom->nmax; memory->sfree(REBO_numneigh); memory->sfree(REBO_firstneigh); - memory->sfree(closestdistsq); memory->sfree(nC); memory->sfree(nH); REBO_numneigh = (int *) memory->smalloc(maxlocal*sizeof(int),"AIREBO:numneigh"); REBO_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *),"AIREBO:firstneigh"); - closestdistsq = (double *) - memory->smalloc(maxlocal*sizeof(double),"AIREBO:closestdistsq"); nC = (double *) memory->smalloc(maxlocal*sizeof(double),"AIREBO:nC"); nH = (double *) memory->smalloc(maxlocal*sizeof(double),"AIREBO:nH"); } - inum = list->inum; + allnum = list->inum + list->gnum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - // initialize ghost atom references to -1 and closest distance = cutljrebosq - - for (i = nlocal; i < nall; i++) { - REBO_numneigh[i] = -1; - closestdistsq[i] = cutljrebosq; - } - - // store all REBO neighs of owned atoms + // store all REBO neighs of owned and ghost atoms // scan full neighbor list of I - // if J is ghost and within LJ cutoff: - // flag it via REBO_numneigh so its REBO neighbors will be stored below - // REBO requires neighbors of neighbors of i,j in each i,j LJ interaction npage = 0; int npnt = 0; - for (ii = 0; ii < inum; ii++) { + for (ii = 0; ii < allnum; ii++) { i = ilist[ii]; if (pgsize - npnt < oneatom) { @@ -367,82 +366,6 @@ void PairAIREBO::REBO_neigh() delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; - if (rsq < rcmaxsq[itype][jtype]) { - neighptr[n++] = j; - if (jtype == 0) - nC[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS); - else - nH[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS); - } - if (j >= nlocal && rsq < closestdistsq[j]) { - REBO_numneigh[j] = i; - closestdistsq[j] = rsq; - } - } - - REBO_firstneigh[i] = neighptr; - REBO_numneigh[i] = n; - npnt += n; - if (npnt >= pgsize) - error->one("Neighbor list overflow, boost neigh_modify one or page"); - } - - // store REBO neighs of ghost atoms I that have been flagged in REBO_numneigh - // find by scanning full neighbor list of owned atom M = closest neigh of I - - for (i = nlocal; i < nall; i++) { - - if (pgsize - npnt < oneatom) { - npnt = 0; - npage++; - if (npage == maxpage) add_pages(npage); - } - neighptr = &pages[npage][npnt]; - n = 0; - - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = map[type[i]]; - nC[i] = nH[i] = 0.0; - m = REBO_numneigh[i]; - if (m < 0) { - REBO_firstneigh[i] = neighptr; - REBO_numneigh[i] = 0; - continue; - } - - jtype = map[type[m]]; - delx = xtmp - x[m][0]; - dely = ytmp - x[m][1]; - delz = ztmp - x[m][2]; - rsq = delx*delx + dely*dely + delz*delz; - - // add M as neigh of I if close enough - - if (rsq < rcmaxsq[itype][jtype]) { - neighptr[n++] = m; - if (jtype == 0) - nC[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS); - else - nH[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS); - } - - jlist = firstneigh[m]; - jnum = numneigh[m]; - - // add M's neighbors as neighs of I if close enough - // skip I itself - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - if (j == i) continue; - jtype = map[type[j]]; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - if (rsq < rcmaxsq[itype][jtype]) { neighptr[n++] = j; if (jtype == 0) diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h index 09c632f6be..4ee0c7797e 100644 --- a/src/MANYBODY/pair_airebo.h +++ b/src/MANYBODY/pair_airebo.h @@ -36,6 +36,8 @@ class PairAIREBO : public Pair { double memory_usage(); private: + double time1,time2; + int me; int ljflag,torflag; // 0/1 if LJ,torsion terms included int maxlocal; // size of numneigh, firstneigh arrays diff --git a/src/PERI/pair_peri_lps.cpp b/src/PERI/pair_peri_lps.cpp index a67d28c7ac..fc92015c16 100644 --- a/src/PERI/pair_peri_lps.cpp +++ b/src/PERI/pair_peri_lps.cpp @@ -409,11 +409,6 @@ double PairPeriLPS::init_one(int i, int j) { if (setflag[i][j] == 0) error->all("All pair coeffs are not set"); - cutsq[i][j] = cut[i][j] * cut[i][j]; - cutsq[j][i] = cutsq[i][j]; - - // set other j,i parameters - bulkmodulus[j][i] = bulkmodulus[i][j]; shearmodulus[j][i] = shearmodulus[i][j]; s00[j][i] = s00[i][j]; diff --git a/src/PERI/pair_peri_pmb.cpp b/src/PERI/pair_peri_pmb.cpp index ccd2a2d627..afb1f017a8 100644 --- a/src/PERI/pair_peri_pmb.cpp +++ b/src/PERI/pair_peri_pmb.cpp @@ -346,11 +346,6 @@ double PairPeriPMB::init_one(int i, int j) { if (setflag[i][j] == 0) error->all("All pair coeffs are not set"); - cutsq[i][j] = cut[i][j] * cut[i][j]; - cutsq[j][i] = cutsq[i][j]; - - // set other j,i parameters - kspring[j][i] = kspring[i][j]; alpha[j][i] = alpha[i][j]; s00[j][i] = s00[i][j]; diff --git a/src/USER-REAXC/fix_reax_c.cpp b/src/USER-REAXC/fix_reax_c.cpp index 1e440dd765..c3af1de47e 100644 --- a/src/USER-REAXC/fix_reax_c.cpp +++ b/src/USER-REAXC/fix_reax_c.cpp @@ -152,5 +152,5 @@ void FixReaxC::unpack_comm(int n, int first, double *buf) m = 0; last = first + n; for (i = first; i < last; i++) - num_bonds[i] = buf[m++]; + num_bonds[i] = static_cast (buf[m++]); } diff --git a/src/USER-REAXC/reaxc_types.h b/src/USER-REAXC/reaxc_types.h index 0be09f9dcd..be0378de39 100644 --- a/src/USER-REAXC/reaxc_types.h +++ b/src/USER-REAXC/reaxc_types.h @@ -30,7 +30,6 @@ #include "string.h" #include "sys/time.h" #include "time.h" -#include "zlib.h" /************* SOME DEFS - crucial for reax_types.h *********/ diff --git a/src/atom.cpp b/src/atom.cpp index 3d323c3c75..bfb78d6721 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -55,6 +55,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) firstgroupname = NULL; sortfreq = 1000; + nextsort = 0; userbinsize = 0.0; maxbin = maxnext = 0; binhead = NULL; @@ -1357,6 +1358,10 @@ void Atom::sort() { int i,m,n,ix,iy,iz,ibin,empty; + // set next timestep for sorting to take place + + nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq; + // re-setup sort bins if needed if (domain->box_change) setup_sort_bins(); @@ -1440,10 +1445,6 @@ void Atom::sort() //int flagall; //MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); //if (flagall) error->all("Atom sort did not operate correctly"); - - // set next timestep for sorting to take place - - nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq; } /* ---------------------------------------------------------------------- diff --git a/src/comm.cpp b/src/comm.cpp index 975799f373..4a3c942a45 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -77,7 +77,7 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp) allocate_swap(maxswap); sendlist = (int **) memory->smalloc(maxswap*sizeof(int *),"comm:sendlist"); - maxsendlist = new int[maxswap]; + memory->create(maxsendlist,maxswap,"comm:maxsendlist"); for (int i = 0; i < maxswap; i++) { maxsendlist[i] = BUFMIN; memory->create(sendlist[i],BUFMIN,"comm:sendlist[i]"); @@ -101,7 +101,7 @@ Comm::~Comm() if (sendlist) for (int i = 0; i < maxswap; i++) memory->destroy(sendlist[i]); memory->sfree(sendlist); - delete [] maxsendlist; + memory->destroy(maxsendlist); memory->destroy(buf_send); memory->destroy(buf_recv); @@ -1225,7 +1225,7 @@ void Comm::free_swap() memory->destroy(slabhi); memory->destroy(firstrecv); memory->destroy(pbc_flag); - memory->destroy_2d_int_array(pbc); + memory->destroy(pbc); } /* ---------------------------------------------------------------------- diff --git a/src/neigh_derive.cpp b/src/neigh_derive.cpp index 022e011880..c33fac9e62 100644 --- a/src/neigh_derive.cpp +++ b/src/neigh_derive.cpp @@ -161,6 +161,7 @@ void Neighbor::half_from_full_newton(NeighList *list) build skip list for subset of types from parent list iskip and ijskip flag which atom types and type pairs to skip this is for half and full lists + if ghostflag, also store neighbors of ghost atoms & set inum,gnum correctly ------------------------------------------------------------------------- */ void Neighbor::skip_from(NeighList *list) @@ -169,6 +170,7 @@ void Neighbor::skip_from(NeighList *list) int *neighptr,*jlist; int *type = atom->type; + int nlocal = atom->nlocal; int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; @@ -178,7 +180,8 @@ void Neighbor::skip_from(NeighList *list) int *ilist_skip = list->listskip->ilist; int *numneigh_skip = list->listskip->numneigh; int **firstneigh_skip = list->listskip->firstneigh; - int inum_skip = list->listskip->inum; + int num_skip = list->listskip->inum; + if (list->ghostflag) num_skip += list->listskip->gnum; int *iskip = list->iskip; int **ijskip = list->ijskip; @@ -191,7 +194,7 @@ void Neighbor::skip_from(NeighList *list) // skip I atom entirely if iskip is set for type[I] // skip I,J pair if ijskip is set for type[I],type[J] - for (ii = 0; ii < inum_skip; ii++) { + for (ii = 0; ii < num_skip; ii++) { i = ilist_skip[ii]; itype = type[i]; if (iskip[itype]) continue; @@ -226,6 +229,14 @@ void Neighbor::skip_from(NeighList *list) } list->inum = inum; + if (list->ghostflag) { + int num = 0; + for (i = 0; i < inum; i++) + if (ilist[i] < nlocal) num++; + else break; + list->inum = num; + list->gnum = inum - num; + } } /* ---------------------------------------------------------------------- @@ -486,6 +497,7 @@ void Neighbor::copy_from(NeighList *list) NeighList *listcopy = list->listcopy; list->inum = listcopy->inum; + list->gnum = listcopy->gnum; list->ilist = listcopy->ilist; list->numneigh = listcopy->numneigh; list->firstneigh = listcopy->firstneigh; diff --git a/src/neigh_full.cpp b/src/neigh_full.cpp index de5d576dfc..a173fefeb2 100644 --- a/src/neigh_full.cpp +++ b/src/neigh_full.cpp @@ -55,6 +55,8 @@ void Neighbor::full_nsq(NeighList *list) int npage = 0; int npnt = 0; + // loop over owned atoms, storing neighbors + for (i = 0; i < nlocal; i++) { if (pgsize - npnt < oneatom) { @@ -101,6 +103,109 @@ void Neighbor::full_nsq(NeighList *list) } list->inum = inum; + list->gnum = 0; +} + +/* ---------------------------------------------------------------------- + N^2 search for all neighbors + include neighbors of ghost atoms + every neighbor pair appears in list of both atoms i and j +------------------------------------------------------------------------- */ + +void Neighbor::full_nsq_ghost(NeighList *list) +{ + int i,j,n,itype,jtype,which; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr; + + int **special = atom->special; + int **nspecial = atom->nspecial; + int *tag = atom->tag; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + + int inum = 0; + int npage = 0; + int npnt = 0; + + // loop over owned & ghost atoms, storing neighbors + + for (i = 0; i < nall; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms, owned and ghost + // skip i = j + + if (i < nlocal) { + for (j = 0; j < nall; j++) { + if (i == j) continue; + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq <= cutneighsq[itype][jtype]) { + if (molecular) which = find_special(special[i],nspecial[i],tag[j]); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } + } + } else { + for (j = 0; j < nall; j++) { + if (i == j) continue; + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq <= cutneighghostsq[itype][jtype]) { + if (molecular) which = find_special(special[i],nspecial[i],tag[j]); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + + list->inum = atom->nlocal; + list->gnum = inum - atom->nlocal; } /* ---------------------------------------------------------------------- @@ -114,12 +219,10 @@ void Neighbor::full_bin(NeighList *list) double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - // bin local & ghost atoms + // bin owned & ghost atoms bin_atoms(); - // loop over each atom, storing neighbors - int **special = atom->special; int **nspecial = atom->nspecial; int *tag = atom->tag; @@ -144,6 +247,8 @@ void Neighbor::full_bin(NeighList *list) int npage = 0; int npnt = 0; + // loop over owned atoms, storing neighbors + for (i = 0; i < nlocal; i++) { if (pgsize - npnt < oneatom) { @@ -195,6 +300,135 @@ void Neighbor::full_bin(NeighList *list) } list->inum = inum; + list->gnum = 0; +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction for all neighbors + include neighbors of ghost atoms + every neighbor pair appears in list of both atoms i and j +------------------------------------------------------------------------- */ + +void Neighbor::full_bin_ghost(NeighList *list) +{ + int i,j,k,n,itype,jtype,ibin,which; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int xbin,ybin,zbin,xbin2,ybin2,zbin2; + int *neighptr; + + // bin owned & ghost atoms + + bin_atoms(); + + int **special = atom->special; + int **nspecial = atom->nspecial; + int *tag = atom->tag; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + int nstencil = list->nstencil; + int *stencil = list->stencil; + int **stencilxyz = list->stencilxyz; + + int inum = 0; + int npage = 0; + int npnt = 0; + + // loop over owned & ghost atoms, storing neighbors + + for (i = 0; i < nall; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in surrounding bins in stencil including self + // when i is a ghost atom, must check if stencil bin is out of bounds + // skip i = j + + if (i < nlocal) { + ibin = coord2bin(x[i]); + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (i == j) continue; + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq <= cutneighsq[itype][jtype]) { + if (molecular) which = find_special(special[i],nspecial[i],tag[j]); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } + } + } + + } else { + ibin = coord2bin(x[i],xbin,ybin,zbin); + for (k = 0; k < nstencil; k++) { + xbin2 = xbin + stencilxyz[k][0]; + ybin2 = ybin + stencilxyz[k][1]; + zbin2 = zbin + stencilxyz[k][2]; + if (xbin2 < 0 || xbin2 >= mbinx || + ybin2 < 0 || ybin2 >= mbiny || + zbin2 < 0 || zbin2 >= mbinz) continue; + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (i == j) continue; + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq <= cutneighghostsq[itype][jtype]) { + if (molecular) which = find_special(special[i],nspecial[i],tag[j]); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } + } + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + + list->inum = atom->nlocal; + list->gnum = inum - atom->nlocal; } /* ---------------------------------------------------------------------- @@ -297,4 +531,5 @@ void Neighbor::full_multi(NeighList *list) } list->inum = inum; + list->gnum = 0; } diff --git a/src/neigh_list.cpp b/src/neigh_list.cpp index 6eda7d474e..909c286b0b 100644 --- a/src/neigh_list.cpp +++ b/src/neigh_list.cpp @@ -30,10 +30,10 @@ enum{NSQ,BIN,MULTI}; // also in neighbor.cpp NeighList::NeighList(LAMMPS *lmp, int size) : Pointers(lmp) { - maxlocal = 0; + maxatoms = 0; pgsize = size; - inum = 0; + inum = gnum = 0; ilist = NULL; numneigh = NULL; firstneigh = NULL; @@ -59,6 +59,7 @@ NeighList::NeighList(LAMMPS *lmp, int size) : Pointers(lmp) maxstencil = 0; stencil = NULL; + stencilxyz = NULL; maxstencil_multi = 0; nstencil_multi = NULL; @@ -88,6 +89,8 @@ NeighList::~NeighList() memory->destroy(ijskip); if (maxstencil) memory->destroy(stencil); + if (ghostflag) memory->destroy(stencilxyz); + if (maxstencil_multi) { for (int i = 1; i <= atom->ntypes; i++) { memory->destroy(stencil_multi[i]); @@ -106,22 +109,27 @@ NeighList::~NeighList() void NeighList::grow(int nmax) { - // skip if grow not needed + // skip if grow not needed by this list + // each list stores own maxatoms, b/c list->grow() called at different times + // if list does not store neighbors of ghosts, compare nmax to maxatoms + // else compare nlocal+nghost to maxatoms + // if reset list size, set it to nmax - if (nmax <= maxlocal) return; - maxlocal = nmax; + if (!ghostflag && atom->nlocal <= maxatoms) return; + if (ghostflag && atom->nlocal+atom->nghost <= maxatoms) return; + maxatoms = nmax; memory->destroy(ilist); memory->destroy(numneigh); memory->sfree(firstneigh); - memory->create(ilist,maxlocal,"neighlist:ilist"); - memory->create(numneigh,maxlocal,"neighlist:numneigh"); + memory->create(ilist,maxatoms,"neighlist:ilist"); + memory->create(numneigh,maxatoms,"neighlist:numneigh"); firstneigh = (int **) - memory->smalloc(maxlocal*sizeof(int *),"neighlist:firstneigh"); + memory->smalloc(maxatoms*sizeof(int *),"neighlist:firstneigh"); if (dnum) firstdouble = (double **) - memory->smalloc(maxlocal*sizeof(double *),"neighlist:firstdouble"); + memory->smalloc(maxatoms*sizeof(double *),"neighlist:firstdouble"); } /* ---------------------------------------------------------------------- @@ -138,6 +146,10 @@ void NeighList::stencil_allocate(int smax, int style) maxstencil = smax; memory->destroy(stencil); memory->create(stencil,maxstencil,"neighlist:stencil"); + if (ghostflag) { + memory->destroy(stencilxyz); + memory->create(stencilxyz,maxstencil,3,"neighlist:stencilxyz"); + } } } else { @@ -220,6 +232,7 @@ void NeighList::print_attributes() printf(" %d = build flag\n",buildflag); printf(" %d = grow flag\n",growflag); printf(" %d = stencil flag\n",stencilflag); + printf(" %d = ghost flag\n",ghostflag); printf("\n"); printf(" %d = pair\n",rq->pair); printf(" %d = fix\n",rq->fix); @@ -237,6 +250,7 @@ void NeighList::print_attributes() printf("\n"); printf(" %d = occasional\n",rq->occasional); printf(" %d = dnum\n",rq->dnum); + printf(" %d = ghost\n",rq->ghost); printf(" %d = copy\n",rq->copy); printf(" %d = skip\n",rq->skip); printf(" %d = otherlist\n",rq->otherlist); @@ -246,24 +260,26 @@ void NeighList::print_attributes() /* ---------------------------------------------------------------------- return # of bytes of allocated memory - if growflag = 0, maxlocal & maxpage will also be 0 + if growflag = 0, maxatoms & maxpage will also be 0 if stencilflag = 0, maxstencil * maxstencil_multi will also be 0 ------------------------------------------------------------------------- */ bigint NeighList::memory_usage() { bigint bytes = 0; - bytes += memory->usage(ilist,maxlocal); - bytes += memory->usage(numneigh,maxlocal); - bytes += maxlocal * sizeof(int *); + bytes += memory->usage(ilist,maxatoms); + bytes += memory->usage(numneigh,maxatoms); + bytes += maxatoms * sizeof(int *); bytes += memory->usage(pages,maxpage,pgsize); if (dnum) { - bytes += maxlocal * sizeof(double *); + bytes += maxatoms * sizeof(double *); bytes += memory->usage(dpages,maxpage,dnum*pgsize); } if (maxstencil) bytes += memory->usage(stencil,maxstencil); + if (ghostflag) bytes += memory->usage(stencilxyz,maxstencil,3); + if (maxstencil_multi) { bytes += memory->usage(stencil_multi,atom->ntypes,maxstencil_multi); bytes += memory->usage(distsq_multi,atom->ntypes,maxstencil_multi); diff --git a/src/neigh_list.h b/src/neigh_list.h index 4f14f9c020..dfb13b8d9a 100644 --- a/src/neigh_list.h +++ b/src/neigh_list.h @@ -26,10 +26,12 @@ class NeighList : protected Pointers { int buildflag; // 1 if pair_build invoked every reneigh int growflag; // 1 if stores atom-based arrays & pages int stencilflag; // 1 if stores stencil arrays + int ghostflag; // 1 if it stores neighbors of ghosts // data structs to store neighbor pairs I,J and associated values int inum; // # of I atoms neighbors are stored for + int gnum; // # of ghost atoms neighbors are stored for 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 @@ -45,7 +47,7 @@ class NeighList : protected Pointers { // iskip,ijskip are just ptrs to corresponding request int *iskip; // iskip[i] = 1 if atoms of type I are not in list - int **ijskip; // ijskip[i][j] =1 if pairs of type I,J are not in list + int **ijskip; // ijskip[i][j] = 1 if pairs of type I,J are not in list // settings and pointers for related neighbor lists and fixes @@ -64,6 +66,7 @@ class NeighList : protected Pointers { int maxstencil; // max size of stencil int nstencil; // # of bins in stencil int *stencil; // list of bin offsets + int **stencilxyz; // bin offsets in xyz dims int maxstencil_multi; // max sizes of stencils int *nstencil_multi; // # bins in each type-based multi stencil @@ -80,7 +83,7 @@ class NeighList : protected Pointers { bigint memory_usage(); private: - int maxlocal; // size of allocated atom arrays + int maxatoms; // size of allocated atom arrays }; } diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp index bd55b6e79f..ea41d7e9bc 100644 --- a/src/neigh_request.cpp +++ b/src/neigh_request.cpp @@ -37,6 +37,7 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp) gran = granhistory = 0; respainner = respamiddle = respaouter = 0; half_from_full = 0; + ghost = 0; // default is use newton_pair setting in force @@ -50,6 +51,10 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp) dnum = 0; + // default is no neighbors of ghosts + + ghost = 0; + // default is no copy or skip copy = 0; @@ -96,6 +101,7 @@ int NeighRequest::identical(NeighRequest *other) if (newton != other->newton) same = 0; if (occasional != other->occasional) same = 0; if (dnum != other->dnum) same = 0; + if (ghost != other->ghost) same = 0; if (copy != other->copy) same = 0; if (same_skip(other) == 0) same = 0; @@ -122,6 +128,7 @@ int NeighRequest::same_kind(NeighRequest *other) if (respaouter != other->respaouter) same = 0; if (half_from_full != other->half_from_full) same = 0; if (newton != other->newton) same = 0; + if (ghost != other->ghost) same = 0; return same; } @@ -169,4 +176,5 @@ void NeighRequest::copy_request(NeighRequest *other) newton = other->newton; dnum = other->dnum; + ghost = other->ghost; } diff --git a/src/neigh_request.h b/src/neigh_request.h index 2efc042c6e..7dd7f73a67 100644 --- a/src/neigh_request.h +++ b/src/neigh_request.h @@ -31,7 +31,7 @@ class NeighRequest : protected Pointers { int compute; int command; - // kind of list requested, one flag is 1, others are 0 + // kind of list requested // set by requesting class int half; // 1 if half neigh list @@ -63,6 +63,10 @@ class NeighRequest : protected Pointers { int dnum; + // 1 if also need neighbors of ghosts + + int ghost; + // set by neighbor and pair_hybrid after all requests are made // these settings do not change kind value diff --git a/src/neigh_stencil.cpp b/src/neigh_stencil.cpp index 8c41021782..c72974a43e 100644 --- a/src/neigh_stencil.cpp +++ b/src/neigh_stencil.cpp @@ -366,6 +366,28 @@ void Neighbor::stencil_full_bin_2d(NeighList *list, /* ---------------------------------------------------------------------- */ +void Neighbor::stencil_full_ghost_bin_2d(NeighList *list, + int sx, int sy, int sz) +{ + int i,j; + int *stencil = list->stencil; + int **stencilxyz = list->stencilxyz; + int nstencil = 0; + + for (j = -sy; j <= sy; j++) + for (i = -sx; i <= sx; i++) + if (bin_distance(i,j,0) < cutneighmaxsq) { + stencilxyz[nstencil][0] = i; + stencilxyz[nstencil][1] = j; + stencilxyz[nstencil][2] = 0; + stencil[nstencil++] = j*mbinx + i; + } + + list->nstencil = nstencil; +} + +/* ---------------------------------------------------------------------- */ + void Neighbor::stencil_full_bin_3d(NeighList *list, int sx, int sy, int sz) { @@ -382,6 +404,28 @@ void Neighbor::stencil_full_bin_3d(NeighList *list, list->nstencil = nstencil; } +/* ---------------------------------------------------------------------- */ + +void Neighbor::stencil_full_ghost_bin_3d(NeighList *list, + int sx, int sy, int sz) +{ + int i,j,k; + int *stencil = list->stencil; + int **stencilxyz = list->stencilxyz; + int nstencil = 0; + + for (k = -sz; k <= sz; k++) + for (j = -sy; j <= sy; j++) + for (i = -sx; i <= sx; i++) + if (bin_distance(i,j,k) < cutneighmaxsq) { + stencilxyz[nstencil][0] = i; + stencilxyz[nstencil][1] = j; + stencilxyz[nstencil][2] = k; + stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i; + } + + list->nstencil = nstencil; +} /* ---------------------------------------------------------------------- */ diff --git a/src/neighbor.cpp b/src/neighbor.cpp index cf5dad179c..8e5c5209e2 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -73,6 +73,7 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp) build_once = 0; cutneighsq = NULL; + cutneighghostsq = NULL; cuttype = NULL; cuttypesq = NULL; fixchecklist = NULL; @@ -139,6 +140,7 @@ Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp) Neighbor::~Neighbor() { memory->destroy(cutneighsq); + memory->destroy(cutneighghostsq); delete [] cuttype; delete [] cuttypesq; delete [] fixchecklist; @@ -226,6 +228,7 @@ void Neighbor::init() n = atom->ntypes; if (cutneighsq == NULL) { memory->create(cutneighsq,n+1,n+1,"neigh:cutneighsq"); + memory->create(cutneighghostsq,n+1,n+1,"neigh:cutneighghostsq"); cuttype = new double[n+1]; cuttypesq = new double[n+1]; } @@ -248,6 +251,11 @@ void Neighbor::init() cuttypesq[i] = MAX(cuttypesq[i],cut*cut); cutneighmin = MIN(cutneighmin,cut); cutneighmax = MAX(cutneighmax,cut); + + if (force->pair && force->pair->ghostneigh) { + cut = force->pair->cutghost[i][j] + skin; + cutneighghostsq[i][j] = cut*cut; + } } } cutneighmaxsq = cutneighmax * cutneighmax; @@ -463,7 +471,7 @@ void Neighbor::init() // skip: point this list at request->otherlist, copy skip info from request // half_from_full: point this list at preceeding full list // granhistory: set preceeding list's listgranhistory to this list - // also set precedding list's ptr to FixShearHistory + // also set preceeding list's ptr to FixShearHistory // respaouter: point this list at preceeding 1/2 inner/middle lists // pair and half: if there is a full non-occasional non-skip list // change this list to half_from_full and point at the full list @@ -551,11 +559,14 @@ void Neighbor::init() else stencil_create[i] = NULL; } - // set each list's build/grow/stencil flags based on neigh request + // set each list's build/grow/stencil/ghost flags based on neigh request // buildflag = 1 if its pair_build() invoked every reneighbor // growflag = 1 if it stores atom-based arrays and pages // stencilflag = 1 if it stores stencil arrays + // ghostflag = 1 if it stores neighbors of ghosts + // anyghostlist = 1 if any non-occasional list stores neighbors of ghosts + anyghostlist = 0; for (i = 0; i < nlist; i++) { lists[i]->buildflag = 1; if (pair_build[i] == NULL) lists[i]->buildflag = 0; @@ -567,6 +578,10 @@ void Neighbor::init() lists[i]->stencilflag = 1; if (style == NSQ) lists[i]->stencilflag = 0; if (stencil_create[i] == NULL) lists[i]->stencilflag = 0; + + lists[i]->ghostflag = 0; + if (requests[i]->ghost) lists[i]->ghostflag = 1; + if (requests[i]->ghost && !requests[i]->occasional) anyghostlist = 1; } #ifdef NEIGH_LIST_DEBUG @@ -773,8 +788,8 @@ int Neighbor::request(void *requestor) determine which pair_build function each neigh list needs based on settings of neigh request copy -> copy_from function - skip -> granular function if gran with granhistory - respa function if respaouter + skip -> granular function if gran with granhistory, + respa function if respaouter, skip_from function for everything else half_from_full, half, full, gran, respaouter -> choose by newton and rq->newton and tri settings @@ -804,10 +819,8 @@ void Neighbor::choose_build(int index, NeighRequest *rq) } else if (rq->half) { if (style == NSQ) { if (rq->newton == 0) { - if (newton_pair == 0) - pb = &Neighbor::half_nsq_no_newton; - else if (newton_pair == 1) - pb = &Neighbor::half_nsq_newton; + if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton; + else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton; } else if (rq->newton == 1) { pb = &Neighbor::half_nsq_newton; } else if (rq->newton == 2) { @@ -821,9 +834,7 @@ void Neighbor::choose_build(int index, NeighRequest *rq) } else if (rq->newton == 1) { if (triclinic == 0) pb = &Neighbor::half_bin_newton; else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri; - } else if (rq->newton == 2) { - pb = &Neighbor::half_bin_no_newton; - } + } else if (rq->newton == 2) pb = &Neighbor::half_bin_no_newton; } else if (style == MULTI) { if (rq->newton == 0) { if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton; @@ -832,15 +843,20 @@ void Neighbor::choose_build(int index, NeighRequest *rq) } else if (rq->newton == 1) { if (triclinic == 0) pb = &Neighbor::half_multi_newton; else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri; - } else if (rq->newton == 2) { - pb = &Neighbor::half_multi_no_newton; - } + } else if (rq->newton == 2) pb = &Neighbor::half_multi_no_newton; } } else if (rq->full) { - if (style == NSQ) pb = &Neighbor::full_nsq; - else if (style == BIN) pb = &Neighbor::full_bin; - else if (style == MULTI) pb = &Neighbor::full_multi; + if (style == NSQ) { + if (rq->ghost == 0) pb = &Neighbor::full_nsq; + else if (rq->ghost == 1) pb = &Neighbor::full_nsq_ghost; + } else if (style == BIN) { + if (rq->ghost == 0) pb = &Neighbor::full_bin; + else if (rq->ghost == 1) pb = &Neighbor::full_bin_ghost; + } else if (style == MULTI) { + if (rq->ghost == 0) pb = &Neighbor::full_multi; + else error->all("Neighbor multi not yet enabled for ghost neighbors"); + } } else if (rq->gran) { if (style == NSQ) { @@ -865,6 +881,11 @@ void Neighbor::choose_build(int index, NeighRequest *rq) error->all("Neighbor multi not yet enabled for rRESPA"); } + // general error check + + if (rq->ghost && !rq->full) + error->all("Neighbors of ghost atoms only allowed for full neighbor lists"); + pair_build[index] = pb; } @@ -961,8 +982,14 @@ void Neighbor::choose_stencil(int index, NeighRequest *rq) } else if (rq->full) { if (style == BIN) { - if (dimension == 2) sc = &Neighbor::stencil_full_bin_2d; - else if (dimension == 3) sc = &Neighbor::stencil_full_bin_3d; + if (dimension == 2) { + if (rq->ghost) sc = &Neighbor::stencil_full_ghost_bin_2d; + else sc = &Neighbor::stencil_full_bin_2d; + } + else if (dimension == 3) { + if (rq->ghost) sc = &Neighbor::stencil_full_ghost_bin_3d; + else sc = &Neighbor::stencil_full_bin_3d; + } } else if (style == MULTI) { if (dimension == 2) sc = &Neighbor::stencil_full_multi_2d; else if (dimension == 3) sc = &Neighbor::stencil_full_multi_3d; @@ -1119,14 +1146,16 @@ void Neighbor::build() } } - // if necessary, extend atom arrays in pairwise lists - // only for lists with growflag set and which are used every reneighbor + // if any lists store neighbors of ghosts: + // invoke grow() on all in case nlocal+nghost is now too big + // else only invoke grow() if nlocal has exceeded previous list size + // only for lists with growflag set and which are perpetual - if (atom->nlocal > maxlocal) { + if (anyghostlist || atom->nlocal > maxlocal) { maxlocal = atom->nmax; for (i = 0; i < nglist; i++) lists[glist[i]]->grow(maxlocal); } - + // extend atom bin list if necessary if (style != NSQ && atom->nmax > maxbin) { @@ -1612,6 +1641,42 @@ int Neighbor::coord2bin(double *x) return (iz-mbinzlo)*mbiny*mbinx + (iy-mbinylo)*mbinx + (ix-mbinxlo); } +/* ---------------------------------------------------------------------- + same as coord2bin, but also return ix,iy,iz offsets in each dim +------------------------------------------------------------------------- */ + +int Neighbor::coord2bin(double *x, int &ix, int &iy, int &iz) +{ + if (x[0] >= bboxhi[0]) + ix = static_cast ((x[0]-bboxhi[0])*bininvx) + nbinx; + else if (x[0] >= bboxlo[0]) { + ix = static_cast ((x[0]-bboxlo[0])*bininvx); + ix = MIN(ix,nbinx-1); + } else + ix = static_cast ((x[0]-bboxlo[0])*bininvx) - 1; + + if (x[1] >= bboxhi[1]) + iy = static_cast ((x[1]-bboxhi[1])*bininvy) + nbiny; + else if (x[1] >= bboxlo[1]) { + iy = static_cast ((x[1]-bboxlo[1])*bininvy); + iy = MIN(iy,nbiny-1); + } else + iy = static_cast ((x[1]-bboxlo[1])*bininvy) - 1; + + if (x[2] >= bboxhi[2]) + iz = static_cast ((x[2]-bboxhi[2])*bininvz) + nbinz; + else if (x[2] >= bboxlo[2]) { + iz = static_cast ((x[2]-bboxlo[2])*bininvz); + iz = MIN(iz,nbinz-1); + } else + iz = static_cast ((x[2]-bboxlo[2])*bininvz) - 1; + + ix -= mbinxlo; + iy -= mbinylo; + iz -= mbinzlo; + return iz*mbiny*mbinx + iy*mbinx + ix; +} + /* ---------------------------------------------------------------------- test if atom pair i,j is excluded from neighbor list due to type, group, molecule settings from neigh_modify command diff --git a/src/neighbor.h b/src/neighbor.h index 71995d30a2..4cf5edec78 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -86,6 +86,7 @@ class Neighbor : protected Pointers { int *fixchecklist; // which fixes to check double **cutneighsq; // neighbor cutneigh sq for each type pair + double **cutneighghostsq; // neighbor cutnsq for each ghost type pair double cutneighmaxsq; // cutneighmax squared double *cuttypesq; // cuttype squared @@ -130,6 +131,9 @@ class Neighbor : protected Pointers { int special_flag[4]; // flags for 1-2, 1-3, 1-4 neighbors + int anyghostlist; // 1 if any non-occasional list + // stores neighbors of ghosts + int exclude; // 0 if no type/group exclusions, 1 if yes int nex_type; // # of entries in type exclusion list @@ -155,6 +159,7 @@ class Neighbor : protected Pointers { void bin_atoms(); // bin all atoms double bin_distance(int, int, int); // distance between binx int coord2bin(double *); // mapping atom coord to a bin + int coord2bin(double *, int &, int &, int&); // ditto int exclusion(int, int, int, int, int *, int *); // test for pair exclusion void choose_build(int, class NeighRequest *); @@ -210,7 +215,9 @@ class Neighbor : protected Pointers { void half_multi_newton_tri(class NeighList *); void full_nsq(class NeighList *); + void full_nsq_ghost(class NeighList *); void full_bin(class NeighList *); + void full_bin_ghost(class NeighList *); void full_multi(class NeighList *); void half_from_full_no_newton(class NeighList *); @@ -252,7 +259,9 @@ class Neighbor : protected Pointers { void stencil_half_multi_3d_newton_tri(class NeighList *, int, int, int); void stencil_full_bin_2d(class NeighList *, int, int, int); + void stencil_full_ghost_bin_2d(class NeighList *, int, int, int); void stencil_full_bin_3d(class NeighList *, int, int, int); + void stencil_full_ghost_bin_3d(class NeighList *, int, int, int); void stencil_full_multi_2d(class NeighList *, int, int, int); void stencil_full_multi_3d(class NeighList *, int, int, int); diff --git a/src/pair.cpp b/src/pair.cpp index d48a004997..da81910642 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -55,6 +55,8 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp) respa_enable = 0; one_coeff = 0; no_virial_compute = 0; + ghostneigh = 0; + nextra = 0; pvector = NULL; diff --git a/src/pair.h b/src/pair.h index 564987f8f5..871c0688e4 100644 --- a/src/pair.h +++ b/src/pair.h @@ -29,7 +29,7 @@ class Pair : protected Pointers { double *eatom,**vatom; // accumulated per-atom energy/virial double cutforce; // max cutoff for all atom pairs - double **cutsq; // max cutoff sq for each atom pair + double **cutsq; // cutoff sq for each atom pair int **setflag; // 0/1 = whether each i,j has been set int comm_forward; // size of forward communication (0 if none) @@ -39,6 +39,8 @@ class Pair : protected Pointers { int respa_enable; // 1 if inner/middle/outer rRESPA routines int one_coeff; // 1 if allows only one coeff * * call int no_virial_compute; // 1 if does not invoke virial_compute() + int ghostneigh; // 1 if pair style needs neighbors of ghosts + double **cutghost; // cutoff for each ghost pair int tail_flag; // pair_modify flag for LJ tail correction double etail,ptail; // energy/pressure tail corrections diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index 00a2c5d65d..3e9f03ad91 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -51,6 +51,7 @@ PairHybrid::~PairHybrid() if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); + memory->destroy_2d_double_array(cutghost); memory->destroy_2d_int_array(nmap); memory->destroy_3d_int_array(map); } @@ -157,6 +158,7 @@ void PairHybrid::allocate() setflag[i][j] = 0; cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + cutghost = memory->create_2d_double_array(n+1,n+1,"pair:cutghost"); nmap = memory->create_2d_int_array(n+1,n+1,"pair:nmap"); map = memory->create_3d_int_array(n+1,n+1,nstyles,"pair:map"); @@ -187,6 +189,7 @@ void PairHybrid::settings(int narg, char **arg) if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); + memory->destroy_2d_double_array(cutghost); memory->destroy_2d_int_array(nmap); memory->destroy_3d_int_array(map); } @@ -262,6 +265,7 @@ void PairHybrid::settings(int narg, char **arg) // single_enable = 0 if any sub-style = 0 // respa_enable = 1 if any sub-style is set // no_virial_compute = 1 if any sub-style is set + // ghostneigh = 1 if any sub-style is set for (m = 0; m < nstyles; m++) if (styles[m]->single_enable == 0) single_enable = 0; @@ -269,6 +273,8 @@ void PairHybrid::settings(int narg, char **arg) if (styles[m]->respa_enable) respa_enable = 1; for (m = 0; m < nstyles; m++) if (styles[m]->no_virial_compute) no_virial_compute = 1; + for (m = 0; m < nstyles; m++) + if (styles[m]->ghostneigh) ghostneigh = 1; } /* ---------------------------------------------------------------------- @@ -455,12 +461,14 @@ double PairHybrid::init_one(int i, int j) } // call init/mixing for all sub-styles of I,J - // set cutsq in sub-style just as pair::init_one() does + // set cutsq in sub-style just as Pair::init() does via call to init_one() + // set cutghost for I,J and J,I just as sub-style does // sum tail corrections for I,J // return max cutoff of all sub-styles assigned to I,J // if no sub-styles assigned to I,J (pair_coeff none), cutmax = 0.0 returned double cutmax = 0.0; + cutghost[i][j] = cutghost[j][i] = 0.0; if (tail_flag) etail_ij = ptail_ij = 0.0; nmap[j][i] = nmap[i][j]; @@ -470,6 +478,9 @@ double PairHybrid::init_one(int i, int j) double cut = styles[map[i][j][k]]->init_one(i,j); styles[map[i][j][k]]->cutsq[i][j] = styles[map[i][j][k]]->cutsq[j][i] = cut*cut; + if (styles[map[i][j][k]]->ghostneigh) + cutghost[i][j] = cutghost[j][i] = + MAX(cutghost[i][j],styles[map[i][j][k]]->cutghost[i][j]); if (tail_flag) { etail_ij += styles[map[i][j][k]]->etail_ij; ptail_ij += styles[map[i][j][k]]->ptail_ij; @@ -491,7 +502,7 @@ void PairHybrid::modify_requests() // loop over pair requests only // if list is skip list and not copy, look for non-skip list of same kind - // if one exists, point at that one + // if one exists, point at that one via otherlist // else make new non-skip request of same kind and point at that one // don't bother to set ID for new request, since pair hybrid ignores list // only exception is half_from_full: