--- trunk/src/brains/Thermo.cpp	2005/01/12 22:41:40	246
+++ trunk/src/brains/Thermo.cpp	2006/04/25 02:09:01	945
@@ -1,4 +1,4 @@
- /*
+/*
  * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
  *
  * The University of Notre Dame grants you ("Licensee") a
@@ -53,7 +53,7 @@ namespace oopse {
 
 namespace oopse {
 
-double Thermo::getKinetic() {
+  double Thermo::getKinetic() {
     SimInfo::MoleculeIterator miter;
     std::vector<StuntDouble*>::iterator iiter;
     Molecule* mol;
@@ -68,30 +68,30 @@ double Thermo::getKinetic() {
     double kinetic_global = 0.0;
     
     for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
-        for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL; 
-               integrableObject = mol->nextIntegrableObject(iiter)) {
+      for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL; 
+	   integrableObject = mol->nextIntegrableObject(iiter)) {
+        
+	double mass = integrableObject->getMass();
+	Vector3d vel = integrableObject->getVel();
+        
+	kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
+        
+	if (integrableObject->isDirectional()) {
+	  angMom = integrableObject->getJ();
+	  I = integrableObject->getI();
 
-            double mass = integrableObject->getMass();
-            Vector3d vel = integrableObject->getVel();
-
-            kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
-
-            if (integrableObject->isDirectional()) {
-                angMom = integrableObject->getJ();
-                I = integrableObject->getI();
-
-                if (integrableObject->isLinear()) {
-                    i = integrableObject->linearAxis();
-                    j = (i + 1) % 3;
-                    k = (i + 2) % 3;
-                    kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
-                } else {                        
-                    kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1) 
-                                    + angMom[2]*angMom[2]/I(2, 2);
-                }
-            }
+	  if (integrableObject->isLinear()) {
+	    i = integrableObject->linearAxis();
+	    j = (i + 1) % 3;
+	    k = (i + 2) % 3;
+	    kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
+	  } else {                        
+	    kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1) 
+	      + angMom[2]*angMom[2]/I(2, 2);
+	  }
+	}
             
-        }
+      }
     }
     
 #ifdef IS_MPI
@@ -105,49 +105,49 @@ double Thermo::getKinetic() {
     kinetic = kinetic * 0.5 / OOPSEConstant::energyConvert;
 
     return kinetic;
-}
+  }
 
-double Thermo::getPotential() {
+  double Thermo::getPotential() {
     double potential = 0.0;
     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
-    double potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] + 
-                                             curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
+    double shortRangePot_local =  curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
 
     // Get total potential for entire system from MPI.
 
 #ifdef IS_MPI
 
-    MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM,
+    MPI_Allreduce(&shortRangePot_local, &potential, 1, MPI_DOUBLE, MPI_SUM,
                   MPI_COMM_WORLD);
+    potential += curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
 
 #else
 
-    potential = potential_local;
+    potential = shortRangePot_local + curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
 
 #endif // is_mpi
 
     return potential;
-}
+  }
 
-double Thermo::getTotalE() {
+  double Thermo::getTotalE() {
     double total;
 
     total = this->getKinetic() + this->getPotential();
     return total;
-}
+  }
 
-double Thermo::getTemperature() {
+  double Thermo::getTemperature() {
     
     double temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* OOPSEConstant::kb );
     return temperature;
-}
+  }
 
-double Thermo::getVolume() { 
+  double Thermo::getVolume() { 
     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
     return curSnapshot->getVolume();
-}
+  }
 
-double Thermo::getPressure() {
+  double Thermo::getPressure() {
 
     // Relies on the calculation of the full molecular pressure tensor
 
@@ -160,9 +160,26 @@ double Thermo::getPressure() {
     pressure = OOPSEConstant::pressureConvert * (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0;
 
     return pressure;
-}
+  }
 
-Mat3x3d Thermo::getPressureTensor() {
+  double Thermo::getPressure(int direction) {
+
+    // Relies on the calculation of the full molecular pressure tensor
+
+	  
+    Mat3x3d tensor;
+    double pressure;
+
+    tensor = getPressureTensor();
+
+    pressure = OOPSEConstant::pressureConvert * tensor(direction, direction);
+
+    return pressure;
+  }
+
+
+
+  Mat3x3d Thermo::getPressureTensor() {
     // returns pressure tensor in units amu*fs^-2*Ang^-1
     // routine derived via viral theorem description in:
     // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
@@ -175,13 +192,13 @@ Mat3x3d Thermo::getPressureTensor() {
     Molecule* mol;
     StuntDouble* integrableObject;    
     for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
-        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; 
-               integrableObject = mol->nextIntegrableObject(j)) {
+      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; 
+	   integrableObject = mol->nextIntegrableObject(j)) {
 
-            double mass = integrableObject->getMass();
-            Vector3d vcom = integrableObject->getVel();
-            p_local += mass * outProduct(vcom, vcom);         
-        }
+	double mass = integrableObject->getMass();
+	Vector3d vcom = integrableObject->getVel();
+	p_local += mass * outProduct(vcom, vcom);         
+      }
     }
     
 #ifdef IS_MPI
@@ -197,9 +214,9 @@ Mat3x3d Thermo::getPressureTensor() {
     pressureTensor =  (p_global + OOPSEConstant::energyConvert* tau)/volume;
 
     return pressureTensor;
-}
+  }
 
-void Thermo::saveStat(){
+  void Thermo::saveStat(){
     Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
     Stats& stat = currSnapshot->statData;
     
@@ -210,9 +227,15 @@ void Thermo::saveStat(){
     stat[Stats::PRESSURE] = getPressure();
     stat[Stats::VOLUME] = getVolume();      
 
+    Mat3x3d tensor =getPressureTensor();
+    stat[Stats::PRESSURE_TENSOR_X] = tensor(0, 0);      
+    stat[Stats::PRESSURE_TENSOR_Y] = tensor(1, 1);      
+    stat[Stats::PRESSURE_TENSOR_Z] = tensor(2, 2);      
+
+
     /**@todo need refactorying*/
     //Conserved Quantity is set by integrator and time is set by setTime
     
-}
+  }
 
 } //end namespace oopse