--- trunk/src/restraints/RestraintForceManager.cpp	2009/11/25 20:02:06	1390
+++ trunk/src/restraints/RestraintForceManager.cpp	2014/10/16 19:13:51	2024
@@ -35,11 +35,18 @@
  *                                                                      
  * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
  * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
- * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
- * [4]  Vardeman & Gezelter, in progress (2009).                        
+ * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
  
+#ifdef IS_MPI
+#include <mpi.h>
+#endif
+
+#include "config.h"
 #include <cmath>
+
 #include "restraints/RestraintForceManager.hpp"
 #include "restraints/MolecularRestraint.hpp"
 #include "restraints/ObjectRestraint.hpp"
@@ -49,11 +56,7 @@
 #include "utils/StringUtils.hpp"
 #include "selection/SelectionEvaluator.hpp"
 #include "selection/SelectionManager.hpp"
-#ifdef IS_MPI
-#include <mpi.h>
-#endif
 
-
 namespace OpenMD {
 
   RestraintForceManager::RestraintForceManager(SimInfo* info): ForceManager(info) {
@@ -89,6 +92,8 @@ namespace OpenMD {
     int nRestraintStamps = simParam->getNRestraintStamps();
     std::vector<RestraintStamp*> stamp = simParam->getRestraintStamps();
 
+    std::vector<int> stuntDoubleIndex;
+
     for (int i = 0; i < nRestraintStamps; i++){
 
       std::string myType = toUpperCopy(stamp[i]->getType());
@@ -96,7 +101,6 @@ namespace OpenMD {
       if (myType.compare("MOLECULAR")==0){
 
         int molIndex;
-        std::vector<Vector3d> ref;
         Vector3d refCom;
 
         if (!stamp[i]->haveMolIndex()) {
@@ -132,8 +136,11 @@ namespace OpenMD {
           // this proc doesn't have the molecule.  Do a quick check to
           // make sure another processor is supposed to have it.
          
-          int myrank = MPI::COMM_WORLD.Get_rank();
+          int myrank;
+          MPI_Comm_rank( MPI_COMM_WORLD, &myrank);
+
           if (info_->getMolToProc(molIndex) == myrank) {
+         
             // If we were supposed to have it but got a null, then freak out.
 #endif
 
@@ -149,6 +156,16 @@ namespace OpenMD {
 #endif
         }
         
+
+#ifdef IS_MPI
+        // only handle this molecular restraint if this processor owns the
+        // molecule
+        int myrank;
+        MPI_Comm_rank( MPI_COMM_WORLD, &myrank);
+        if (info_->getMolToProc(molIndex) == myrank) {
+
+#endif
+
         MolecularRestraint* rest = new MolecularRestraint();
 
         std::string restPre("mol_");
@@ -184,7 +201,9 @@ namespace OpenMD {
         restraints_.push_back(rest);
         mol->addProperty(new RestraintData("Restraint", rest));
         restrainedMols_.push_back(mol);
-        
+#ifdef IS_MPI
+        }
+#endif        
       } else if (myType.compare("OBJECT") == 0) {
         
         std::string objectSelection;
@@ -206,12 +225,18 @@ namespace OpenMD {
         evaluator.loadScriptString(objectSelection);
         seleMan.setSelectionSet(evaluator.evaluate());        
         int selectionCount = seleMan.getSelectionCount();
-        
+
+#ifdef IS_MPI    
+        MPI_Allreduce(MPI_IN_PLACE, &selectionCount, 1, MPI_INT, MPI_SUM, 
+                      MPI_COMM_WORLD);
+#endif
+                
         sprintf(painCave.errMsg,
                 "Restraint Info: The specified restraint objectSelection,\n"
                 "\t\t%s\n"
                 "\twill result in %d integrable objects being\n"
                 "\trestrained.\n", objectSelection.c_str(), selectionCount);
+        painCave.severity = OPENMD_INFO;
         painCave.isFatal = 0;
         simError();                     
 
@@ -220,7 +245,8 @@ namespace OpenMD {
 
         for (sd = seleMan.beginSelected(selei); sd != NULL; 
              sd = seleMan.nextSelected(selei)) {
-          
+          stuntDoubleIndex.push_back(sd->getGlobalIntegrableObjectIndex());
+
           ObjectRestraint* rest = new ObjectRestraint();
           
           if (stamp[i]->haveDisplacementSpringConstant()) {
@@ -260,16 +286,14 @@ namespace OpenMD {
     // are times when it won't use restraints at all, so only open the
     // restraint file if we are actually using restraints:
 
-    if (simParam->getUseRestraints()) {
+    if (simParam->getUseRestraints()) { 
       std::string refFile = simParam->getRestraint_file();
-      RestReader* rr = new RestReader(info, refFile);
-
+      RestReader* rr = new RestReader(info, refFile, stuntDoubleIndex);
       rr->readReferenceStructure();
     }
 
     restOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".rest";
     restOut = new RestWriter(info_, restOutput_.c_str(), restraints_);
-    
     if(!restOut){
       sprintf(painCave.errMsg, "Restraint error: Failed to create RestWriter\n");
       painCave.isFatal = 1;
@@ -292,22 +316,20 @@ namespace OpenMD {
     currRestTime_ = currSnapshot_->getTime();
   }
 
-  void RestraintForceManager::calcForces(bool needPotential, bool needStress){
+  void RestraintForceManager::calcForces(){
 
-    ForceManager::calcForces(needPotential, needStress);    
-    RealType restPot_local, restPot;
+    ForceManager::calcForces();    
+    RealType restPot(0.0);
 
-    restPot_local = doRestraints(1.0);
+    restPot = doRestraints(1.0);
 
 #ifdef IS_MPI    
-    MPI::COMM_WORLD.Allreduce(&restPot_local, &restPot, 1, 
-                              MPI::REALTYPE, MPI::SUM);
-#else
-    restPot = restPot_local;
+    MPI_Allreduce(MPI_IN_PLACE, &restPot, 1, MPI_REALTYPE, MPI_SUM, 
+                  MPI_COMM_WORLD);
 #endif
+
     currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
-    currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += restPot;
-    currSnapshot_->statData[Stats::VHARM] = restPot;
+    currSnapshot_->setRestraintPotential(restPot);
 
     //write out forces and current positions of restrained molecules    
     if (currSnapshot_->getTime() >= currRestTime_){
@@ -322,7 +344,6 @@ namespace OpenMD {
     Molecule::IntegrableObjectIterator ioi;
     MolecularRestraint* mRest;
     StuntDouble* sd;
-    RealType pTot;
 
     std::vector<StuntDouble*>::const_iterator ro;
     ObjectRestraint* oRest;
@@ -372,9 +393,10 @@ namespace OpenMD {
 
       std::vector<Vector3d> struc;
       std::vector<Vector3d> forces;
+
       
       for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL; 
-          sd = (*rm)->nextIntegrableObject(ioi)) {        
+          sd = (*rm)->nextIntegrableObject(ioi)) {
         struc.push_back(sd->getPos());
       }
 
@@ -390,12 +412,11 @@ namespace OpenMD {
         index++;
       }
       
-      unscaledPotential_ += mRest->getUnscaledPotential();
+      unscaledPotential_ += mRest->getUnscaledPotential();      
 
-      restInfo = mRest->getRestraintInfo();
-
       // only collect data on restraints that we're going to print:
       if (mRest->getPrintRestraint()) 
+        restInfo = mRest->getRestraintInfo();
         restInfo_.push_back(restInfo);
     }
 
@@ -432,7 +453,6 @@ namespace OpenMD {
       
       // phew.  At this point, we should have the pointer to the
       // correct Object restraint in the variable oRest.
-
       oRest->setScaleFactor(scalingFactor);
       
       Vector3d pos = (*ro)->getPos();
@@ -456,10 +476,9 @@ namespace OpenMD {
 
       unscaledPotential_ += oRest->getUnscaledPotential();
 
-      restInfo = oRest->getRestraintInfo();
-
       // only collect data on restraints that we're going to print:
       if (oRest->getPrintRestraint()) 
+        restInfo = oRest->getRestraintInfo();      
         restInfo_.push_back(restInfo);
     }