--- trunk/src/restraints/ThermoIntegrationForceManager.cpp	2006/06/19 01:36:06	990
+++ trunk/src/restraints/ThermoIntegrationForceManager.cpp	2009/09/07 16:31:51	1360
@@ -39,60 +39,53 @@
  * such damages.
  */
  
-#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>
+#endif
+
 namespace oopse {
   
   ThermoIntegrationForceManager::ThermoIntegrationForceManager(SimInfo* info): 
-    ForceManager(info){
-      currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
-      simParam = info_->getSimParams();
+    RestraintForceManager(info){
+    currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
+    simParam = info_->getSimParams();
     
-      if (simParam->haveThermodynamicIntegrationLambda()){
-	tIntLambda_ = simParam->getThermodynamicIntegrationLambda();
-      }
-      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->haveThermodynamicIntegrationK()){
-	tIntK_ = simParam->getThermodynamicIntegrationK();
-      }
-      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->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_);
-
+    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();      
+    }
+    
+    // build the scaling factor used to modulate the forces and torques
+    factor_ = pow(tIntLambda_, tIntK_);
+  }
   
   ThermoIntegrationForceManager::~ThermoIntegrationForceManager(){
   }
@@ -131,7 +124,7 @@ namespace oopse {
         }
       }
     }
-  
+    
     // set vraw to be the unmodulated potential
     lrPot_ = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
     curSnapshot->statData[Stats::VRAW] = lrPot_;
@@ -145,24 +138,31 @@ namespace oopse {
     tempTau *= factor_;
     curSnapshot->statData.setTau(tempTau);
 
-//     sprintf(painCave.errMsg, "Before Calc_Restraint_Forces\n");
-//     painCave.isFatal = 0;
-//     simError();
-  
-    // do crystal restraint forces for thermodynamic integration
-    if (simParam->getUseSolidThermInt()) {
-      
-      lrPot_ += restraint_->Calc_Restraint_Forces();
-      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
-      
-      vHarm_ = restraint_->getVharm();
-      curSnapshot->statData[Stats::VHARM] = vHarm_;
+    // now, on to the applied restraining potentials (if needed):
+    RealType restPot_local = 0.0;
+    RealType vHarm_local = 0.0;
+    
+    if (simParam->getUseRestraints()) {
+      // do restraints from RestraintForceManager:
+      //restPot_local = doRestraints(1.0 - factor_);     
+      restPot_local = doRestraints(1.0 - factor_);      
+      vHarm_local = getUnscaledPotential();
     }
+      
+#ifdef IS_MPI
+    RealType restPot;
+    MPI::COMM_WORLD.Allreduce(&restPot_local, &restPot, 1, 
+                              MPI::REALTYPE, MPI::SUM);
+    MPI::COMM_WORLD.Allreduce(&vHarm_local, &vHarm_, 1, 
+                              MPI::REALTYPE, MPI::SUM);         
+    lrPot_ += restPot;
+#else
+    lrPot_ += restPot_local;
+    vHarm_ = vHarm_local;
+#endif
 
-//     sprintf(painCave.errMsg, "After Calc_Restraint_Forces\n");
-//     painCave.isFatal = 0;
-//     simError();
-    
-  }
-  
+    // give the final values to stats
+    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
+    curSnapshot->statData[Stats::VHARM] = vHarm_;
+  }  
 }