--- trunk/src/restraints/ThermoIntegrationForceManager.cpp	2005/03/10 15:10:24	417
+++ trunk/src/restraints/ThermoIntegrationForceManager.cpp	2006/08/31 20:12:12	1028
@@ -1,4 +1,4 @@
- /*
+/*
  * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
  *
  * The University of Notre Dame grants you ("Licensee") a
@@ -42,56 +42,64 @@
 #include <cmath>
 #include "restraints/ThermoIntegrationForceManager.hpp"
 #include "integrators/Integrator.hpp"
+#include "math/SquareMatrix3.hpp"
 #include "primitives/Molecule.hpp"
 #include "utils/simError.h"
 #include "utils/OOPSEConstant.hpp"
 #include "utils/StringUtils.hpp"
 
+#ifdef IS_MPI
+#include <mpi.h>
+#define TAKE_THIS_TAG_REAL 2
+#endif //is_mpi
+
 namespace oopse {
   
   ThermoIntegrationForceManager::ThermoIntegrationForceManager(SimInfo* info): 
-  ForceManager(info){
-    currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
-    simParam = info_->getSimParams();
+    ForceManager(info){
+      currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
+      simParam = info_->getSimParams();
     
-    if (simParam->haveThermIntLambda()){
-      tIntLambda_ = simParam->getThermIntLambda();
-    }
-    else{
-      tIntLambda_ = 1.0;
-      sprintf(painCave.errMsg,
-              "ThermoIntegration error: the transformation parameter (lambda) was\n"
-              "\tnot specified. OOPSE will use a default value of %f. To set\n"
-              "\tlambda, use the thermodynamicIntegrationLambda variable.\n",
-              tIntLambda_);
-      painCave.isFatal = 0;
-      simError();
-    }
+      if (simParam->haveThermodynamicIntegrationLambda()){
+	tIntLambda_ = simParam->getThermodynamicIntegrationLambda();
+      }
+      else{
+	tIntLambda_ = 1.0;
+	sprintf(painCave.errMsg,
+		"ThermoIntegration error: the transformation parameter\n"
+		"\t(lambda) was not specified. OOPSE will use a default\n"
+		"\tvalue of %f. To set lambda, use the \n"
+		"\tthermodynamicIntegrationLambda variable.\n",
+		tIntLambda_);
+	painCave.isFatal = 0;
+	simError();
+      }
     
-    if (simParam->haveThermIntK()){
-      tIntK_ = simParam->getThermIntK();
-    }
-    else{
-      tIntK_ = 1.0;
-      sprintf(painCave.errMsg,
-              "ThermoIntegration Warning: the tranformation parameter exponent\n"
-              "\t(k) was not specified. OOPSE will use a default value of %f.\n"
-              "\tTo set k, use the thermodynamicIntegrationK variable.\n",
-              tIntK_);
-      painCave.isFatal = 0;
-      simError();      
-    }
+      if (simParam->haveThermodynamicIntegrationK()){
+	tIntK_ = simParam->getThermodynamicIntegrationK();
+      }
+      else{
+	tIntK_ = 1.0;
+	sprintf(painCave.errMsg,
+		"ThermoIntegration Warning: the tranformation parameter\n"
+		"\texponent (k) was not specified. OOPSE will use a default\n"
+		"\tvalue of %f. To set k, use the thermodynamicIntegrationK\n"
+		"\tvariable.\n",
+		tIntK_);
+	painCave.isFatal = 0;
+	simError();      
+      }
     
-    if (simParam->getUseSolidThermInt()) {
-      // build a restraint object
-      restraint_ =  new Restraints(info_, tIntLambda_, tIntK_);
+      if (simParam->getUseSolidThermInt()) {
+	// build a restraint object
+	restraint_ =  new Restraints(info_, tIntLambda_, tIntK_);
       
-    }
+      }
     
-    // build the scaling factor used to modulate the forces and torques
-    factor_ = pow(tIntLambda_, tIntK_);
-    
-  }
+      // build the scaling factor used to modulate the forces and torques
+      factor_ = pow(tIntLambda_, tIntK_);
+
+    }
   
   ThermoIntegrationForceManager::~ThermoIntegrationForceManager(){
   }
@@ -105,12 +113,13 @@ namespace oopse {
     StuntDouble* integrableObject;
     Vector3d frc;
     Vector3d trq;
+    Mat3x3d tempTau;
     
-    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
-    
     // perform the standard calcForces first
     ForceManager::calcForces(needPotential, needStress);
     
+    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
+
     // now scale forces and torques of all the integrableObjects
       
     for (mol = info_->beginMolecule(mi); mol != NULL; 
@@ -128,26 +137,77 @@ namespace oopse {
           integrableObject->setTrq(trq);
         }
       }
-      
-      // set vraw to be the unmodulated potential
-      lrPot_ = curSnapshot->statData[Stats::POTENTIAL_ENERGY];
-      curSnapshot->statData[Stats::VRAW] = lrPot_;
-      
-      // modulate the potential and update the snapshot
-      lrPot_ *= factor_;
-      curSnapshot->statData[Stats::POTENTIAL_ENERGY] = lrPot_;
     }
+  
+    // set vraw to be the unmodulated potential
+    lrPot_ = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
+    curSnapshot->statData[Stats::VRAW] = lrPot_;
     
-    // do crystal restraint forces for thermodynamic integration
+    // modulate the potential and update the snapshot
+    lrPot_ *= factor_;
+    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
+    
+    // scale the pressure tensor
+    tempTau = curSnapshot->statData.getTau();
+    tempTau *= factor_;
+    curSnapshot->statData.setTau(tempTau);
+#ifndef IS_MPI
+    // do the single processor crystal restraint forces for 
+    // thermodynamic integration
     if (simParam->getUseSolidThermInt()) {
       
       lrPot_ += restraint_->Calc_Restraint_Forces();
-      curSnapshot->statData[Stats::POTENTIAL_ENERGY] = lrPot_;
+      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
       
       vHarm_ = restraint_->getVharm();
       curSnapshot->statData[Stats::VHARM] = vHarm_;
     }
-    
+#else
+    double tempLRPot = 0.0;
+    double tempVHarm = 0.0;
+    MPI_Status ierr;
+    int nproc;
+    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
+
+    // do the MPI crystal restraint forces for each processor
+    if (simParam->getUseSolidThermInt()) {
+      tempLRPot = restraint_->Calc_Restraint_Forces();
+      tempVHarm = restraint_->getVharm();
+    }
+
+    // master receives and accumulates the restraint info
+    if (worldRank == 0) {
+      for(int i = 0 ; i < nproc; ++i) {
+	if (i == worldRank) {
+	  lrPot_ += tempLRPot;
+	  vHarm_ += tempVHarm;
+	} else {
+	  MPI_Recv(&tempLRPot, 1, MPI_REALTYPE, i, 
+		   TAKE_THIS_TAG_REAL, MPI_COMM_WORLD, &ierr);
+	  MPI_Recv(&tempVHarm, 1, MPI_REALTYPE, i, 
+		   TAKE_THIS_TAG_REAL, MPI_COMM_WORLD, &ierr);
+	  lrPot_ += tempLRPot;
+	  vHarm_ += tempVHarm;
+	}
+      }
+
+      // give the final values to stats
+      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
+      curSnapshot->statData[Stats::VHARM] = vHarm_;
+
+    } else {
+      // pack up and send the appropriate info to the master node
+      for(int j = 1; j < nproc; ++j) {
+     	if (worldRank == j) {
+
+	  MPI_Send(&tempLRPot, 1, MPI_REALTYPE, 0, 
+		   TAKE_THIS_TAG_REAL, MPI_COMM_WORLD);
+	  MPI_Send(&tempVHarm, 1, MPI_REALTYPE, 0, 
+		   TAKE_THIS_TAG_REAL, MPI_COMM_WORLD);
+	}
+      }
+    }
+#endif //is_mpi
   }
   
 }