--- trunk/src/restraints/ThermoIntegrationForceManager.cpp	2005/03/12 19:05:16	437
+++ trunk/src/restraints/ThermoIntegrationForceManager.cpp	2009/09/07 16:31:51	1360
@@ -1,4 +1,4 @@
- /*
+/*
  * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
  *
  * The University of Notre Dame grants you ("Licensee") a
@@ -39,59 +39,52 @@
  * 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){
+    RestraintForceManager(info){
     currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
     simParam = info_->getSimParams();
     
-    if (simParam->haveThermIntLambda()){
-      tIntLambda_ = simParam->getThermIntLambda();
+    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",
+              "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();
+    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",
+              "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_);
-      
-    }
-    
     // 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,16 +138,31 @@ namespace oopse {
     tempTau *= factor_;
     curSnapshot->statData.setTau(tempTau);
 
-    // 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
+
+    // give the final values to stats
+    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
+    curSnapshot->statData[Stats::VHARM] = vHarm_;
+  }  
 }