diff --git a/src/Configuration.cpp b/src/Configuration.cpp
index 0df5f36bc5b567354d634450903ee9f51c22186a..d0017a18ba17d6cc5fa3d3fe5ffad0dc9e31631d 100644
--- a/src/Configuration.cpp
+++ b/src/Configuration.cpp
@@ -348,6 +348,12 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
 	    }
 	}
 
+	std::map<std::string,float> grid_mean_dict;
+	for (const auto& pair : part_grid_dictionary)
+	{
+	    grid_mean_dict.insert({pair.first, pair.second.mean()});
+	}
+
 	// Then assign grid addresses to particles
 	for (int i = 0; i < numParts; i++)
         {
@@ -358,7 +364,7 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
 	    {
 		part[i].pmf[j] = &(part_grid_dictionary.find( std::string(partGridFile[i][j]) )->second);
 		part[i].pmf_scale[j] = partGridFileScale[i][j];
-		part[i].meanPmf[j] = part[i].pmf[j]->mean(); // TODO: review how this is used and decide whether to scale
+		part[i].meanPmf[j] = grid_mean_dict.find( std::string(partGridFile[i][j]) )->second * part[i].pmf_scale[j];
 	    }
 		if (partForceXGridFile[i].length() != 0) {
 			part[i].forceXGrid = new BaseGrid(partForceXGridFile[i].val());