--- branches/development/src/brains/Thermo.cpp	2012/08/29 18:13:11	1787
+++ branches/development/src/brains/Thermo.cpp	2013/04/29 17:53:48	1867
@@ -35,7 +35,7 @@
  *                                                                      
  * [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).          
+ * [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).
  */
@@ -324,7 +324,6 @@ namespace OpenMD {
       Molecule* mol;
       Atom* atom;
       RealType charge;
-      RealType moment(0.0);
       Vector3d ri(0.0);
       Vector3d dipoleVector(0.0);
       Vector3d nPos(0.0);
@@ -440,7 +439,6 @@ namespace OpenMD {
     RealType kinetic;
     RealType potential;
     RealType eatom;
-    RealType AvgE_a_ = 0;
     // Convective portion of the heat flux
     Vector3d heatFluxJc = V3Zero;
 
@@ -621,8 +619,8 @@ namespace OpenMD {
   }        
    
   /**
-   * Return intertia tensor for entire system and angular momentum
-   * Vector.
+   * \brief Return inertia tensor for entire system and angular momentum
+   *  Vector.
    *
    *
    *
@@ -702,7 +700,68 @@ namespace OpenMD {
     
     return;
   }
+
+
+  Mat3x3d Thermo::getBoundingBox(){
+    
+    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
+    
+    if (!(snap->hasBoundingBox)) {
+      
+      SimInfo::MoleculeIterator i;
+      Molecule::RigidBodyIterator ri;
+      Molecule::AtomIterator ai;
+      Molecule* mol;
+      RigidBody* rb;
+      Atom* atom;
+      Vector3d pos, bMax, bMin;
+      int index = 0;
+      
+      for (mol = info_->beginMolecule(i); mol != NULL; 
+           mol = info_->nextMolecule(i)) {
+        
+        //change the positions of atoms which belong to the rigidbodies
+        for (rb = mol->beginRigidBody(ri); rb != NULL; 
+             rb = mol->nextRigidBody(ri)) {          
+          rb->updateAtoms();
+        }
+        
+        for(atom = mol->beginAtom(ai); atom != NULL;
+            atom = mol->nextAtom(ai)) {
+          
+          pos = atom->getPos();
+
+          if (index == 0) {
+            bMax = pos;
+            bMin = pos;
+          } else {
+            for (int i = 0; i < 3; i++) {
+              bMax[i] = max(bMax[i], pos[i]);
+              bMin[i] = min(bMin[i], pos[i]);
+            }
+          }
+          index++;
+        }
+      }
+      
+#ifdef IS_MPI
+      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bMax[0], 3, MPI::REALTYPE, 
+                                MPI::MAX);
 
+      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bMin[0], 3, MPI::REALTYPE, 
+                                MPI::MIN);
+#endif
+      Mat3x3d bBox = Mat3x3d(0.0);
+      for (int i = 0; i < 3; i++) {           
+        bBox(i,i) = bMax[i] - bMin[i];
+      }
+      snap->setBoundingBox(bBox);
+    }
+    
+    return snap->getBoundingBox();    
+  }
+  
+  
   // Returns the angular momentum of the system
   Vector3d Thermo::getAngularMomentum(){
     Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
@@ -828,9 +887,9 @@ namespace OpenMD {
           data[0] = pos1.x();
           data[1] = pos1.y();
           data[2] = pos1.z();          
-          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
+          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc1);
         } else {
-          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
+          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc1);
           pos1 = Vector3d(data);
         }
 
@@ -839,10 +898,10 @@ namespace OpenMD {
           pos2 = sd2->getPos();
           data[0] = pos2.x();
           data[1] = pos2.y();
-          data[2] = pos2.z();          
-          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
+          data[2] = pos2.z();  
+          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc2);
         } else {
-          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
+          MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc2);
           pos2 = Vector3d(data);
         }
 #else
@@ -898,10 +957,15 @@ namespace OpenMD {
       // Compute surface Mesh
       surfaceMesh_->computeHull(localSites_);
       snap->setHullVolume(surfaceMesh_->getVolume());
+
+      delete surfaceMesh_;
     }
+
     return snap->getHullVolume();
 #else
     return 0.0;
 #endif
   }
+
+
 }