--- trunk/src/flucq/FluctuatingChargeConstraints.cpp	2013/07/19 21:25:45	1908
+++ trunk/src/flucq/FluctuatingChargeConstraints.cpp	2014/02/26 14:14:50	1969
@@ -40,27 +40,34 @@
  * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
  
-#include "FluctuatingChargeConstraints.hpp"
-#include "primitives/Molecule.hpp"
-
 #ifdef IS_MPI
 #include <mpi.h>
 #endif
 
+#include "FluctuatingChargeConstraints.hpp"
+#include "primitives/Molecule.hpp"
+
 namespace OpenMD {
 
   FluctuatingChargeConstraints::FluctuatingChargeConstraints(SimInfo* info) : 
-    info_(info), constrainRegions_(false), hasFlucQ_(false) {
-    
-    if (info_->usesFluctuatingCharges()) {
-      if (info_->getNFluctuatingCharges() > 0) {
-        hasFlucQ_ = true;
+    info_(info), constrainRegions_(false), hasFlucQ_(false), initialized_(false) {
+  }
+
+  void FluctuatingChargeConstraints::initialize(){
+    if(info_->usesFluctuatingCharges()){
+      if(info_->getNFluctuatingCharges() > 0){
+	hasFlucQ_ = true;
       }
     }
+    initialized_ = true;
   }
 
+
   void FluctuatingChargeConstraints::setConstrainRegions(bool cr) {
     constrainRegions_ = cr;
+
+    if (!initialized_) initialize();
+
     regionKeys_.clear();
     regionForce_.clear();
     regionCharges_.clear();
@@ -75,10 +82,13 @@ namespace OpenMD {
            mol = info_->nextMolecule(i)) {
         reg = mol->getRegion();
         if (reg >= 0) regions.insert(reg);
-        
       }
+      // resize the keys vector to the largest found value for regions.
+      regionKeys_.resize( *(regions.end()) );
+      int which = 0;
       for (std::set<int>::iterator r=regions.begin(); r!=regions.end(); ++r) {
-        regionKeys_.push_back( (*r) );
+	regionKeys_[ (*r) ] = which;
+	which++;
       }
       regionForce_.resize( regionKeys_.size() );
       regionCharges_.resize( regionKeys_.size() );
@@ -87,15 +97,15 @@ namespace OpenMD {
 
 
   void FluctuatingChargeConstraints::applyConstraints() {
+    if (!initialized_) initialize();
     if (!hasFlucQ_) return;
-    
+
     SimInfo::MoleculeIterator i;
     Molecule::FluctuatingChargeIterator  j;
     Molecule* mol;
     Atom* atom;
     
     RealType totalFrc, totalMolFrc, regionFrc, constrainedFrc;
-    
     // accumulate the total system fluctuating charge forces
     totalFrc = 0.0;
     if (constrainRegions_) {
@@ -114,8 +124,10 @@ namespace OpenMD {
         RealType frc = atom->getFlucQFrc();
         totalFrc += frc;
         if (constrainRegions_) {
-          regionForce_[regionKeys_[region]] += frc;
-          regionCharges_[regionKeys_[region]] += 1;
+	  if (region >= 0) {	    
+	    regionForce_[regionKeys_[region]] += frc;
+	    regionCharges_[regionKeys_[region]] += 1;
+	  }
         }
       }
     }
@@ -123,14 +135,14 @@ namespace OpenMD {
 #ifdef IS_MPI
     // in parallel, we need to add up the contributions from all
     // processors:
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalFrc, 1, MPI::REALTYPE, 
-                              MPI::SUM);
+    MPI_Allreduce(MPI_IN_PLACE, &totalFrc, 1, MPI_REALTYPE, 
+                  MPI_SUM, MPI_COMM_WORLD);
 
     if (constrainRegions_) {
-      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &regionForce_[0], 
-                                regionForce_.size(), MPI::REALTYPE, MPI::SUM);
-      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &regionCharges_[0], 
-                                regionCharges_.size(), MPI::INT, MPI::SUM);
+      MPI_Allreduce(MPI_IN_PLACE, &regionForce_[0], 
+                    regionForce_.size(), MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
+      MPI_Allreduce(MPI_IN_PLACE, &regionCharges_[0], 
+                    regionCharges_.size(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
     }
 
 #endif
@@ -140,7 +152,7 @@ namespace OpenMD {
     
     // do the same in the regions:
     if (constrainRegions_) {
-      for (int i = 0; i < regionForce_.size(); ++i)  {
+      for (unsigned int i = 0; i < regionForce_.size(); ++i)  {
         regionForce_[ i ] /= regionCharges_[ i ];
       }
     }
@@ -150,7 +162,10 @@ namespace OpenMD {
       
       if (constrainRegions_) {
         int region = mol->getRegion();
-        regionFrc = regionForce_[regionKeys_[region]];
+	if (region >= 0) 
+	  regionFrc = regionForce_[regionKeys_[region]];
+	else
+	  regionFrc = 0.0;
       } else {
         regionFrc = 0.0;
       }
@@ -168,9 +183,9 @@ namespace OpenMD {
 
       for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
            atom = mol->nextFluctuatingCharge(j)) {
-        //constrainedFrc = atom->getFlucQFrc() - totalFrc - totalMolFrc;
+        constrainedFrc = atom->getFlucQFrc() - totalFrc - totalMolFrc;
 
-        constrainedFrc = atom->getFlucQFrc() - totalMolFrc;
+        //constrainedFrc = atom->getFlucQFrc() - totalMolFrc;
 
         if (constrainRegions_) 
           constrainedFrc -= regionFrc;