Added average keyword

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

doc/compute_msd.txt

@@ -15,14 +15,15 @@ compute ID group-ID msd keyword values ... :pre
 ID, group-ID are documented in "compute"_compute.html command :ulb,l
 msd = style name of this compute command :l
 zero or more keyword/value pairs may be appended :l
-keyword = {com} :l
-  {com} value = {yes} or {no} :pre
+keyword = {com} or {average} :l
+  {com} value = {yes} or {no}
+  {average} value = {yes} or {no} :pre
 :ule
 
 [Examples:]
 
 compute 1 all msd
-compute 1 upper msd com yes :pre
+compute 1 upper msd com yes average yes :pre
 
 [Description:]
 
@@ -41,13 +42,19 @@ averaged over atoms in the group.
 The slope of the mean-squared displacement (MSD) versus time is
 proportional to the diffusion coefficient of the diffusing atoms.
 
-The displacement of an atom is from its original position at the time
-the compute command was issued.  The value of the displacement will be
+The displacement of an atom is from its reference position. This is
+normally the original position at the time
+the compute command was issued, unless the {average} keyword is set to {yes}.  
+The value of the displacement will be
 0.0 for atoms not in the specified compute group.
 
 If the {com} option is set to {yes} then the effect of any drift in
-the center-of-mass of the group of atoms is subtracted out before xhe
-displacment of each atom is calcluated.
+the center-of-mass of the group of atoms is subtracted out before the
+displacment of each atom is calculated.
+
+If the {average} option is set to {yes} then the reference position of 
+an atom is a running average of the atom position over all the calls
+to the msd calculation since the compute was invoked.
 
 IMPORTANT NOTE: Initial coordinates are stored in "unwrapped" form, by
 using the image flags associated with each atom.  See the "dump
@@ -62,7 +69,13 @@ to be continuous when running from a "restart file"_read_restart.html,
 then you should use the same ID for this compute, as in the original
 run.  This is so that the fix this compute creates to store per-atom
 quantities will also have the same ID, and thus be initialized
-correctly with time=0 atom coordinates from the restart file.
+correctly with atom reference positions from the restart file.
+When {average} is set to yes, then the atom reference positions 
+are restored correctly, but not the number of samples used
+obtain them. As a result, the reference positions from the restart
+file are combined with subsequent positions as if they were from a 
+single sample, instead of many, which will change the values of msd
+somewhat.
 
 [Output info:]
 
@@ -86,4 +99,4 @@ msd/chunk"_compute_msd_chunk.html
 
 [Default:]
 
-The option default is com = no.
+The option default are com = no, average = no.

src/compute_msd.cpp

@@ -38,6 +38,7 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) :
   // optional args
 
   comflag = 0;
+  avflag = 0;
 
   int iarg = 3;
   while (iarg < narg) {
@@ -47,10 +48,16 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) :
       else if (strcmp(arg[iarg+1],"yes") == 0) comflag = 1;
       else error->all(FLERR,"Illegal compute msd command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"average") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal compute msd command");
+      if (strcmp(arg[iarg+1],"no") == 0) avflag = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) avflag = 1;
+      else error->all(FLERR,"Illegal compute msd command");
+      iarg += 2;
     } else error->all(FLERR,"Illegal compute msd command");
   }
 
-  // create a new fix STORE style
+  // create a new fix STORE style for reference positions
   // id = compute-ID + COMPUTE_STORE, fix group = compute group
 
   int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
@@ -97,6 +104,11 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) :
           xoriginal[i][2] -= cm[2];
         }
     }
+
+    // initialize counter for average positions if requested
+
+    if (avflag) naverage = 1;
+
   }
 
   // displacement vector
