--- trunk/OOPSE/libmdtools/OOPSEMinimizer.cpp	2004/02/24 15:44:45	1064
+++ trunk/OOPSE/libmdtools/OOPSEMinimizer.cpp	2004/05/01 18:52:38	1144
@@ -4,12 +4,11 @@ OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, Forc
 OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, ForceFields* the_ff ,
                                               MinimizerParameterSet * param)
                      :RealIntegrator(theInfo, the_ff), bVerbose(false), bShake(true){
+  dumpOut = NULL;
+  statOut = NULL;
 
-  atoms = info->atoms;
-
   tStats = new Thermo(info);
-  dumpOut = new DumpWriter(info);
-  statOut = new StatWriter(info);
+
   
   paramSet = param;
 
@@ -18,12 +17,15 @@ OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, Forc
   curX = getCoor();
   curG.resize(ndim);
 
+  preMove();
 }
 
 OOPSEMinimizer::~OOPSEMinimizer(){
   delete tStats;
-  delete dumpOut;
-  delete statOut;
+  if(dumpOut)
+    delete dumpOut;
+  if(statOut)
+    delete statOut;
   delete paramSet;
 }
 
@@ -53,12 +55,12 @@ void OOPSEMinimizer::calcEnergyGradient(vector<double>
 
   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];
@@ -69,7 +71,7 @@ void OOPSEMinimizer::calcEnergyGradient(vector<double>
 
     }
     else{
-      atoms[i]->getFrc(force);
+      integrableObjects[i]->getFrc(force);
 
       grad[index++] = -force[0];
       grad[index++] = -force[1];
@@ -98,22 +100,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]);
 
     }
     
@@ -136,17 +139,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];
@@ -424,9 +427,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;      
   }
 }
@@ -483,7 +486,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,
@@ -517,6 +520,7 @@ int OOPSEMinimizer::doLineSearch(vector<double>& direc
   double mu;
   double eta;
   double ftol;  
+  double lsTol;
 
   xa.resize(ndim);
   xb.resize(ndim);
@@ -532,6 +536,7 @@ int OOPSEMinimizer::doLineSearch(vector<double>& direc
   ga = curG;
   c = a + stepSize; 
   ftol = paramSet->getFTol();
+  lsTol = paramSet->getLineSearchTol();
         
   //calculate the derivative at a = 0 
   for (size_t i = 0; i < ndim; i++)
@@ -598,6 +603,12 @@ int OOPSEMinimizer::doLineSearch(vector<double>& direc
       eta = 3 *(fa -fc) /(c - a) + slopeA + slopeC;
       mu = sqrt(eta * eta - slopeA * slopeC);      
       b = a + (c - a) * (1 - (slopeC + mu - eta) /(slopeC - slopeA + 2 * mu));       
+
+      if (b < lsTol){
+        if (bVerbose)
+          cout << "stepSize is less than line search tolerance" << endl;
+        break;        
+      }
     //}
 
     // Take a trial step to this new point - new coords in xb 
@@ -684,6 +695,9 @@ void OOPSEMinimizer::minimize(){
   
   if (bVerbose)
     printMinimizerInfo();
+
+  dumpOut = new DumpWriter(info);
+  statOut = new StatWriter(info);
   
   init();