diff --git a/src/balance.cpp b/src/balance.cpp index ce1e71f4e9..fa2d2c2ad9 100644 --- a/src/balance.cpp +++ b/src/balance.cpp @@ -32,7 +32,7 @@ enum{NONE,UNIFORM,USER,DYNAMIC}; enum{X,Y,Z}; enum{EXPAND,CONTRACT}; -#define BALANCE_DEBUG 1 +//#define BALANCE_DEBUG 1 /* ---------------------------------------------------------------------- */ diff --git a/src/domain.cpp b/src/domain.cpp index e35f768ab7..5f52528f6d 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -1025,6 +1025,43 @@ void Domain::unmap(double *x, int image, double *y) } } +/* ---------------------------------------------------------------------- + adjust image flags due to triclinic box flip + flip operation is changing box vectors A,B,C to new A',B',C' + A' = A (A does not change) + B' = B + mA (B shifted by A) + C' = C + pB + nA (C shifted by B and/or A) + this requires the image flags change from (a,b,c) to (a',b',c') + so that x_unwrap for each atom is same before/after + x_unwrap_before = xlocal + aA + bB + cC + x_unwrap_after = xlocal + a'A' + b'B' + c'C' + this requires: + c' = c + b' = b - cp + a' = a - (b-cp)m - cn = a - b'm - cn + in other words, for xy flip, change in x flag depends on current y flag + this is b/c the xy flip dramatically changes which tiled image of + simulation box an unwrapped point maps to +------------------------------------------------------------------------- */ + +void Domain::image_flip(int m, int n, int p) +{ + int *image = atom->image; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + int xbox = (image[i] & 1023) - 512; + int ybox = (image[i] >> 10 & 1023) - 512; + int zbox = (image[i] >> 20) - 512; + + ybox -= p*zbox; + xbox -= m*ybox + n*zbox; + + image[i] = ((zbox + 512 & 1023) << 20) | + ((ybox + 512 & 1023) << 10) | (xbox + 512 & 1023); + } +} + /* ---------------------------------------------------------------------- create a lattice delete it if style = none diff --git a/src/domain.h b/src/domain.h index bf33b6ca5c..ad638bcb59 100644 --- a/src/domain.h +++ b/src/domain.h @@ -93,16 +93,17 @@ class Domain : protected Pointers { virtual void set_local_box(); virtual void reset_box(); virtual void pbc(); - void remap(double *, int &); - void remap(double *); - void remap_near(double *, double *); - void unmap(double *, int); - void unmap(double *, int, double *); int minimum_image_check(double, double, double); void minimum_image(double &, double &, double &); void minimum_image(double *); void closest_image(const double * const, const double * const, double * const); + void remap(double *, int &); + void remap(double *); + void remap_near(double *, double *); + void unmap(double *, int); + void unmap(double *, int, double *); + void image_flip(int, int, int); void set_lattice(int, char **); void add_region(int, char **); void delete_region(int, char **); diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index 95f4ad1166..3570b6d137 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -635,13 +635,10 @@ void FixDeform::init() /* ---------------------------------------------------------------------- box flipped on previous step - do not allow box flip in a non-periodic dimension (can tilt but not flip) - this is b/c the box length will be changed (dramatically) by flip - perform irregular comm to migrate atoms to new procs reset box tilts for flipped config and create new box in domain - remap to put far-away atoms back into new box - perform irregular on atoms in lamda coords to get atoms to new procs - force reneighboring on next timestep + remap() puts atoms outside the new box back into the new box + image_tilt() adjusts image flags due to box shape change induced by flip + perform irregular on atoms in lamda coords to migrate atoms to new procs ------------------------------------------------------------------------- */ void FixDeform::pre_exchange() @@ -658,6 +655,7 @@ void FixDeform::pre_exchange() int *image = atom->image; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); + domain->image_flip(flipxy,flipxz,flipyz); domain->x2lamda(atom->nlocal); irregular->migrate_atoms(); @@ -828,12 +826,13 @@ void FixDeform::end_of_step() if (varflag) modify->addstep_compute(update->ntimestep + nevery); - // if any tilt ratios exceed 0.5, set flip = 1 & compute new tilt_flip values - // do not flip in x or y if non-periodic - // when yz flips and xy is non-zero, xz must also change - // this is to keep the edge vectors of the flipped shape matrix - // an integer combination of the edge vectors of the unflipped shape matrix - // flip is performed on next timestep before reneighboring + // if any tilt ratios exceed 0.5, set flip = 1 and compute new tilt values + // do not flip in x or y if non-periodic (can tilt but not flip) + // this is b/c the box length would be changed (dramatically) by flip + // if yz tilt exceeded, adjust C vector by one B vector + // if xz tilt exceeded, adjust C vector by one A vector + // if xy tilt exceeded, adjust B vector by one A vector + // flip is performed on next timestep, before reneighboring in pre-exchange() if (triclinic) { double xprd = set[0].hi_target - set[0].lo_target; @@ -850,36 +849,40 @@ void FixDeform::end_of_step() set[4].tilt_flip = set[4].tilt_target; set[5].tilt_flip = set[5].tilt_target; + flipxy = flipxz = flipyz = 0; + if (domain->yperiodic) { if (set[3].tilt_flip*yprdinv < -0.5) { set[3].tilt_flip += yprd; set[4].tilt_flip += set[5].tilt_flip; - flip = 1; + flipyz = 1; } else if (set[3].tilt_flip*yprdinv > 0.5) { set[3].tilt_flip -= yprd; set[4].tilt_flip -= set[5].tilt_flip; - flipy = 1; + flipyz = -1; } } if (domain->xperiodic) { if (set[4].tilt_flip*xprdinv < -0.5) { set[4].tilt_flip += xprd; - flip = 1; + flipxz = 1; } if (set[4].tilt_flip*xprdinv > 0.5) { set[4].tilt_flip -= xprd; - flip = 1; + flipxz = -1; } if (set[5].tilt_flip*xprdinv < -0.5) { set[5].tilt_flip += xprd; - flip = 1; + flipxy = 1; } if (set[5].tilt_flip*xprdinv > 0.5) { set[5].tilt_flip -= xprd; - flip = 1; + flipxy = -1; } } + flip = 0; + if (flipxy || flipxz || flipyz) flip = 1; if (flip) next_reneighbor = update->ntimestep + 1; } } diff --git a/src/fix_deform.h b/src/fix_deform.h index 7e5b7acb06..33a86b375b 100644 --- a/src/fix_deform.h +++ b/src/fix_deform.h @@ -38,7 +38,7 @@ class FixDeform : public Fix { private: int triclinic,scaleflag; - int flip,flipx,flipy; + int flip,flipxy,flipxz,flipyz; double *h_rate,*h_ratelo; int varflag; // 1 if VARIABLE option is used, 0 if not int kspace_flag; // 1 if KSpace invoked, 0 if not diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index 1dff0347d6..34d80ac251 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -85,7 +85,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (domain->xz != 0.0) scalexz = 1; } - // Set fixed-point to default = center of cell + // set fixed-point to default = center of cell fixedpoint[0] = 0.5*(domain->boxlo[0]+domain->boxhi[0]); fixedpoint[1] = 0.5*(domain->boxlo[1]+domain->boxhi[1]); @@ -116,7 +116,8 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) t_stop = atof(arg[iarg+2]); t_period = atof(arg[iarg+3]); if (t_start < 0.0 || t_stop <= 0.0) - error->all(FLERR,"Target temperature for fix nvt/npt/nph cannot be 0.0"); + error->all(FLERR, + "Target temperature for fix nvt/npt/nph cannot be 0.0"); iarg += 4; } else if (strcmp(arg[iarg],"iso") == 0) { @@ -2137,12 +2138,16 @@ void FixNH::nh_omega_dot() } /* ---------------------------------------------------------------------- - if any tilt ratios exceed 0.5, set flip = 1 & compute new tilt_flip values - do not flip in x or y if non-periodic - when yz flips and xy is non-zero, xz must also change - this is to keep the edge vectors of the flipped shape matrix - an integer combination of the edge vectors of the unflipped shape matrix - perform irregular on atoms in lamda coords to get atoms to new procs + if any tilt ratios exceed 0.5, set flip = 1 and compute new tilt values + do not flip in x or y if non-periodic (can tilt but not flip) + this is b/c the box length would be changed (dramatically) by flip + if yz tilt exceeded, adjust C vector by one B vector + if xz tilt exceeded, adjust C vector by one A vector + if xy tilt exceeded, adjust B vector by one A vector + if any flip occurs, create new box in domain + remap() puts atoms outside the new box back into the new box + image_tilt() adjusts image flags due to box shape change induced by flip + perform irregular on atoms in lamda coords to migrate atoms to new procs ------------------------------------------------------------------------- */ void FixNH::pre_exchange() @@ -2156,37 +2161,41 @@ void FixNH::pre_exchange() double xtiltmax = (0.5+DELTAFLIP)*xprd; double ytiltmax = (0.5+DELTAFLIP)*yprd; - int flip = 0; + int flipxy,flipxz,flipyz; + flipxy = flipxz = flipyz = 0; if (domain->yperiodic) { if (domain->yz < -ytiltmax) { - flip = 1; domain->yz += yprd; domain->xz += domain->xy; + flipyz = 1; } else if (domain->yz >= ytiltmax) { - flip = 1; domain->yz -= yprd; domain->xz -= domain->xy; + flipyz = -1; } } if (domain->xperiodic) { if (domain->xz < -xtiltmax) { - flip = 1; domain->xz += xprd; + flipxz = 1; } else if (domain->xz >= xtiltmax) { - flip = 1; domain->xz -= xprd; + flipxz = -1; } if (domain->xy < -xtiltmax) { - flip = 1; domain->xy += xprd; + flipxy = 1; } else if (domain->xy >= xtiltmax) { - flip = 1; domain->xy -= xprd; + flipxy = -1; } } + int flip = 0; + if (flipxy || flipxz || flipyz) flip = 1; + if (flip) { domain->set_global_box(); domain->set_local_box(); @@ -2195,7 +2204,8 @@ void FixNH::pre_exchange() int *image = atom->image; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); - + domain->image_flip(flipxy,flipxz,flipyz); + domain->x2lamda(atom->nlocal); irregular->migrate_atoms(); domain->lamda2x(atom->nlocal); diff --git a/src/lattice.cpp b/src/lattice.cpp index 1eac7d2da2..3f7e7bf33b 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -313,6 +313,15 @@ Lattice::~Lattice() int Lattice::orthogonal() { + double a = orientx[0]*orienty[0] + orientx[1]*orienty[1] + + orientx[2]*orienty[2]; + double b = orienty[0]*orientz[0] + orienty[1]*orientz[1] + + orienty[2]*orientz[2]; + double c = orientx[0]*orientz[0] + orientx[1]*orientz[1] + + orientx[2]*orientz[2]; + + printf("ABC %g %g %g\n",a,b,c); + if (orientx[0]*orienty[0] + orientx[1]*orienty[1] + orientx[2]*orienty[2]) return 0; if (orienty[0]*orientz[0] + orienty[1]*orientz[1] + diff --git a/src/set.cpp b/src/set.cpp index 2a4c25731b..4b2d7f0b04 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -497,7 +497,7 @@ void Set::set(int keyword) if (yimageflag) ybox = yimage; if (zimageflag) zbox = zimage; atom->image[i] = ((zbox + 512 & 1023) << 20) | - ((ybox + 512 & 1023) << 10) | (xbox + 512 & 1023); + ((ybox + 512 & 1023) << 10) | (xbox + 512 & 1023); // set dipole moment