From efea2a18357d0f1d43a16c0be0f65c8a8371bce2 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Sat, 8 Aug 2009 18:45:22 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3021
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/domain.cpp     | 46 ++++++++++++++++++++++++++++++++++++++
 src/domain.h       |  1 +
 src/fix_indent.cpp | 55 +++++++++++++++++++++++++++-------------------
 3 files changed, 80 insertions(+), 22 deletions(-)

diff --git a/src/domain.cpp b/src/domain.cpp
index c2b700ea2e..2ffa88d88b 100644
--- a/src/domain.cpp
+++ b/src/domain.cpp
@@ -676,6 +676,52 @@ void Domain::remap(double *x, int &image)
   if (triclinic) lamda2x(coord,x);
 }
 
+/* ----------------------------------------------------------------------
+   remap the point into the periodic box no matter how far away
+   resulting coord must satisfy lo <= coord < hi
+   MAX is important since coord - prd < lo can happen when coord = hi
+   for triclinic, point is converted to lamda coords (0-1) before doing remap
+------------------------------------------------------------------------- */
+
+void Domain::remap(double *x)
+{
+  double *lo,*hi,*period,*coord;
+  double lamda[3];
+
+  if (triclinic == 0) {
+    lo = boxlo;
+    hi = boxhi;
+    period = prd;
+    coord = x;
+  } else {
+    lo = boxlo_lamda;
+    hi = boxhi_lamda;
+    period = prd_lamda;
+    x2lamda(x,lamda);
+    coord = lamda;
+  }
+
+  if (xperiodic) {
+    while (coord[0] < lo[0]) coord[0] += period[0];
+    while (coord[0] >= hi[0]) coord[0] -= period[0];
+    coord[0] = MAX(coord[0],lo[0]);
+  }
+
+  if (yperiodic) {
+    while (coord[1] < lo[1]) coord[1] += period[1];
+    while (coord[1] >= hi[1]) coord[1] -= period[1];
+    coord[1] = MAX(coord[1],lo[1]);
+  }
+
+  if (zperiodic) {
+    while (coord[2] < lo[2]) coord[2] += period[2];
+    while (coord[2] >= hi[2]) coord[2] -= period[2];
+    coord[2] = MAX(coord[2],lo[2]);
+  }
+
+  if (triclinic) lamda2x(coord,x);
+}
+
 /* ----------------------------------------------------------------------
    unmap the point via image flags
    x overwritten with result, don't reset image flag
diff --git a/src/domain.h b/src/domain.h
index 318a01fe06..5732670d66 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -90,6 +90,7 @@ class Domain : protected Pointers {
   void reset_box();
   void pbc();
   void remap(double *, int &);
+  void remap(double *);
   void unmap(double *, int);
   void unmap(double *, int, double *);
   void minimum_image(double &, double &, double &);
diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp
index 7e4eeee9c7..ffa2adc533 100644
--- a/src/fix_indent.cpp
+++ b/src/fix_indent.cpp
@@ -171,12 +171,15 @@ void FixIndent::post_force(int vflag)
 
   if (istyle == SPHERE) {
 
-    // x1,y1,z1 = current position of indenter from original x0,y0,z0
+    // ctr = current indenter center from original x0,y0,z0
+    // remap into periodic box
 
+    double ctr[3];
     double delta = (update->ntimestep - update->beginstep) * update->dt;
-    double x1 = x0 + delta*vx;
-    double y1 = y0 + delta*vy;
-    double z1 = z0 + delta*vz;
+    ctr[0] = x0 + delta*vx;
+    ctr[1] = y0 + delta*vy;
+    ctr[2] = z0 + delta*vz;
+    domain->remap(ctr);
 
     double **x = atom->x;
     double **f = atom->f;
@@ -187,9 +190,10 @@ void FixIndent::post_force(int vflag)
 
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
-	delx = x[i][0] - x1;
-	dely = x[i][1] - y1;
-	delz = x[i][2] - z1;
+	delx = x[i][0] - ctr[0];
+	dely = x[i][1] - ctr[1];
+	delz = x[i][2] - ctr[2];
+	domain->minimum_image(delx,dely,delz);
 	r = sqrt(delx*delx + dely*dely + delz*delz);
 	if (side == OUTSIDE) {
 	  dr = r - r0;
@@ -215,21 +219,27 @@ void FixIndent::post_force(int vflag)
 
   } else if (istyle == CYLINDER) {
 
-    // c1new,c2new = current coords of indenter axis from original c1,c2
+    // ctr = current indenter axis from original c1,c2
+    // remap into periodic box
+    // 3rd coord is just near box for remap(), since isn't used
 	      
+    double ctr[3];
     double delta = (update->ntimestep - update->beginstep) * update->dt;
-    double c1new,c2new;
     if (cdim == 0) {
-      c1new = c1 + delta*vy;
-      c2new = c2 + delta*vz;
+      ctr[0] = domain->boxlo[0];
+      ctr[1] = c1 + delta*vy;
+      ctr[2] = c2 + delta*vz;
     } else if (cdim == 1) {
-      c1new = c1 + delta*vx;
-      c2new = c2 + delta*vz;
+      ctr[0] = c1 + delta*vx;
+      ctr[1] = domain->boxlo[1];
+      ctr[2] = c2 + delta*vz;
     } else {
-      c1new = c1 + delta*vx;
-      c2new = c2 + delta*vy;
+      ctr[0] = c1 + delta*vx;
+      ctr[1] = c2 + delta*vy;
+      ctr[2] = domain->boxlo[2];
     }
-    
+    domain->remap(ctr);
+
     double **x = atom->x;
     double **f = atom->f;
     int *mask = atom->mask;
@@ -241,17 +251,18 @@ void FixIndent::post_force(int vflag)
       if (mask[i] & groupbit) {
 	if (cdim == 0) {
 	  delx = 0;
-	  dely = x[i][1] - c1new;
-	  delz = x[i][2] - c2new;
+	  dely = x[i][1] - ctr[1];
+	  delz = x[i][2] - ctr[2];
 	} else if (cdim == 1) {
-	  delx = x[i][0] - c1new;
+	  delx = x[i][0] - ctr[0];
 	  dely = 0;
-	  delz = x[i][2] - c2new;
+	  delz = x[i][2] - ctr[2];
 	} else {
-	  delx = x[i][0] - c1new;
-	  dely = x[i][1] - c2new;
+	  delx = x[i][0] - ctr[0];
+	  dely = x[i][1] - ctr[1];
 	  delz = 0;
 	}
+	domain->minimum_image(delx,dely,delz);
 	r = sqrt(delx*delx + dely*dely + delz*delz);
 	if (side == OUTSIDE) {
 	  dr = r - r0;