--- branches/development/src/flucq/FluctuatingChargePropagator.cpp	2012/06/18 18:23:20	1756
+++ trunk/src/flucq/FluctuatingChargePropagator.cpp	2014/04/14 18:32:51	1981
@@ -35,7 +35,7 @@
  *                                                                      
  * [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).          
+ * [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).
  */
@@ -48,8 +48,6 @@
 #include "optimization/StatusFunction.hpp"
 #include "optimization/OptimizationFactory.hpp"
 
-
-
 #ifdef IS_MPI
 #include <mpi.h>
 #endif
@@ -57,76 +55,72 @@ namespace OpenMD {
 using namespace QuantLib;
 namespace OpenMD {
 
-  FluctuatingChargePropagator::FluctuatingChargePropagator(SimInfo* info, 
-                                                           ForceManager* fm) : 
-    info_(info), forceMan_(fm), hasFlucQ_(false) {
+  FluctuatingChargePropagator::FluctuatingChargePropagator(SimInfo* info) : 
+    info_(info), hasFlucQ_(false), forceMan_(NULL), initialized_(false) {
     
+    Globals* simParams = info_->getSimParams();
+    fqParams_ = simParams->getFluctuatingChargeParameters();    
+  }
+
+  FluctuatingChargePropagator::~FluctuatingChargePropagator() {
+  }
+
+  void FluctuatingChargePropagator::setForceManager(ForceManager* forceMan) {
+    forceMan_ = forceMan;
+  }
+
+  void FluctuatingChargePropagator::initialize() {
     if (info_->usesFluctuatingCharges()) {
       if (info_->getNFluctuatingCharges() > 0) {
-        
         hasFlucQ_ = true;
-        Globals* simParams = info_->getSimParams();
-        fqParams_ = simParams->getFluctuatingChargeParameters();
 	fqConstraints_ = new FluctuatingChargeConstraints(info_);
-        
+	fqConstraints_->setConstrainRegions(fqParams_->getConstrainRegions());
       }
     }
-  }
+    
+    if (!hasFlucQ_) {
+      initialized_ = true;
+      return;
+    }
 
-  void FluctuatingChargePropagator::initialize() {
-
-    if (!hasFlucQ_) return;
-
     SimInfo::MoleculeIterator i;
     Molecule::FluctuatingChargeIterator  j;
     Molecule* mol;
     Atom* atom;
     
-    for (mol = info_->beginMolecule(i); mol != NULL; 
-         mol = info_->nextMolecule(i)) {
-      for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
-           atom = mol->nextFluctuatingCharge(j)) {
-        atom->setFlucQPos(0.0);
-        atom->setFlucQVel(0.0);
-      }
-    }
-
-    std::cerr << "doing a minimization\n";
-    
-    fqConstraints_ = new FluctuatingChargeConstraints(info_);
-    FluctuatingChargeObjectiveFunction flucQobjf(info_, forceMan_, fqConstraints_);    
-    DynamicVector<RealType> initCoords = flucQobjf.setInitialCoords();
-    Problem problem(flucQobjf, *(new NoConstraint()), *(new NoStatus()), initCoords);
-    EndCriteria endCriteria(1000, 100, 1e-5, 1e-5, 1e-5);       
-    OptimizationMethod* minim = OptimizationFactory::getInstance()->createOptimization("SD", info_);
-
-    DumpStatusFunction dsf(info_);                           // we want a dump file written every iteration
-
-    minim->minimize(problem, endCriteria);
-    cerr << "Finished minimization\n";
-    for (mol = info_->beginMolecule(i); mol != NULL; 
-         mol = info_->nextMolecule(i)) {
-      for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
-           atom = mol->nextFluctuatingCharge(j)) {
-	cerr << atom->getType() << "\tQ Pos: " << atom->getFlucQPos() << "\n";
-      }
-    }
-    // std::cerr << "after minim\n";
+    // For single-minima flucq, this ensures a net neutral system, but
+    // for multiple minima, this is no longer the right thing to do:
+    //
     // for (mol = info_->beginMolecule(i); mol != NULL; 
     //      mol = info_->nextMolecule(i)) {
     //   for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
     //        atom = mol->nextFluctuatingCharge(j)) {
-    //     cerr << "q = " << atom->getFlucQPos(0.0) << "\n";
+    //     atom->setFlucQPos(0.0);
+    //     atom->setFlucQVel(0.0);
     //   }
     // }
+    
+    FluctuatingChargeObjectiveFunction flucQobjf(info_, forceMan_, 
+                                                 fqConstraints_);
 
+    DynamicVector<RealType> initCoords = flucQobjf.setInitialCoords();
+    Problem problem(flucQobjf, *(new NoConstraint()), *(new NoStatus()), 
+                    initCoords);
 
-  }
+    EndCriteria endCriteria(1000, 100, 1e-5, 1e-5, 1e-5);       
 
+    OptimizationMethod* minim = OptimizationFactory::getInstance()->createOptimization("SD", info_);
 
+    DumpStatusFunction dsf(info_);  // we want a dump file written
+                                    // every iteration
+    minim->minimize(problem, endCriteria);
+    cerr << "back from minim\n";
+    initialized_ = true;
+  }
+
   void FluctuatingChargePropagator::applyConstraints() {
+    if (!initialized_) initialize();
     if (!hasFlucQ_) return;
-
     fqConstraints_->applyConstraints();
   }
 }