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

src/domain.cpp

@@ -678,7 +678,7 @@ void Domain::remap(double *x, int &image)
 
 /* ----------------------------------------------------------------------
    unmap the point via image flags
-   don't reset image flag
+   x overwritten with result, don't reset image flag
    for triclinic, use h[] to add in tilt factors in other dims as needed
 ------------------------------------------------------------------------- */
 
@@ -699,6 +699,29 @@ void Domain::unmap(double *x, int image)
   }
 }
 
+/* ----------------------------------------------------------------------
+   unmap the point via image flags
+   result returned in y, don't reset image flag
+   for triclinic, use h[] to add in tilt factors in other dims as needed
+------------------------------------------------------------------------- */
+
+void Domain::unmap(double *x, int image, double *y)
+{
+  int xbox = (image & 1023) - 512;
+  int ybox = (image >> 10 & 1023) - 512;
+  int zbox = (image >> 20) - 512;
+
+  if (triclinic == 0) {
+    y[0] = x[0] + xbox*xprd;
+    y[1] = x[1] + ybox*yprd;
+    y[2] = x[2] + zbox*zprd;
+  } else {
+    y[0] = x[0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
+    y[1] = x[1] + h[1]*ybox + h[3]*zbox;
+    y[2] = x[2] + h[2]*zbox;
+  }
+}
+
 /* ----------------------------------------------------------------------
    create a lattice
    delete it if style = none

src/domain.h

@@ -91,6 +91,7 @@ class Domain : protected Pointers {
   void pbc();
   void remap(double *, int &);
   void unmap(double *, int);
+  void unmap(double *, int, double *);
   void minimum_image(double &, double &, double &);
   void minimum_image(double *);
   void set_lattice(int, char **);

src/fix_msd.cpp

@@ -72,14 +72,8 @@ FixMSD::FixMSD(LAMMPS *lmp, int narg, char **arg) :
   int xbox,ybox,zbox;
 
   for (int i = 0; i < nlocal; i++) {
-    if (mask[i] & groupbit) {
+    if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
-      xbox = (image[i] & 1023) - 512;
+    else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
-      ybox = (image[i] >> 10 & 1023) - 512;
-      zbox = (image[i] >> 20) - 512;
-      xoriginal[i][0] = x[i][0] + xbox*xprd;
-      xoriginal[i][1] = x[i][1] + ybox*yprd;
-      xoriginal[i][2] = x[i][2] + zbox*zprd;
-    } else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
   }
 
   // nmsd = # of atoms in group
@@ -149,6 +143,7 @@ void FixMSD::end_of_step()
   int *image = atom->image;
   int nlocal = atom->nlocal;
 
+  double *h = domain->h;
   double xprd = domain->xprd;
   double yprd = domain->yprd;
   double zprd = domain->zprd;
@@ -157,19 +152,39 @@ void FixMSD::end_of_step()
   double msd[4];
   msd[0] = msd[1] = msd[2] = msd[3] = 0.0;
 
-  for (int i = 0; i < nlocal; i++)
+  // dx,dy,dz = displacement of atom from original position
-    if (mask[i] & groupbit) {
+  // for triclinic, need to unwrap current atom coord via h matrix
-      xbox = (image[i] & 1023) - 512;
+
-      ybox = (image[i] >> 10 & 1023) - 512;
+  if (domain->triclinic == 0) {
-      zbox = (image[i] >> 20) - 512;
+    for (int i = 0; i < nlocal; i++)
-      dx = x[i][0] + xbox*xprd - xoriginal[i][0];
+      if (mask[i] & groupbit) {
-      dy = x[i][1] + ybox*yprd - xoriginal[i][1];
+	xbox = (image[i] & 1023) - 512;
-      dz = x[i][2] + zbox*zprd - xoriginal[i][2];
+	ybox = (image[i] >> 10 & 1023) - 512;
-      msd[0] += dx*dx;
+	zbox = (image[i] >> 20) - 512;
-      msd[1] += dy*dy;
+	dx = x[i][0] + xbox*xprd - xoriginal[i][0];
-      msd[2] += dz*dz;
+	dy = x[i][1] + ybox*yprd - xoriginal[i][1];
-      msd[3] += dx*dx + dy*dy + dz*dz;
+	dz = x[i][2] + zbox*zprd - xoriginal[i][2];
-    }
+	msd[0] += dx*dx;
+	msd[1] += dy*dy;
+	msd[2] += dz*dz;
+	msd[3] += dx*dx + dy*dy + dz*dz;
+      }
+
+  } else {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+	xbox = (image[i] & 1023) - 512;
+	ybox = (image[i] >> 10 & 1023) - 512;
+	zbox = (image[i] >> 20) - 512;
+	dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox - xoriginal[i][0];
+	dy = x[i][1] + h[1]*ybox + h[3]*zbox - xoriginal[i][1];
+	dz = x[i][2] + h[2]*zbox - xoriginal[i][2];
+	msd[0] += dx*dx;
+	msd[1] += dy*dy;
+	msd[2] += dz*dz;
+	msd[3] += dx*dx + dy*dy + dz*dz;
+      }
+  }
 
   double msd_all[4];
   MPI_Allreduce(msd,msd_all,4,MPI_DOUBLE,MPI_SUM,world);