--- trunk/src/flucq/FluctuatingChargeConstraints.cpp	2013/07/19 21:25:45	1908
+++ trunk/src/flucq/FluctuatingChargeConstraints.cpp	2013/07/31 19:30:46	1920
@@ -50,17 +50,24 @@ namespace OpenMD {
 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;
+	  }
         }
       }
     }
@@ -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;
       }