--- trunk/src/applications/staticProps/ObjectCount.cpp	2010/10/19 18:40:54	1513
+++ trunk/src/applications/staticProps/ObjectCount.cpp	2015/03/07 21:41:51	2071
@@ -35,8 +35,9 @@
  *                                                                      
  * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
  * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
- * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
- * [4]  Vardeman & Gezelter, in progress (2009).                        
+ * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
 
 #include <algorithm>
@@ -51,7 +52,7 @@ namespace OpenMD {
   ObjectCount::ObjectCount(SimInfo* info, const std::string& filename, 
                            const std::string& sele)
     : StaticAnalyser(info, filename), selectionScript_(sele), 
-      evaluator_(info), seleMan_(info)   {
+      seleMan_(info), evaluator_(info) {
     
     setOutputName(getPrefix(filename) + ".counts");
     
@@ -62,60 +63,55 @@ namespace OpenMD {
     }            
   }
 
-void ObjectCount::process() {
-  Molecule* mol;
-  RigidBody* rb;
-  SimInfo::MoleculeIterator mi;
-  Molecule::RigidBodyIterator rbIter;
+  void ObjectCount::process() {
+    Molecule* mol;
+    RigidBody* rb;
+    SimInfo::MoleculeIterator mi;
+    Molecule::RigidBodyIterator rbIter;
   
-  DumpReader reader(info_, dumpFilename_);    
-  int nFrames = reader.getNFrames();
-  for (int i = 0; i < nFrames; i += step_) {
-    reader.readFrame(i);
-    currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
+    counts_.clear();
+    counts_.resize(10, 0);
+    DumpReader reader(info_, dumpFilename_);    
+    int nFrames = reader.getNFrames();
+    unsigned long int nsum = 0;
+    unsigned long int n2sum = 0;
+
+    for (int i = 0; i < nFrames; i += step_) {
+      reader.readFrame(i);
+      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
     
-    for (mol = info_->beginMolecule(mi); mol != NULL; 
-         mol = info_->nextMolecule(mi)) {
-      //change the positions of atoms which belong to the rigidbodies
-      for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
-           rb = mol->nextRigidBody(rbIter)) {
-        rb->updateAtoms();        
-      }      
-    }
+      for (mol = info_->beginMolecule(mi); mol != NULL; 
+	   mol = info_->nextMolecule(mi)) {
+	//change the positions of atoms which belong to the rigidbodies
+	for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
+	     rb = mol->nextRigidBody(rbIter)) {
+	  rb->updateAtoms();        
+	}      
+      }
     
-    if (evaluator_.isDynamic()) {
+      if (evaluator_.isDynamic()) {
 	seleMan_.setSelectionSet(evaluator_.evaluate());
-    }
+      }
         
-    int count = 0;
-    int k; 
+      unsigned int count = seleMan_.getSelectionCount();
 
-    for (StuntDouble* sd = seleMan_.beginSelected(k); 
-         sd != NULL; sd = seleMan_.nextSelected(k)) {
-      count++;
+      if (counts_.size() <= count)  {
+	counts_.resize(count, 0);
+      }
+
+      counts_[count]++;
+
+      nsum += count;
+      n2sum += count * count;
     }
-    
-    if (counts_.size() < count) counts_.resize(count);
-    counts_[count]++;
-  }
    
-  int nProcessed = nFrames /step_;
+    int nProcessed = nFrames /step_;
 
-  vector<int>::iterator it;
-  unsigned long int n;
-  unsigned long int nsum = 0;
-  unsigned long int n2sum = 0;
-  for (it = counts_.begin(); it != counts_.end(); ++it) {
-    n = (*it);
-    nsum += n;
-    n2sum += n*n;
+    nAvg = nsum / nProcessed;
+    n2Avg = n2sum / nProcessed;
+    sDev = sqrt(n2Avg - nAvg*nAvg);
+    writeCounts();   
   }
-
-  nAvg = nsum / nProcessed;
-  n2Avg = n2sum / nProcessed;
-  sDev = sqrt(n2Avg - nAvg*nAvg);
-  writeCounts();   
-}
   
   void ObjectCount::writeCounts() {
     std::ofstream ofs(outputFilename_.c_str(), std::ios::binary);
@@ -125,14 +121,15 @@ void ObjectCount::process() {
       ofs << "# <N> = "<< nAvg << "\n";
       ofs << "# <N^2> = " << n2Avg << "\n";
       ofs << "# sqrt(<N^2> - <N>^2)  = " << sDev << "\n";
-      ofs << "# N\tcount(N)\n";
+      ofs << "# N\tcounts[N]\n";
+      for (unsigned int i = 0; i < counts_.size(); ++i) {
+        ofs << i << "\t" << counts_[i] << "\n";
+      }
       
-      for (int i = 0; i < counts_.size(); ++i) {
-        ofs << i <<"\t" << counts_[i]<< std::endl;
-      }        
     } else {
       
-      sprintf(painCave.errMsg, "ObjectCount: unable to open %s\n", outputFilename_.c_str());
+      sprintf(painCave.errMsg, "ObjectCount: unable to open %s\n", 
+	      outputFilename_.c_str());
       painCave.isFatal = 1;
       simError();  
     }