git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8008 f3b2605a-c512-4ea7-a41b-209d697bcdaa

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

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

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

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
+  remap() puts atoms outside the new box back into the new box
-  perform irregular on atoms in lamda coords to get atoms to new procs
+  image_tilt() adjusts image flags due to box shape change induced by flip
-  force reneighboring on next timestep
+  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
+  // 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
+  // do not flip in x or y if non-periodic (can tilt but not flip)
-  // when yz flips and xy is non-zero, xz must also change
+  //   this is b/c the box length would be changed (dramatically) by flip
-  // this is to keep the edge vectors of the flipped shape matrix
+  // if yz tilt exceeded, adjust C vector by one B vector
-  //   an integer combination of the edge vectors of the unflipped shape matrix
+  // if xz tilt exceeded, adjust C vector by one A vector
-  // flip is performed on next timestep before reneighboring
+  // 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;
     }
   }

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

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
+  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
+  do not flip in x or y if non-periodic (can tilt but not flip)
-  when yz flips and xy is non-zero, xz must also change
+    this is b/c the box length would be changed (dramatically) by flip
-  this is to keep the edge vectors of the flipped shape matrix
+  if yz tilt exceeded, adjust C vector by one B vector
-    an integer combination of the edge vectors of the unflipped shape matrix
+  if xz tilt exceeded, adjust C vector by one A vector
-  perform irregular on atoms in lamda coords to get atoms to new procs
+  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);

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] + 

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