--- trunk/OOPSE/libmdtools/OOPSEMinimizer.cpp	2004/04/12 20:32:20	1097
+++ trunk/OOPSE/libmdtools/OOPSEMinimizer.cpp	2004/07/16 16:29:44	1330
@@ -1,13 +1,16 @@
 #include <math.h>
 #include "OOPSEMinimizer.hpp"
+#include "ShakeMin.hpp"
+#include "Integrator.cpp"
 
 OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, ForceFields* the_ff ,
-                                              MinimizerParameterSet * param)
-                     :RealIntegrator(theInfo, the_ff), bVerbose(false), bShake(true){
+                                              MinimizerParameterSet * param) :
+              RealIntegrator(theInfo, the_ff), bShake(true), bVerbose(false) {
+  dumpOut = NULL;
+  statOut = NULL;
 
   tStats = new Thermo(info);
-  dumpOut = new DumpWriter(info);
-  statOut = new StatWriter(info);
+
   
   paramSet = param;
 
@@ -16,13 +19,16 @@ OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, Forc
   curX = getCoor();
   curG.resize(ndim);
 
-  preMove();
+  shakeAlgo = new ShakeMinFramework(theInfo);
+  shakeAlgo->doPreConstraint();
 }
 
 OOPSEMinimizer::~OOPSEMinimizer(){
   delete tStats;
-  delete dumpOut;
-  delete statOut;
+  if(dumpOut)
+    delete dumpOut;
+  if(statOut)
+    delete statOut;
   delete paramSet;
 }
 
@@ -34,19 +40,25 @@ void OOPSEMinimizer::calcEnergyGradient(vector<double>
   double force[3];
   double dAtomGrad[6];
   int shakeStatus;
+
+  status = 1;
   
   setCoor(x);
 
-  if (nConstrained && bShake){
-    shakeStatus = shakeR();
+  if (bShake){
+    shakeStatus = shakeAlgo->doShakeR();
+    if(shakeStatus < 0)
+      status = -1;
   }
-      
+
   calcForce(1, 1);
-  
-  if (nConstrained && bShake){
-    shakeStatus |= shakeF();
+
+  if (bShake){
+    shakeStatus = shakeAlgo->doShakeF();
+    if(shakeStatus < 0)
+      status = -1;
   }
-  
+
   x = getCoor();
   
 
@@ -80,7 +92,6 @@ void OOPSEMinimizer::calcEnergyGradient(vector<double>
  
   energy = tStats->getPotential();
 
-  status = shakeStatus;
 }
 
 /**
@@ -159,6 +170,7 @@ vector<double> OOPSEMinimizer::getCoor(){
 
 }
 
+/*
 int OOPSEMinimizer::shakeR(){
   int i, j;
   int done;
@@ -398,7 +410,9 @@ int OOPSEMinimizer::shakeF(){
     return 1;
 }
 
-    //calculate the value of object function
+*/
+
+//calculate the value of object function
 void OOPSEMinimizer::calcF(){
   calcEnergyGradient(curX, curG, curF, egEvalStatus);
 }
@@ -483,7 +497,7 @@ void OOPSEMinimizer::printMinimizerInfo(){
 
 /**
  * In thoery, we need to find the minimum along the search direction
- * However, function evaluation is too expensive. I
+ * However, function evaluation is too expensive. 
  * At the very begining of the problem, we check the search direction and make sure
  * it is a descent direction
  * we will compare the energy of two end points,
@@ -536,6 +550,7 @@ int OOPSEMinimizer::doLineSearch(vector<double>& direc
   lsTol = paramSet->getLineSearchTol();
         
   //calculate the derivative at a = 0 
+  slopeA = 0;
   for (size_t i = 0; i < ndim; i++)
     slopeA += curG[i]*direction[i];
 
@@ -692,6 +707,9 @@ void OOPSEMinimizer::minimize(){
   
   if (bVerbose)
     printMinimizerInfo();
+
+  dumpOut = new DumpWriter(info);
+  statOut = new StatWriter(info);
   
   init();
 
@@ -707,6 +725,9 @@ void OOPSEMinimizer::minimize(){
 
     stepStatus = step();
 
+    if (bShake)
+      shakeAlgo->doPreConstraint();
+    
     if (stepStatus < 0){
       saveResult();
       minStatus = MIN_LSERROR;
@@ -719,11 +740,6 @@ void OOPSEMinimizer::minimize(){
       writeOut(curX, curIter); 
     }
 
-    //if (curIter == nextResetIter){
-    //  reset();
-    //  nextResetIter += resetFrq;      
-    //}
-
     convgStatus = checkConvg();
 
     if (convgStatus > 0){