From 39e30f3196e53ac02bbf1d4bfa4df3aa2f2e71cb Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 15 Feb 2011 18:10:20 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5637 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/domain.cpp | 27 +++++++++++ src/domain.h | 2 + src/neighbor.cpp | 113 +++++++++++++++++++++++++++++++++-------------- src/neighbor.h | 12 ++--- 4 files changed, 116 insertions(+), 38 deletions(-) diff --git a/src/domain.cpp b/src/domain.cpp index 1bd781bc46..d92e9c1dc3 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -1138,3 +1138,30 @@ void Domain::bbox(double *lo, double *hi, double *bboxlo, double *bboxhi) bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]); bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]); } + +/* ---------------------------------------------------------------------- + compute 8 corner pts of triclinic box + 8 are ordered with x changing fastest, then y, finally z + could be more efficient if just coded with xy,yz,xz explicitly +------------------------------------------------------------------------- */ + +void Domain::box_corners() +{ + corners[0][0] = 0.0; corners[0][1] = 0.0; corners[0][2] = 0.0; + lamda2x(corners[0],corners[0]); + corners[1][0] = 1.0; corners[1][1] = 0.0; corners[1][2] = 0.0; + lamda2x(corners[1],corners[1]); + corners[2][0] = 0.0; corners[2][1] = 1.0; corners[2][2] = 0.0; + lamda2x(corners[2],corners[2]); + corners[3][0] = 1.0; corners[3][1] = 1.0; corners[3][2] = 0.0; + lamda2x(corners[3],corners[3]); + corners[4][0] = 0.0; corners[4][1] = 0.0; corners[4][2] = 1.0; + lamda2x(corners[4],corners[4]); + corners[5][0] = 1.0; corners[5][1] = 0.0; corners[5][2] = 1.0; + lamda2x(corners[5],corners[5]); + corners[6][0] = 0.0; corners[6][1] = 1.0; corners[6][2] = 1.0; + lamda2x(corners[6],corners[6]); + corners[7][0] = 1.0; corners[7][1] = 1.0; corners[7][2] = 1.0; + lamda2x(corners[7],corners[7]); +} + diff --git a/src/domain.h b/src/domain.h index dd2c004a8a..a630416708 100644 --- a/src/domain.h +++ b/src/domain.h @@ -54,6 +54,7 @@ class Domain : protected Pointers { // boxlo/hi = same as if untilted double boxlo_lamda[3],boxhi_lamda[3]; // lamda box = (0,1) double boxlo_bound[3],boxhi_bound[3]; // bounding box of tilted domain + double corners[8][3]; // 8 corner points // orthogonal box & triclinic box double minxlo,minxhi; // minimum size of global box @@ -111,6 +112,7 @@ class Domain : protected Pointers { void lamda2x(double *, double *); void x2lamda(double *, double *); void bbox(double *, double *, double *, double *); + void box_corners(); }; } diff --git a/src/neighbor.cpp b/src/neighbor.cpp index f1df926a23..c8b25246d7 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -208,10 +208,10 @@ void Neighbor::init() // cutneigh = force cutoff + skin if cutforce > 0, else cutneigh = 0 triggersq = 0.25*skin*skin; - shrinkcheck = 0; + boxcheck = 0; if (domain->box_change && (domain->xperiodic || domain->yperiodic || (dimension == 3 && domain->zperiodic))) - shrinkcheck = 1; + boxcheck = 1; n = atom->ntypes; if (cutneighsq == NULL) { @@ -1002,20 +1002,46 @@ int Neighbor::decide() /* ---------------------------------------------------------------------- if any atom moved trigger distance (half of neighbor skin) return 1 - shrink trigger distance if periodic box dimension has decreased + shrink trigger distance if box size has changed + conservative shrink procedure: + compute distance each of 8 corners of box has moved since last reneighbor + reduce skin distance by sum of 2 largest of the 8 values + new trigger = 1/2 of reduced skin distance + for orthogonal box, only need 2 lo/hi corners + for triclinic, need all 8 corners since deformations can displace all 8 ------------------------------------------------------------------------- */ int Neighbor::check_distance() { - double delx,dely,delz,rsq,deltasq; + double delx,dely,delz,rsq; + double delta,deltasq,delta1,delta2; - if (shrinkcheck) { - double delta = 0.0; - if (domain->xperiodic) delta = MIN(delta,domain->xprd-xprdhold); - if (domain->yperiodic) delta = MIN(delta,domain->yprd-yprdhold); - if (domain->zperiodic) delta = MIN(delta,domain->zprd-zprdhold); - delta = MAX(0.5*skin+delta,0.0); - deltasq = delta*delta; + if (boxcheck) { + if (triclinic == 0) { + delx = bboxlo[0] - boxlo_hold[0]; + dely = bboxlo[1] - boxlo_hold[1]; + delz = bboxlo[2] - boxlo_hold[2]; + delta1 = sqrt(delx*delx + dely*dely + delz*delz); + delx = bboxhi[0] - boxhi_hold[0]; + dely = bboxhi[1] - boxhi_hold[1]; + delz = bboxhi[2] - boxhi_hold[2]; + delta2 = sqrt(delx*delx + dely*dely + delz*delz); + delta = 0.5 * (skin - (delta1+delta2)); + deltasq = delta*delta; + } else { + domain->box_corners(); + delta1 = delta2 = 0.0; + for (int i = 0; i < 8; i++) { + delx = corners[i][0] - corners_hold[i][0]; + dely = corners[i][1] - corners_hold[i][1]; + delz = corners[i][2] - corners_hold[i][2]; + delta = sqrt(delx*delx + dely*dely + delz*delz); + if (delta > delta1) delta1 = delta; + else if (delta > delta2) delta2 = delta; + } + delta = 0.5 * (skin - (delta1+delta2)); + deltasq = delta*delta; + } } else deltasq = triggersq; double **x = atom->x; @@ -1065,9 +1091,23 @@ void Neighbor::build() xhold[i][1] = x[i][1]; xhold[i][2] = x[i][2]; } - xprdhold = domain->xprd; - yprdhold = domain->yprd; - zprdhold = domain->zprd; + if (boxcheck) { + if (triclinic == 0) { + boxlo_hold[0] = bboxlo[0]; + boxlo_hold[1] = bboxlo[1]; + boxlo_hold[2] = bboxlo[2]; + boxhi_hold[0] = bboxhi[0]; + boxhi_hold[1] = bboxhi[1]; + boxhi_hold[2] = bboxhi[2]; + } else { + domain->box_corners(); + for (i = 0; i < 8; i++) { + corners_hold[i][0] = corners[i][0]; + corners_hold[i][1] = corners[i][1]; + corners_hold[i][2] = corners[i][2]; + } + } + } } // if necessary, extend atom arrays in pairwise lists @@ -1139,8 +1179,8 @@ void Neighbor::build_one(int i) setup neighbor binning parameters bin numbering in each dimension is global: 0 = 0.0 to binsize, 1 = binsize to 2*binsize, etc - nbin-1,nbin,etc = bbox-binsize to binsize, bbox to bbox+binsize, etc - -1,-2,etc = -binsize to 0.0, -2*size to -size, etc + nbin-1,nbin,etc = bbox-binsize to bbox, bbox to bbox+binsize, etc + -1,-2,etc = -binsize to 0.0, -2*binsize to -binsize, etc code will work for any binsize since next(xyz) and stencil extend as far as necessary binsize = 1/2 of cutoff is roughly optimal @@ -1190,6 +1230,7 @@ void Neighbor::setup_bins() hi[1] = domain->subhi_lamda[1] + cutghost[1]; hi[2] = domain->subhi_lamda[2] + cutghost[2]; domain->bbox(lo,hi,bsubboxlo,bsubboxhi); + corners = domain->corners; } bbox[0] = bboxhi[0] - bboxlo[0]; @@ -1526,7 +1567,10 @@ void Neighbor::bin_atoms() /* ---------------------------------------------------------------------- convert atom coords into local bin # for orthogonal, only ghost atoms will have coord >= bboxhi or coord < bboxlo - take special care to insure ghosts are put in correct bins + take special care to insure ghosts are in correct bins even w/ roundoff + hi ghost atoms = nbin,nbin+1,etc + owned atoms = 0 to nbin-1 + lo ghost atoms = -1,-2,etc this is necessary so that both procs on either side of PBC treat a pair of atoms straddling the PBC in a consistent way for triclinic, doesn't matter since stencil & neigh list built differently @@ -1537,27 +1581,30 @@ int Neighbor::coord2bin(double *x) int ix,iy,iz; if (x[0] >= bboxhi[0]) - ix = static_cast ((x[0]-bboxhi[0])*bininvx) + nbinx - mbinxlo; - else if (x[0] >= bboxlo[0]) - ix = static_cast ((x[0]-bboxlo[0])*bininvx) - mbinxlo; - else - ix = static_cast ((x[0]-bboxlo[0])*bininvx) - mbinxlo - 1; + 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 - mbinylo; - else if (x[1] >= bboxlo[1]) - iy = static_cast ((x[1]-bboxlo[1])*bininvy) - mbinylo; - else - iy = static_cast ((x[1]-bboxlo[1])*bininvy) - mbinylo - 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 - mbinzlo; - else if (x[2] >= bboxlo[2]) - iz = static_cast ((x[2]-bboxlo[2])*bininvz) - mbinzlo; - else - iz = static_cast ((x[2]-bboxlo[2])*bininvz) - mbinzlo - 1; + 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; - return (iz*mbiny*mbinx + iy*mbinx + ix); + return (iz-mbinzlo)*mbiny*mbinx + (iy-mbinylo)*mbinx + (ix-mbinxlo); } /* ---------------------------------------------------------------------- diff --git a/src/neighbor.h b/src/neighbor.h index 3185e6563c..25df310bce 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -91,10 +91,11 @@ class Neighbor : protected Pointers { double triggersq; // trigger = build when atom moves this dist - double **xhold; // atom coords at last neighbor build - int maxhold; // size of xhold array - double xprdhold,yprdhold,zprdhold; // box size at last neighbor build - int shrinkcheck; + double **xhold; // atom coords at last neighbor build + int maxhold; // size of xhold array + int boxcheck; // 1 if need to store box size + double boxlo_hold[3],boxhi_hold[3]; // box size at last neighbor build + double corners_hold[8][3]; // box corners at last neighbor build int nbinx,nbiny,nbinz; // # of global bins int *bins; // ptr to next atom in each bin @@ -119,7 +120,8 @@ class Neighbor : protected Pointers { int triclinic; // 0 if domain is orthog, 1 if triclinic int newton_pair; // 0 if newton off, 1 if on for pairwise - double *bboxlo,*bboxhi; // copy of full domain bounding box + double *bboxlo,*bboxhi; // ptrs to full domain bounding box + double (*corners)[3]; // ptr to 8 corners of triclinic box double inner[2],middle[2]; // rRESPA cutoffs for extra lists double cut_inner_sq; // outer cutoff for inner neighbor list