--- trunk/OOPSE/libmdtools/OOPSEMinimizer.cpp	2004/03/01 20:01:50	1074
+++ trunk/OOPSE/libmdtools/OOPSEMinimizer.cpp	2004/07/16 16:29:44	1330
@@ -1,15 +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;
 
-  atoms = info->atoms;
-
   tStats = new Thermo(info);
-  dumpOut = new DumpWriter(info);
-  statOut = new StatWriter(info);
+
   
   paramSet = param;
 
@@ -18,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;
 }
 
@@ -36,30 +40,36 @@ 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();
   
 
   index = 0;
 
-  for(int i = 0; i < nAtoms; i++){
+  for(int i = 0; i < integrableObjects.size(); i++){
 
-    if(atoms[i]->isDirectional()){
-      dAtom = (DirectionalAtom*) atoms[i];
-      dAtom->getGrad(dAtomGrad);
+    if (integrableObjects[i]->isDirectional()) {
 
+      integrableObjects[i]->getGrad(dAtomGrad);
+
       //gradient is equal to -f
       grad[index++] = -dAtomGrad[0];
       grad[index++] = -dAtomGrad[1];
@@ -70,7 +80,7 @@ void OOPSEMinimizer::calcEnergyGradient(vector<double>
 
     }
     else{
-      atoms[i]->getFrc(force);
+      integrableObjects[i]->getFrc(force);
 
       grad[index++] = -force[0];
       grad[index++] = -force[1];
@@ -82,7 +92,6 @@ void OOPSEMinimizer::calcEnergyGradient(vector<double>
  
   energy = tStats->getPotential();
 
-  status = shakeStatus;
 }
 
 /**
@@ -99,22 +108,23 @@ void OOPSEMinimizer::setCoor(vector<double>& x){
 
   index = 0;
   
-  for(int i = 0; i < nAtoms; i++){
+  for(int i = 0; i < integrableObjects.size(); i++){
     
     position[0] = x[index++];
     position[1] = x[index++];
     position[2] = x[index++];
 
-    atoms[i]->setPos(position);
+    integrableObjects[i]->setPos(position);
 
-    if (atoms[i]->isDirectional()){
-      dAtom = (DirectionalAtom*) atoms[i];
+    if (integrableObjects[i]->isDirectional()){
  
       eulerAngle[0] = x[index++];
       eulerAngle[1] = x[index++];
       eulerAngle[2] = x[index++];
 
-       dAtom->setEuler(eulerAngle[0], eulerAngle[1], eulerAngle[2]);
+      integrableObjects[i]->setEuler(eulerAngle[0], 
+                                     eulerAngle[1], 
+                                     eulerAngle[2]);
 
     }
     
@@ -137,17 +147,17 @@ vector<double> OOPSEMinimizer::getCoor(){
 
   index = 0;
   
-  for(int i = 0; i < nAtoms; i++){
-    atoms[i]->getPos(position);
+  for(int i = 0; i < integrableObjects.size(); i++){
+    integrableObjects[i]->getPos(position);
 
     x[index++] = position[0];
     x[index++] = position[1];
     x[index++] = position[2];
 
-    if (atoms[i]->isDirectional()){
-      dAtom = (DirectionalAtom*) atoms[i];
-      dAtom->getEulerAngles(eulerAngle);
+    if (integrableObjects[i]->isDirectional()){
 
+      integrableObjects[i]->getEulerAngles(eulerAngle);
+      
       x[index++] = eulerAngle[0];
       x[index++] = eulerAngle[1];
       x[index++] = eulerAngle[2];
@@ -160,6 +170,7 @@ vector<double> OOPSEMinimizer::getCoor(){
 
 }
 
+/*
 int OOPSEMinimizer::shakeR(){
   int i, j;
   int done;
@@ -399,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);
 }
@@ -425,9 +438,9 @@ void OOPSEMinimizer::calcDim(){
 
   ndim = 0;
 
-  for(int i = 0; i < nAtoms; i++){
+  for(int i = 0; i < integrableObjects.size(); i++){
     ndim += 3;
-    if (atoms[i]->isDirectional())
+    if (integrableObjects[i]->isDirectional())
       ndim += 3;      
   }
 }
@@ -484,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,
@@ -537,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];
 
@@ -693,6 +707,9 @@ void OOPSEMinimizer::minimize(){
   
   if (bVerbose)
     printMinimizerInfo();
+
+  dumpOut = new DumpWriter(info);
+  statOut = new StatWriter(info);
   
   init();
 
@@ -708,6 +725,9 @@ void OOPSEMinimizer::minimize(){
 
     stepStatus = step();
 
+    if (bShake)
+      shakeAlgo->doPreConstraint();
+    
     if (stepStatus < 0){
       saveResult();
       minStatus = MIN_LSERROR;
@@ -720,11 +740,6 @@ void OOPSEMinimizer::minimize(){
       writeOut(curX, curIter); 
     }
 
-    //if (curIter == nextResetIter){
-    //  reset();
-    //  nextResetIter += resetFrq;      
-    //}
-
     convgStatus = checkConvg();
 
     if (convgStatus > 0){