forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@378 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
ff12c6b9b7
commit
2fa8f3d114
26
src/comm.cpp
26
src/comm.cpp
|
|
@ -171,17 +171,19 @@ void Comm::set_procs()
|
|||
|
||||
void Comm::setup()
|
||||
{
|
||||
// for triclinic, tilted planes are neigh cutoff apart
|
||||
// but cutneigh is distance between those planes along xyz (non-tilted) axes
|
||||
// cutcomm = distance at which ghost atoms need to be acquired
|
||||
// for orthogonal, cutcomm = box coords = cutneigh in all 3 dims
|
||||
// for triclinic, cutneigh = distance between tilted planes in box coords
|
||||
// cutcomm = lamda coords = distance between those same planes
|
||||
// can be different for each dim
|
||||
|
||||
double cutneigh[3];
|
||||
double *prd,*prd_border,*sublo,*subhi;
|
||||
|
||||
if (triclinic == 0) {
|
||||
prd = prd_border = domain->prd;
|
||||
sublo = domain->sublo;
|
||||
subhi = domain->subhi;
|
||||
cutneigh[0] = cutneigh[1] = cutneigh[2] = neighbor->cutneigh;
|
||||
cutcomm[0] = cutcomm[1] = cutcomm[2] = neighbor->cutneigh;
|
||||
} else {
|
||||
prd = domain->prd;
|
||||
prd_border = domain->prd_lamda;
|
||||
|
|
@ -190,19 +192,19 @@ void Comm::setup()
|
|||
double *h_inv = domain->h_inv;
|
||||
double length;
|
||||
length = sqrt(h_inv[0]*h_inv[0] + h_inv[5]*h_inv[5] + h_inv[4]*h_inv[4]);
|
||||
cutneigh[0] = neighbor->cutneigh * length;
|
||||
cutcomm[0] = neighbor->cutneigh * length;
|
||||
length = sqrt(h_inv[1]*h_inv[1] + h_inv[3]*h_inv[3]);
|
||||
cutneigh[1] = neighbor->cutneigh * length;
|
||||
cutcomm[1] = neighbor->cutneigh * length;
|
||||
length = h_inv[2];
|
||||
cutneigh[2] = neighbor->cutneigh * length;
|
||||
cutcomm[2] = neighbor->cutneigh * length;
|
||||
}
|
||||
|
||||
// need = # of procs I need atoms from in each dim
|
||||
// for 2d, don't communicate in z
|
||||
|
||||
need[0] = static_cast<int> (cutneigh[0] * procgrid[0] / prd_border[0]) + 1;
|
||||
need[1] = static_cast<int> (cutneigh[1] * procgrid[1] / prd_border[1]) + 1;
|
||||
need[2] = static_cast<int> (cutneigh[2] * procgrid[2] / prd_border[2]) + 1;
|
||||
need[0] = static_cast<int> (cutcomm[0] * procgrid[0] / prd_border[0]) + 1;
|
||||
need[1] = static_cast<int> (cutcomm[1] * procgrid[1] / prd_border[1]) + 1;
|
||||
need[2] = static_cast<int> (cutcomm[2] * procgrid[2] / prd_border[2]) + 1;
|
||||
if (force->dimension == 2) need[2] = 0;
|
||||
|
||||
// if non-periodic, do not communicate further than procgrid-1 away
|
||||
|
|
@ -247,7 +249,7 @@ void Comm::setup()
|
|||
recvproc[iswap] = procneigh[dim][1];
|
||||
if (ineed < 2) slablo[iswap] = -BIG;
|
||||
else slablo[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
|
||||
slabhi[iswap] = sublo[dim] + cutneigh[dim];
|
||||
slabhi[iswap] = sublo[dim] + cutcomm[dim];
|
||||
if (myloc[dim] == 0) {
|
||||
if (periodicity[dim] == 0)
|
||||
slabhi[iswap] = slablo[iswap] - 1.0;
|
||||
|
|
@ -267,7 +269,7 @@ void Comm::setup()
|
|||
} else {
|
||||
sendproc[iswap] = procneigh[dim][1];
|
||||
recvproc[iswap] = procneigh[dim][0];
|
||||
slablo[iswap] = subhi[dim] - cutneigh[dim];
|
||||
slablo[iswap] = subhi[dim] - cutcomm[dim];
|
||||
if (ineed < 2) slabhi[iswap] = BIG;
|
||||
else slabhi[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
|
||||
if (myloc[dim] == procgrid[dim]-1) {
|
||||
|
|
|
|||
|
|
@ -29,6 +29,7 @@ class Comm : protected Pointers {
|
|||
int need[3]; // procs I need atoms from in each dim
|
||||
int maxforward_fix,maxreverse_fix; // comm sizes called from Fix,Pair
|
||||
int maxforward_pair,maxreverse_pair;
|
||||
double cutcomm[3]; // cutoffs used for acquiring ghost atoms
|
||||
|
||||
Comm(class LAMMPS *);
|
||||
~Comm();
|
||||
|
|
|
|||
|
|
@ -90,19 +90,19 @@ void CreateAtoms::command(int narg, char **arg)
|
|||
} else error->all("Illegal create_atoms command");
|
||||
}
|
||||
|
||||
// convert 8 corners of my sub-box from box coords to lattice coords
|
||||
// min to max = bounding box around the pts in lattice space
|
||||
// convert 8 corners of my subdomain from box coords to lattice coords
|
||||
// for orthogonal, use corner pts of my subbox
|
||||
// for triclinic, use bounding box of my subbox
|
||||
// xyz min to max = bounding box around the domain corners in lattice space
|
||||
|
||||
int triclinic = domain->triclinic;
|
||||
double *bboxlo,*bboxhi;
|
||||
double bboxlo[3],bboxhi[3];
|
||||
|
||||
if (triclinic == 0) {
|
||||
bboxlo = domain->sublo;
|
||||
bboxhi = domain->subhi;
|
||||
} else {
|
||||
bboxlo = domain->sublo_bound;
|
||||
bboxhi = domain->subhi_bound;
|
||||
}
|
||||
bboxlo[0] = domain->sublo[0]; bboxhi[0] = domain->subhi[0];
|
||||
bboxlo[1] = domain->sublo[1]; bboxhi[1] = domain->subhi[1];
|
||||
bboxlo[2] = domain->sublo[2]; bboxhi[2] = domain->subhi[2];
|
||||
} else domain->bbox(domain->sublo_lamda,domain->subhi_lamda,bboxlo,bboxhi);
|
||||
|
||||
double xmin,ymin,zmin,xmax,ymax,zmax;
|
||||
xmin = ymin = zmin = BIG;
|
||||
|
|
|
|||
|
|
@ -219,7 +219,6 @@ void Domain::set_lamda_box()
|
|||
set local subbox params
|
||||
assumes global box is defined and proc assignment has been made
|
||||
for uppermost proc, insure subhi = boxhi (in case round-off occurs)
|
||||
for triclinic, lo/hi_bound is a bbox around all 8 tilted pts of sub-box
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Domain::set_local_box()
|
||||
|
|
@ -242,26 +241,6 @@ void Domain::set_local_box()
|
|||
if (myloc[2] < procgrid[2]-1)
|
||||
subhi[2] = boxlo[2] + (myloc[2]+1) * zprd / procgrid[2];
|
||||
else subhi[2] = boxhi[2];
|
||||
|
||||
} else {
|
||||
sublo[0] = h[0]*sublo_lamda[0] + h[5]*sublo_lamda[1] +
|
||||
h[4]*sublo_lamda[2] + boxlo[0];
|
||||
sublo[1] = h[1]*sublo_lamda[1] + h[3]*sublo_lamda[2] + boxlo[1];
|
||||
sublo[2] = h[2]*sublo_lamda[2] + boxlo[2];
|
||||
|
||||
subhi[0] = sublo[0] + xprd/procgrid[0];
|
||||
subhi[1] = sublo[1] + yprd/procgrid[1];
|
||||
subhi[2] = sublo[2] + zprd/procgrid[2];
|
||||
|
||||
sublo_bound[0] = MIN(sublo[0],sublo[0] + xy/procgrid[1]);
|
||||
sublo_bound[0] = MIN(sublo_bound[0],sublo_bound[0] + xz/procgrid[2]);
|
||||
sublo_bound[1] = MIN(sublo[1],sublo[1] + yz/procgrid[2]);
|
||||
sublo_bound[2] = sublo[2];
|
||||
|
||||
subhi_bound[0] = MAX(subhi[0],subhi[0] + xy/procgrid[1]);
|
||||
subhi_bound[0] = MAX(subhi_bound[0],subhi_bound[0] + xz/procgrid[2]);
|
||||
subhi_bound[1] = MAX(subhi[1],subhi[1] + yz/procgrid[2]);
|
||||
subhi_bound[2] = subhi[2];
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -861,3 +840,64 @@ void Domain::x2lamda(double *x, double *lamda)
|
|||
lamda[1] = h_inv[1]*delta[1] + h_inv[3]*delta[2];
|
||||
lamda[2] = h_inv[2]*delta[2];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
convert 8 lamda corner pts of lo/hi box to box coords
|
||||
return bboxlo/hi = bounding box around 8 corner pts in box coords
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Domain::bbox(double *lo, double *hi, double *bboxlo, double *bboxhi)
|
||||
{
|
||||
double x[3];
|
||||
|
||||
bboxlo[0] = bboxlo[1] = bboxlo[2] = BIG;
|
||||
bboxhi[0] = bboxhi[1] = bboxhi[2] = -BIG;
|
||||
|
||||
x[0] = lo[0]; x[1] = lo[1]; x[2] = lo[2];
|
||||
lamda2x(x,x);
|
||||
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
|
||||
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]);
|
||||
|
||||
x[0] = hi[0]; x[1] = lo[1]; x[2] = lo[2];
|
||||
lamda2x(x,x);
|
||||
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
|
||||
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]);
|
||||
|
||||
x[0] = lo[0]; x[1] = hi[1]; x[2] = lo[2];
|
||||
lamda2x(x,x);
|
||||
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
|
||||
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]);
|
||||
|
||||
x[0] = hi[0]; x[1] = hi[1]; x[2] = lo[2];
|
||||
lamda2x(x,x);
|
||||
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
|
||||
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]);
|
||||
|
||||
x[0] = lo[0]; x[1] = lo[1]; x[2] = hi[2];
|
||||
lamda2x(x,x);
|
||||
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
|
||||
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]);
|
||||
|
||||
x[0] = hi[0]; x[1] = lo[1]; x[2] = hi[2];
|
||||
lamda2x(x,x);
|
||||
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
|
||||
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]);
|
||||
|
||||
x[0] = lo[0]; x[1] = hi[1]; x[2] = hi[2];
|
||||
lamda2x(x,x);
|
||||
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
|
||||
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]);
|
||||
|
||||
x[0] = hi[0]; x[1] = hi[1]; x[2] = hi[2];
|
||||
lamda2x(x,x);
|
||||
bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
|
||||
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]);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -66,10 +66,8 @@ class Domain : protected Pointers {
|
|||
double sublo[3],subhi[3]; // sub-box bounds on this proc
|
||||
|
||||
// triclinic box
|
||||
// sublo = lower-left corner with tilt
|
||||
// subhi = upper-right as if untilted
|
||||
// sublo/hi = undefined
|
||||
double sublo_lamda[3],subhi_lamda[3]; // bounds of subbox in lamda
|
||||
double sublo_bound[3],subhi_bound[3]; // bounding box of tilted subdomain
|
||||
|
||||
// triclinic box
|
||||
double xy,xz,yz; // 3 tilt factors
|
||||
|
|
@ -105,6 +103,7 @@ class Domain : protected Pointers {
|
|||
void x2lamda(int);
|
||||
void lamda2x(double *, double *);
|
||||
void x2lamda(double *, double *);
|
||||
void bbox(double *, double *, double *, double *);
|
||||
};
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -135,7 +135,7 @@ void FixRecenter::init()
|
|||
void FixRecenter::initial_integrate()
|
||||
{
|
||||
// target COM
|
||||
// bounding box around domain works for both orthog and triclinic
|
||||
// bounding box around domain works for both orthogonal and triclinic
|
||||
|
||||
double xtarget,ytarget,ztarget;
|
||||
double *bboxlo,*bboxhi;
|
||||
|
|
|
|||
133
src/neighbor.cpp
133
src/neighbor.cpp
|
|
@ -19,6 +19,7 @@
|
|||
#include "neighbor.h"
|
||||
#include "atom.h"
|
||||
#include "atom_vec.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "domain.h"
|
||||
|
|
@ -868,15 +869,16 @@ void Neighbor::build_full()
|
|||
/* ----------------------------------------------------------------------
|
||||
setup neighbor binning parameters
|
||||
bin numbering is global: 0 = 0.0 to binsize, 1 = binsize to 2*binsize
|
||||
nbin-1 = prd-binsize to binsize
|
||||
nbin = prd to prd+binsize
|
||||
-1 = -binsize to 0.0
|
||||
nbin-1 = bbox-binsize to binsize
|
||||
nbin = bbox to bbox+binsize
|
||||
-1,-2,etc = -binsize to 0.0, -2*size to -size, 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
|
||||
prd must be filled exactly by integer # of bins
|
||||
for orthogonal boxes, prd must be filled exactly by integer # of bins
|
||||
so procs on both sides of PBC see same bin boundary
|
||||
triclinic is an exception, since stencil & neigh list built differently
|
||||
for triclinic, tilted simulation box cannot contain integer # of bins
|
||||
stencil & neigh list built differently to account for this
|
||||
mbinlo = lowest global bin any of my ghost atoms could fall into
|
||||
mbinhi = highest global bin any of my ghost atoms could fall into
|
||||
mbin = number of bins I need in a dimension
|
||||
|
|
@ -886,28 +888,41 @@ void Neighbor::build_full()
|
|||
|
||||
void Neighbor::setup_bins()
|
||||
{
|
||||
// bbox = bounding box of entire domain
|
||||
// bboxlo,bboxhi = bounding box of proc's sub-domain
|
||||
// for triclinic, bounding box surrounds all 8 corner pts
|
||||
// bbox lo/hi = bounding box of entire domain
|
||||
// bbox = size of bbox of entire domain
|
||||
// for triclinic, bbox bounds all 8 corners of tilted box
|
||||
// bsubbox lo/hi = bounding box of my subdomain extended by cutneigh
|
||||
// for triclinic, subdomain with cutneigh extension is first computed
|
||||
// in lamda coords, including tilt factor via cutcomm,
|
||||
// domain->bbox then computes bbox of corner pts transformed to box coords
|
||||
|
||||
double bbox[3];
|
||||
double *bboxlo,*bboxhi;
|
||||
double bbox[3],bsubboxlo[3],bsubboxhi[3];
|
||||
|
||||
if (triclinic == 0) {
|
||||
boxlo = domain->boxlo;
|
||||
boxhi = domain->boxhi;
|
||||
bboxlo = domain->sublo;
|
||||
bboxhi = domain->subhi;
|
||||
bboxlo = domain->boxlo;
|
||||
bboxhi = domain->boxhi;
|
||||
bsubboxlo[0] = domain->sublo[0] - cutneigh;
|
||||
bsubboxlo[1] = domain->sublo[1] - cutneigh;
|
||||
bsubboxlo[2] = domain->sublo[2] - cutneigh;
|
||||
bsubboxhi[0] = domain->subhi[0] + cutneigh;
|
||||
bsubboxhi[1] = domain->subhi[1] + cutneigh;
|
||||
bsubboxhi[2] = domain->subhi[2] + cutneigh;
|
||||
} else {
|
||||
boxlo = domain->boxlo_bound;
|
||||
boxhi = domain->boxhi_bound;
|
||||
bboxlo = domain->sublo_bound;
|
||||
bboxhi = domain->subhi_bound;
|
||||
bboxlo = domain->boxlo_bound;
|
||||
bboxhi = domain->boxhi_bound;
|
||||
double lo[3],hi[3];
|
||||
lo[0] = domain->sublo_lamda[0] - comm->cutcomm[0];
|
||||
lo[1] = domain->sublo_lamda[1] - comm->cutcomm[1];
|
||||
lo[2] = domain->sublo_lamda[2] - comm->cutcomm[2];
|
||||
hi[0] = domain->subhi_lamda[0] + comm->cutcomm[0];
|
||||
hi[1] = domain->subhi_lamda[1] + comm->cutcomm[1];
|
||||
hi[2] = domain->subhi_lamda[2] + comm->cutcomm[2];
|
||||
domain->bbox(lo,hi,bsubboxlo,bsubboxhi);
|
||||
}
|
||||
|
||||
bbox[0] = boxhi[0] - boxlo[0];
|
||||
bbox[1] = boxhi[1] - boxlo[1];
|
||||
bbox[2] = boxhi[2] - boxlo[2];
|
||||
bbox[0] = bboxhi[0] - bboxlo[0];
|
||||
bbox[1] = bboxhi[1] - bboxlo[1];
|
||||
bbox[2] = bboxhi[2] - bboxlo[2];
|
||||
|
||||
// test for too many global bins in any dimension due to huge domain
|
||||
|
||||
|
|
@ -919,6 +934,7 @@ void Neighbor::setup_bins()
|
|||
|
||||
// divide box into bins
|
||||
// optimal size is roughly 1/2 the cutoff
|
||||
// use one bin even if cutoff >> bbox
|
||||
|
||||
nbinx = static_cast<int> (2.0*bbox[0]*cutneighinv);
|
||||
nbiny = static_cast<int> (2.0*bbox[1]*cutneighinv);
|
||||
|
|
@ -940,28 +956,29 @@ void Neighbor::setup_bins()
|
|||
|
||||
// mbinlo/hi = lowest and highest global bins my ghost atoms could be in
|
||||
// coord = lowest and highest values of coords for my ghost atoms
|
||||
// add in SMALL for round-off safety
|
||||
// static_cast(-1.5) = -1, so subract additional -1
|
||||
// add in SMALL for round-off safety
|
||||
|
||||
double coord;
|
||||
int mbinxhi,mbinyhi,mbinzhi;
|
||||
double coord;
|
||||
|
||||
coord = bboxlo[0] - cutneigh - SMALL*bbox[0];
|
||||
mbinxlo = static_cast<int> ((coord-boxlo[0])*bininvx);
|
||||
if (coord < 0.0) mbinxlo = mbinxlo - 1;
|
||||
coord = bboxhi[0] + cutneigh + SMALL*bbox[0];
|
||||
mbinxhi = static_cast<int> ((coord-boxlo[0])*bininvx);
|
||||
coord = bsubboxlo[0] - SMALL*bbox[0];
|
||||
mbinxlo = static_cast<int> ((coord-bboxlo[0])*bininvx);
|
||||
if (coord < bboxlo[0]) mbinxlo = mbinxlo - 1;
|
||||
coord = bsubboxhi[0] + SMALL*bbox[0];
|
||||
mbinxhi = static_cast<int> ((coord-bboxlo[0])*bininvx);
|
||||
|
||||
coord = bboxlo[1] - cutneigh - SMALL*bbox[1];
|
||||
mbinylo = static_cast<int> ((coord-boxlo[1])*bininvy);
|
||||
if (coord < 0.0) mbinylo = mbinylo - 1;
|
||||
coord = bboxhi[1] + cutneigh + SMALL*bbox[1];
|
||||
mbinyhi = static_cast<int> ((coord-boxlo[1])*bininvy);
|
||||
coord = bsubboxlo[1] - SMALL*bbox[1];
|
||||
mbinylo = static_cast<int> ((coord-bboxlo[1])*bininvy);
|
||||
if (coord < bboxlo[1]) mbinylo = mbinylo - 1;
|
||||
coord = bsubboxhi[1] + SMALL*bbox[1];
|
||||
mbinyhi = static_cast<int> ((coord-bboxlo[1])*bininvy);
|
||||
|
||||
coord = bboxlo[2] - cutneigh - SMALL*bbox[2];
|
||||
mbinzlo = static_cast<int> ((coord-boxlo[2])*bininvz);
|
||||
if (coord < 0.0) mbinzlo = mbinzlo - 1;
|
||||
coord = bboxhi[2] + cutneigh + SMALL*bbox[2];
|
||||
mbinzhi = static_cast<int> ((coord-boxlo[2])*bininvz);
|
||||
coord = bsubboxlo[2] - SMALL*bbox[2];
|
||||
mbinzlo = static_cast<int> ((coord-bboxlo[2])*bininvz);
|
||||
if (coord < bboxlo[2]) mbinzlo = mbinzlo - 1;
|
||||
coord = bsubboxhi[2] + SMALL*bbox[2];
|
||||
mbinzhi = static_cast<int> ((coord-bboxlo[2])*bininvz);
|
||||
|
||||
// extend bins by 1 to insure stencil extent is included
|
||||
|
||||
|
|
@ -1415,37 +1432,37 @@ void Neighbor::bin_atoms()
|
|||
|
||||
/* ----------------------------------------------------------------------
|
||||
convert atom coords into local bin #
|
||||
only ghost atoms will have coord >= boxhi or coord < boxlo
|
||||
take special care to insure ghosts are put in correct bins
|
||||
this is necessary so that both procs on either side of PBC
|
||||
treat a pair of atoms straddling the PBC in a consistent way
|
||||
triclinic is an exception, since stencil & neigh list built differently
|
||||
for orthogonal, only ghost atoms will have coord >= bboxhi or coord < bboxlo
|
||||
take special care to insure ghosts are put in correct bins
|
||||
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
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int Neighbor::coord2bin(double *x)
|
||||
{
|
||||
int ix,iy,iz;
|
||||
|
||||
if (x[0] >= boxhi[0])
|
||||
ix = static_cast<int> ((x[0]-boxhi[0])*bininvx) + nbinx - mbinxlo;
|
||||
else if (x[0] >= boxlo[0])
|
||||
ix = static_cast<int> ((x[0]-boxlo[0])*bininvx) - mbinxlo;
|
||||
if (x[0] >= bboxhi[0])
|
||||
ix = static_cast<int> ((x[0]-bboxhi[0])*bininvx) + nbinx - mbinxlo;
|
||||
else if (x[0] >= bboxlo[0])
|
||||
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx) - mbinxlo;
|
||||
else
|
||||
ix = static_cast<int> ((x[0]-boxlo[0])*bininvx) - mbinxlo - 1;
|
||||
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx) - mbinxlo - 1;
|
||||
|
||||
if (x[1] >= boxhi[1])
|
||||
iy = static_cast<int> ((x[1]-boxhi[1])*bininvy) + nbiny - mbinylo;
|
||||
else if (x[1] >= boxlo[1])
|
||||
iy = static_cast<int> ((x[1]-boxlo[1])*bininvy) - mbinylo;
|
||||
if (x[1] >= bboxhi[1])
|
||||
iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy) + nbiny - mbinylo;
|
||||
else if (x[1] >= bboxlo[1])
|
||||
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy) - mbinylo;
|
||||
else
|
||||
iy = static_cast<int> ((x[1]-boxlo[1])*bininvy) - mbinylo - 1;
|
||||
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy) - mbinylo - 1;
|
||||
|
||||
if (x[2] >= boxhi[2])
|
||||
iz = static_cast<int> ((x[2]-boxhi[2])*bininvz) + nbinz - mbinzlo;
|
||||
else if (x[2] >= boxlo[2])
|
||||
iz = static_cast<int> ((x[2]-boxlo[2])*bininvz) - mbinzlo;
|
||||
if (x[2] >= bboxhi[2])
|
||||
iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz) + nbinz - mbinzlo;
|
||||
else if (x[2] >= bboxlo[2])
|
||||
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz) - mbinzlo;
|
||||
else
|
||||
iz = static_cast<int> ((x[2]-boxlo[2])*bininvz) - mbinzlo - 1;
|
||||
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz) - mbinzlo - 1;
|
||||
|
||||
return (iz*mbiny*mbinx + iy*mbinx + ix + 1);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -118,7 +118,7 @@ class Neighbor : protected Pointers {
|
|||
double bininvx,bininvy,bininvz;
|
||||
|
||||
int triclinic; // 0 if domain is orthog, 1 if triclinic
|
||||
double *boxlo,*boxhi; // copies of domain bounds
|
||||
double *bboxlo,*bboxhi; // copy of full domain bounding box
|
||||
|
||||
int **pages; // half neighbor list pages
|
||||
int maxpage; // # of half pages currently allocated
|
||||
|
|
|
|||
Loading…
Reference in New Issue