@@ -120,7 +132,7 @@ ComputeMSD::~ComputeMSD()
 
 void ComputeMSD::init()
 {
-  // set fix which stores original atom coords
+  // set fix which stores reference atom coords
 
   int ifix = modify->find_fix(id_fix);
   if (ifix < 0) error->all(FLERR,"Could not find compute msd fix ID");
@@ -144,8 +156,8 @@ void ComputeMSD::compute_vector()
   if (comflag) group->xcm(igroup,masstotal,cm);
   else cm[0] = cm[1] = cm[2] = 0.0;
 
-  // dx,dy,dz = displacement of atom from original position
-  // original unwrapped position is stored by fix
+  // dx,dy,dz = displacement of atom from reference position
+  // reference unwrapped position is stored by fix
   // relative to center of mass if comflag is set
   // for triclinic, need to unwrap current atom coord via h matrix
 
@@ -167,35 +179,63 @@ void ComputeMSD::compute_vector()
   double msd[4];
   msd[0] = msd[1] = msd[2] = msd[3] = 0.0;
 
+
+
+  double xtmp, ytmp, ztmp;
+
   if (domain->triclinic == 0) {
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
         xbox = (image[i] & IMGMASK) - IMGMAX;
         ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
         zbox = (image[i] >> IMG2BITS) - IMGMAX;
-        dx = x[i][0] + xbox*xprd - cm[0] - xoriginal[i][0];
-        dy = x[i][1] + ybox*yprd - cm[1] - xoriginal[i][1];
-        dz = x[i][2] + zbox*zprd - cm[2] - xoriginal[i][2];
+	xtmp = x[i][0] + xbox*xprd - cm[0];
+	ytmp = x[i][1] + ybox*yprd - cm[1];
+	ztmp = x[i][2] + zbox*zprd - cm[2];
+        dx = xtmp - xoriginal[i][0];
+        dy = ytmp - xoriginal[i][1];
+        dz = ztmp - xoriginal[i][2];
         msd[0] += dx*dx;
         msd[1] += dy*dy;
         msd[2] += dz*dz;
         msd[3] += dx*dx + dy*dy + dz*dz;
-      }
 
+	// use running average position for reference if requested
+
+	if (avflag) {
+	  double navfac = 1.0/(naverage+1);
+	  xoriginal[i][0] = (xoriginal[i][0]*naverage + xtmp)*navfac;
+	  xoriginal[i][1] = (xoriginal[i][1]*naverage + ytmp)*navfac;
+	  xoriginal[i][2] = (xoriginal[i][2]*naverage + ztmp)*navfac;
+	  naverage++;
+	}
+      }
   } else {
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
         xbox = (image[i] & IMGMASK) - IMGMAX;
         ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
         zbox = (image[i] >> IMG2BITS) - IMGMAX;
-        dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox -
-          cm[0] - xoriginal[i][0];
-        dy = x[i][1] + h[1]*ybox + h[3]*zbox - cm[1] - xoriginal[i][1];
-        dz = x[i][2] + h[2]*zbox - cm[2] - xoriginal[i][2];
+        xtmp = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox - cm[0];
+        ytmp = x[i][1] + h[1]*ybox + h[3]*zbox - cm[1];
+        ztmp = x[i][2] + h[2]*zbox - cm[2];
+        dx = xtmp - xoriginal[i][0];
+        dy = ytmp - xoriginal[i][1];
+        dz = ztmp - xoriginal[i][2];
         msd[0] += dx*dx;
         msd[1] += dy*dy;
         msd[2] += dz*dz;
         msd[3] += dx*dx + dy*dy + dz*dz;
+
+	// use running average position for reference if requested
+
+	if (avflag) {
+	  double navfac = 1.0/(naverage+1);
+	  xoriginal[i][0] = (xoriginal[i][0]*naverage + xtmp)*navfac;
+	  xoriginal[i][1] = (xoriginal[i][0]*naverage + xtmp)*navfac;
+	  xoriginal[i][2] = (xoriginal[i][0]*naverage + xtmp)*navfac;
+	  naverage++;
+	}
       }
   }
 

src/compute_msd.h

@@ -33,7 +33,9 @@ class ComputeMSD : public Compute {
   void set_arrays(int);
 
  protected:
-  int comflag;
+  int comflag;   // comflag = 1 if reference moves with center of mass
+  int avflag;    // avflag = 1 if using average position as reference
+  int naverage;  // number of samples for average position
   bigint nmsd;
   double masstotal;
   char *id_fix;