diff --git a/doc/compute_msd.txt b/doc/compute_msd.txt
index 9ada35314caf83184b020ead539b0a11ff766f24..c95f5a8b6b072b5dd1577730f2b3a7a6647fd158 100644
--- a/doc/compute_msd.txt
+++ b/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.
diff --git a/src/compute_msd.cpp b/src/compute_msd.cpp
index 0d35a886202e3b3281cf814e096410e1576de1d3..7200da83f41d791769e607902bf2a44caa5e2abe 100644
--- a/src/compute_msd.cpp
+++ b/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++;
+ }
}
}
diff --git a/src/compute_msd.h b/src/compute_msd.h
index 29d25ceb5385af02087ee699b6448c281ada661c..63abf4c4eb60018ac147046d7d12096e83c20208 100644
--- a/src/compute_msd.h
+++ b/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;