From 970d6053a05f4b363c46d980c8bee25df84d48ae Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 6 Aug 2010 14:48:13 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4450
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/fix_ave_spatial.cpp | 8 +++++---
 src/force.h             | 1 +
 src/update.cpp          | 7 ++++++-
 3 files changed, 12 insertions(+), 4 deletions(-)

diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp
index 0e72fda1ca..eb05df32d4 100644
--- a/src/fix_ave_spatial.cpp
+++ b/src/fix_ave_spatial.cpp
@@ -20,6 +20,7 @@
 #include "fix_ave_spatial.h"
 #include "atom.h"
 #include "update.h"
+#include "force.h"
 #include "domain.h"
 #include "lattice.h"
 #include "modify.h"
@@ -746,9 +747,10 @@ void FixAveSpatial::end_of_step()
   // time average across samples
   // if normflag = ALL, final is total value / total count
   // if normflag = SAMPLE, final is sum of ave / repeat
-  // exception is ALL density: normalized by repeat, not total count
+  // exception is densities: normalized by repeat, not total count
 
   double repeat = nrepeat;
+  double mv2d = force->mv2d;
 
   if (normflag == ALL) {
     MPI_Allreduce(count_many,count_sum,nlayers,MPI_DOUBLE,MPI_SUM,world);
@@ -757,8 +759,8 @@ void FixAveSpatial::end_of_step()
     for (m = 0; m < nlayers; m++) {
       if (count_sum[m] > 0.0)
 	for (j = 0; j < nvalues; j++)
-	  if (which[j] == DENSITY_NUMBER || which[j] == DENSITY_MASS)
-	    values_sum[m][j] /= repeat;
+	  if (which[j] == DENSITY_NUMBER) values_sum[m][j] /= repeat;
+	  else if (which[j] == DENSITY_MASS) values_sum[m][j] *= mv2d/repeat;
 	  else values_sum[m][j] /= count_sum[m];
       count_sum[m] /= repeat;
     }
diff --git a/src/force.h b/src/force.h
index d29e8c6f3a..ef874845f6 100644
--- a/src/force.h
+++ b/src/force.h
@@ -23,6 +23,7 @@ class Force : protected Pointers {
   double boltz;                      // Boltzmann constant (eng/degree-K)
   double mvv2e;                      // conversion of mv^2 to energy
   double ftm2v;                      // conversion of ft/m to velocity
+  double mv2d;                       // conversion of mass/volume to density
   double nktv2p;                     // conversion of NkT/V to pressure
   double qqr2e;                      // conversion of q^2/r to energy
   double qe2f;                       // conversion of qE to force
diff --git a/src/update.cpp b/src/update.cpp
index 7a5f15c47b..2de63f21f8 100644
--- a/src/update.cpp
+++ b/src/update.cpp
@@ -104,6 +104,7 @@ void Update::set_units(const char *style)
     force->boltz = 1.0;
     force->mvv2e = 1.0;
     force->ftm2v = 1.0;
+    force->mv2d = 1.0;
     force->nktv2p = 1.0;
     force->qqr2e = 1.0;
     force->qe2f = 1.0;
@@ -116,6 +117,7 @@ void Update::set_units(const char *style)
     force->boltz = 0.0019872067;
     force->mvv2e = 48.88821291 * 48.88821291;
     force->ftm2v = 1.0 / 48.88821291 / 48.88821291;
+    force->mv2d = 1.0 / 0.602214179;
     force->nktv2p = 68568.415;
     force->qqr2e = 332.06371;
     force->qe2f = 23.060549; 
@@ -127,7 +129,8 @@ void Update::set_units(const char *style)
   } else if (strcmp(style,"metal") == 0) {
     force->boltz = 8.617343e-5;
     force->mvv2e = 1.0364269e-4;
-    force->ftm2v = 1 / 1.0364269e-4;
+    force->ftm2v = 1.0 / 1.0364269e-4;
+    force->mv2d = 1.0 / 0.602214179;
     force->nktv2p = 1.6021765e6;
     force->qqr2e = 14.399645;
     force->qe2f = 1.0;
@@ -140,6 +143,7 @@ void Update::set_units(const char *style)
     force->boltz = 1.3806504e-23;
     force->mvv2e = 1.0;
     force->ftm2v = 1.0;
+    force->mv2d = 1.0;
     force->nktv2p = 1.0;
     force->qqr2e = 8.9876e9;
     force->qe2f = 1.0;
@@ -152,6 +156,7 @@ void Update::set_units(const char *style)
     force->boltz = 1.3806504e-16;
     force->mvv2e = 1.0;
     force->ftm2v = 1.0;
+    force->mv2d = 1.0;
     force->nktv2p = 1.0;
     force->qqr2e = 1.0;
     force->qe2f = 1.0;