--- trunk/src/flucq/FluctuatingChargeConstraints.cpp	2013/06/16 15:15:42	1879
+++ trunk/src/flucq/FluctuatingChargeConstraints.cpp	2014/02/26 14:14:50	1969
@@ -40,63 +40,138 @@
  * [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), 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();
+
+    if (constrainRegions_) {
+      SimInfo::MoleculeIterator i;
+      Molecule* mol;
+      int reg;
+      std::set<int> regions;      
+
+      for (mol = info_->beginMolecule(i); mol != NULL; 
+           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_[ (*r) ] = which;
+	which++;
+      }
+      regionForce_.resize( regionKeys_.size() );
+      regionCharges_.resize( regionKeys_.size() );
+    }
+  }
+
+
   void FluctuatingChargeConstraints::applyConstraints() {
+    if (!initialized_) initialize();
     if (!hasFlucQ_) return;
-    
+
     SimInfo::MoleculeIterator i;
     Molecule::FluctuatingChargeIterator  j;
     Molecule* mol;
     Atom* atom;
     
-    RealType totalFrc, totalMolFrc, constrainedFrc;
-    
+    RealType totalFrc, totalMolFrc, regionFrc, constrainedFrc;
     // accumulate the total system fluctuating charge forces
     totalFrc = 0.0;
+    if (constrainRegions_) {
+      std::fill(regionForce_.begin(), regionForce_.end(), 0.0); 
+      std::fill(regionCharges_.begin(), regionCharges_.end(), 0); 
+    }
 
     for (mol = info_->beginMolecule(i); mol != NULL; 
          mol = info_->nextMolecule(i)) {
 
+      int region = mol->getRegion();
+          
       for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
            atom = mol->nextFluctuatingCharge(j)) {
-        totalFrc += atom->getFlucQFrc();
-      }
 
+        RealType frc = atom->getFlucQFrc();
+        totalFrc += frc;
+        if (constrainRegions_) {
+	  if (region >= 0) {	    
+	    regionForce_[regionKeys_[region]] += frc;
+	    regionCharges_[regionKeys_[region]] += 1;
+	  }
+        }
+      }
     }
 
 #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_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
- 
+
     // divide by the total number of fluctuating charges:
     totalFrc /= info_->getNFluctuatingCharges();
+    
+    // do the same in the regions:
+    if (constrainRegions_) {
+      for (unsigned int i = 0; i < regionForce_.size(); ++i)  {
+        regionForce_[ i ] /= regionCharges_[ i ];
+      }
+    }
 
     for (mol = info_->beginMolecule(i); mol != NULL; 
          mol = info_->nextMolecule(i)) {     
       
-      totalMolFrc = 0.0;
+      if (constrainRegions_) {
+        int region = mol->getRegion();
+	if (region >= 0) 
+	  regionFrc = regionForce_[regionKeys_[region]];
+	else
+	  regionFrc = 0.0;
+      } else {
+        regionFrc = 0.0;
+      }
 
+      totalMolFrc = 0.0;
+      
       // molecular constraints can be done with a second loop.
       if (mol->constrainTotalCharge()) {
         for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
@@ -108,8 +183,13 @@ namespace OpenMD {
 
       for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
            atom = mol->nextFluctuatingCharge(j)) {
-        //constrainedFrc = atom->getFlucQFrc() - totalFrc - totalMolFrc;
-        constrainedFrc = atom->getFlucQFrc() - totalMolFrc;
+        constrainedFrc = atom->getFlucQFrc() - totalFrc - totalMolFrc;
+
+        //constrainedFrc = atom->getFlucQFrc() - totalMolFrc;
+
+        if (constrainRegions_) 
+          constrainedFrc -= regionFrc;
+        
         atom->setFlucQFrc(constrainedFrc);
       }      
     }