--- trunk/src/brains/Thermo.cpp	2005/04/15 22:04:00	507
+++ trunk/src/brains/Thermo.cpp	2008/09/12 20:51:22	1292
@@ -53,7 +53,7 @@ namespace oopse {
 
 namespace oopse {
 
-  double Thermo::getKinetic() {
+  RealType Thermo::getKinetic() {
     SimInfo::MoleculeIterator miter;
     std::vector<StuntDouble*>::iterator iiter;
     Molecule* mol;
@@ -64,18 +64,19 @@ namespace oopse {
     int i;
     int j;
     int k;
-    double kinetic = 0.0;
-    double kinetic_global = 0.0;
+    RealType mass;
+    RealType kinetic = 0.0;
+    RealType 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)) {
-
-	double mass = integrableObject->getMass();
-	Vector3d vel = integrableObject->getVel();
-
+        
+	mass = integrableObject->getMass();
+	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();
@@ -96,7 +97,7 @@ namespace oopse {
     
 #ifdef IS_MPI
 
-    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_DOUBLE, MPI_SUM,
+    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
                   MPI_COMM_WORLD);
     kinetic = kinetic_global;
 
@@ -107,53 +108,53 @@ namespace oopse {
     return kinetic;
   }
 
-  double Thermo::getPotential() {
-    double potential = 0.0;
+  RealType Thermo::getPotential() {
+    RealType potential = 0.0;
     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
-    double potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] + 
-      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
+    RealType 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_REALTYPE, 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 total;
+  RealType Thermo::getTotalE() {
+    RealType total;
 
     total = this->getKinetic() + this->getPotential();
     return total;
   }
 
-  double Thermo::getTemperature() {
+  RealType Thermo::getTemperature() {
     
-    double temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* OOPSEConstant::kb );
+    RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* OOPSEConstant::kb );
     return temperature;
   }
 
-  double Thermo::getVolume() { 
+  RealType Thermo::getVolume() { 
     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
     return curSnapshot->getVolume();
   }
 
-  double Thermo::getPressure() {
+  RealType Thermo::getPressure() {
 
     // Relies on the calculation of the full molecular pressure tensor
 
 
     Mat3x3d tensor;
-    double pressure;
+    RealType pressure;
 
     tensor = getPressureTensor();
 
@@ -162,6 +163,21 @@ namespace oopse {
     return pressure;
   }
 
+  RealType Thermo::getPressure(int direction) {
+
+    // Relies on the calculation of the full molecular pressure tensor
+
+	  
+    Mat3x3d tensor;
+    RealType 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:
@@ -178,27 +194,28 @@ namespace oopse {
       for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; 
 	   integrableObject = mol->nextIntegrableObject(j)) {
 
-	double mass = integrableObject->getMass();
+	RealType mass = integrableObject->getMass();
 	Vector3d vcom = integrableObject->getVel();
 	p_local += mass * outProduct(vcom, vcom);         
       }
     }
     
 #ifdef IS_MPI
-    MPI_Allreduce(p_local.getArrayPointer(), p_global.getArrayPointer(), 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(p_local.getArrayPointer(), p_global.getArrayPointer(), 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
 #else
     p_global = p_local;
 #endif // is_mpi
 
-    double volume = this->getVolume();
+    RealType volume = this->getVolume();
     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
     Mat3x3d tau = curSnapshot->statData.getTau();
 
     pressureTensor =  (p_global + OOPSEConstant::energyConvert* tau)/volume;
-
+    
     return pressureTensor;
   }
 
+
   void Thermo::saveStat(){
     Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
     Stats& stat = currSnapshot->statData;
@@ -210,6 +227,72 @@ namespace oopse {
     stat[Stats::PRESSURE] = getPressure();
     stat[Stats::VOLUME] = getVolume();      
 
+    Mat3x3d tensor =getPressureTensor();
+    stat[Stats::PRESSURE_TENSOR_XX] = tensor(0, 0);      
+    stat[Stats::PRESSURE_TENSOR_XY] = tensor(0, 1);      
+    stat[Stats::PRESSURE_TENSOR_XZ] = tensor(0, 2);      
+    stat[Stats::PRESSURE_TENSOR_YX] = tensor(1, 0);      
+    stat[Stats::PRESSURE_TENSOR_YY] = tensor(1, 1);      
+    stat[Stats::PRESSURE_TENSOR_YZ] = tensor(1, 2);      
+    stat[Stats::PRESSURE_TENSOR_ZX] = tensor(2, 0);      
+    stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);      
+    stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);      
+
+
+    Globals* simParams = info_->getSimParams();
+
+    if (simParams->haveTaggedAtomPair() && 
+        simParams->havePrintTaggedPairDistance()) {
+      if ( simParams->getPrintTaggedPairDistance()) {
+        
+        std::pair<int, int> tap = simParams->getTaggedAtomPair();
+        Vector3d pos1, pos2, rab;
+
+#ifdef IS_MPI        
+
+	int mol1 = info_->getGlobalMolMembership(tap.first);
+	int mol2 = info_->getGlobalMolMembership(tap.second);
+        
+        int proc1 = info_->getMolToProc(mol1);
+        int proc2 = info_->getMolToProc(mol2);
+
+	RealType data[3];
+        if (proc1 == worldRank) {
+          StuntDouble* sd1 = info_->getIOIndexToIntegrableObject(tap.first);
+          pos1 = sd1->getPos();
+          data[0] = pos1.x();
+          data[1] = pos1.y();
+          data[2] = pos1.z();          
+          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
+        } else {
+          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
+          pos1 = Vector3d(data);
+        }
+
+
+        if (proc2 == worldRank) {
+          StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second);
+          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);
+        } else {
+          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
+          pos2 = Vector3d(data);
+        }
+#else
+        StuntDouble* at1 = info_->getIOIndexToIntegrableObject(tap.first);
+        StuntDouble* at2 = info_->getIOIndexToIntegrableObject(tap.second);
+        pos1 = at1->getPos();
+        pos2 = at2->getPos();
+#endif        
+        rab = pos2 - pos1;
+        currSnapshot->wrapVector(rab);
+        stat[Stats::TAGGED_PAIR_DISTANCE] =  rab.length();
+      }
+    }
+      
     /**@todo need refactorying*/
     //Conserved Quantity is set by integrator and time is set by setTime