OOPSE/libmdtools/OOPSEMinimizer.cpp

#include <math.h>
#include "OOPSEMinimizer.hpp"

OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, ForceFields* the_ff ,
                                              MinimizerParameterSet * param)
                     :RealIntegrator(theInfo, the_ff), bVerbose(false), bShake(true){
  dumpOut = NULL;
  statOut = NULL;

  tStats = new Thermo(info);

  
  paramSet = param;

  calcDim();
  
  curX = getCoor();
  curG.resize(ndim);

  preMove();
}

OOPSEMinimizer::~OOPSEMinimizer(){
  delete tStats;
  if(dumpOut)
    delete dumpOut;
  if(statOut)
    delete statOut;
  delete paramSet;
}

void OOPSEMinimizer::calcEnergyGradient(vector<double>& x, vector<double>& grad,
                                                                    double& energy, int& status){
 
  DirectionalAtom* dAtom;
  int index;
  double force[3];
  double dAtomGrad[6];
  int shakeStatus;
  
  setCoor(x);

  if (nConstrained && bShake){
    shakeStatus = shakeR();
  }
      
  calcForce(1, 1);
  
  if (nConstrained && bShake){
    shakeStatus |= shakeF();
  }
  
  x = getCoor();
  

  index = 0;

  for(int i = 0; i < integrableObjects.size(); i++){

    if (integrableObjects[i]->isDirectional()) {

      integrableObjects[i]->getGrad(dAtomGrad);

      //gradient is equal to -f
      grad[index++] = -dAtomGrad[0];
      grad[index++] = -dAtomGrad[1];
      grad[index++] = -dAtomGrad[2];
      grad[index++] = -dAtomGrad[3];
      grad[index++] = -dAtomGrad[4];
      grad[index++] = -dAtomGrad[5];

    }
    else{
      integrableObjects[i]->getFrc(force);

      grad[index++] = -force[0];
      grad[index++] = -force[1];
      grad[index++] = -force[2];

    }
    
  }
 
  energy = tStats->getPotential();

  status = shakeStatus;
}

/**
 *
 */

void OOPSEMinimizer::setCoor(vector<double>& x){

  DirectionalAtom* dAtom;
  int index;
  double position[3];
  double eulerAngle[3];


  index = 0;
  
  for(int i = 0; i < integrableObjects.size(); i++){
    
    position[0] = x[index++];
    position[1] = x[index++];
    position[2] = x[index++];

    integrableObjects[i]->setPos(position);

    if (integrableObjects[i]->isDirectional()){
 
      eulerAngle[0] = x[index++];
      eulerAngle[1] = x[index++];
      eulerAngle[2] = x[index++];

      integrableObjects[i]->setEuler(eulerAngle[0], 
                                     eulerAngle[1], 
                                     eulerAngle[2]);

    }
    
  }
  
}

/**
 *
 */
vector<double> OOPSEMinimizer::getCoor(){
 
  DirectionalAtom* dAtom;
  int index;
  double position[3];
  double eulerAngle[3];
  vector<double> x;

  x.resize(getDim());

  index = 0;
  
  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 (integrableObjects[i]->isDirectional()){

      integrableObjects[i]->getEulerAngles(eulerAngle);
      
      x[index++] = eulerAngle[0];
      x[index++] = eulerAngle[1];
      x[index++] = eulerAngle[2];
       
    }
    
  }

  return x;

}

int OOPSEMinimizer::shakeR(){
  int i, j;
  int done;
  double posA[3], posB[3];
  double velA[3], velB[3];
  double pab[3];
  double rab[3];
  int a, b, ax, ay, az, bx, by, bz;
  double rma, rmb;
  double dx, dy, dz;
  double rpab;
  double rabsq, pabsq, rpabsq;
  double diffsq;
  double gab;
  int iteration;

  for (i = 0; i < nAtoms; i++){
    moving[i] = 0;
    moved[i] = 1;
  }

  iteration = 0;
  done = 0;
  while (!done && (iteration < maxIteration)){
    done = 1;
    for (i = 0; i < nConstrained; i++){
      a = constrainedA[i];
      b = constrainedB[i];

      ax = (a * 3) + 0;
      ay = (a * 3) + 1;
      az = (a * 3) + 2;

      bx = (b * 3) + 0;
      by = (b * 3) + 1;
      bz = (b * 3) + 2;

      if (moved[a] || moved[b]){
        atoms[a]->getPos(posA);
        atoms[b]->getPos(posB);

        for (j = 0; j < 3; j++)
          pab[j] = posA[j] - posB[j];

        //periodic boundary condition

        info->wrapVector(pab);

        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];

        rabsq = constrainedDsqr[i];
        diffsq = rabsq - pabsq;

        // the original rattle code from alan tidesley
        if (fabs(diffsq) > (tol * rabsq * 2)){
          rab[0] = oldPos[ax] - oldPos[bx];
          rab[1] = oldPos[ay] - oldPos[by];
          rab[2] = oldPos[az] - oldPos[bz];

          info->wrapVector(rab);
          
          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];

          rpabsq = rpab * rpab;


          if (rpabsq < (rabsq * -diffsq)){
#ifdef IS_MPI
            a = atoms[a]->getGlobalIndex();
            b = atoms[b]->getGlobalIndex();
#endif //is_mpi
            //cerr << "Waring: constraint failure" << endl;
            gab = sqrt(rabsq/pabsq);
            rab[0] = (posA[0] - posB[0])*gab;
            rab[1]= (posA[1] - posB[1])*gab;
            rab[2] = (posA[2] - posB[2])*gab;
            
            info->wrapVector(rab);
            
            rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
            
          }

          //rma = 1.0 / atoms[a]->getMass();
          //rmb = 1.0 / atoms[b]->getMass();
          rma = 1.0;
          rmb =1.0;

          gab = diffsq / (2.0 * (rma + rmb) * rpab);

          dx = rab[0] * gab;
          dy = rab[1] * gab;
          dz = rab[2] * gab;

          posA[0] += rma * dx;
          posA[1] += rma * dy;
          posA[2] += rma * dz;

          atoms[a]->setPos(posA);

          posB[0] -= rmb * dx;
          posB[1] -= rmb * dy;
          posB[2] -= rmb * dz;

          atoms[b]->setPos(posB);

          moving[a] = 1;
          moving[b] = 1;
          done = 0;
        }
      }
    }

    for (i = 0; i < nAtoms; i++){
      moved[i] = moving[i];
      moving[i] = 0;
    }

    iteration++;
  }

  if (!done){
    cerr << "Waring: can not constraint within maxIteration" << endl;
    return -1;
  }
  else
    return 1;
}


//remove constraint force along the bond direction
int OOPSEMinimizer::shakeF(){
  int i, j;
  int done;
  double posA[3], posB[3];
  double frcA[3], frcB[3];
  double rab[3], fpab[3];
  int a, b, ax, ay, az, bx, by, bz;
  double rma, rmb;
  double rvab;
  double gab;
  double rabsq;
  double rfab;
  int iteration;

  for (i = 0; i < nAtoms; i++){
    moving[i] = 0;
    moved[i] = 1;
  }

  done = 0;
  iteration = 0;
  while (!done && (iteration < maxIteration)){
    done = 1;

    for (i = 0; i < nConstrained; i++){
      a = constrainedA[i];
      b = constrainedB[i];

      ax = (a * 3) + 0;
      ay = (a * 3) + 1;
      az = (a * 3) + 2;

      bx = (b * 3) + 0;
      by = (b * 3) + 1;
      bz = (b * 3) + 2;

      if (moved[a] || moved[b]){

        atoms[a]->getPos(posA);
        atoms[b]->getPos(posB);

        for (j = 0; j < 3; j++)
          rab[j] = posA[j] - posB[j];

        info->wrapVector(rab);

        atoms[a]->getFrc(frcA);
        atoms[b]->getFrc(frcB);

        //rma = 1.0 / atoms[a]->getMass();
        //rmb = 1.0 / atoms[b]->getMass();
        rma = 1.0;
        rmb = 1.0;
        
        
        fpab[0] = frcA[0] * rma - frcB[0] * rmb;
        fpab[1] = frcA[1] * rma - frcB[1] * rmb;
        fpab[2] = frcA[2] * rma - frcB[2] * rmb;


          gab=fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
          
          if (gab < 1.0)
            gab = 1.0;
         
          rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
          rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];

          if (fabs(rfab) > sqrt(rabsq*gab) * 0.00001){

            gab = -rfab /(rabsq*(rma + rmb));
            
            frcA[0] = rab[0] * gab;
            frcA[1] = rab[1] * gab;
            frcA[2] = rab[2] * gab;

            atoms[a]->addFrc(frcA);
            

            frcB[0] = -rab[0] * gab;
            frcB[1] = -rab[1] * gab;
            frcB[2] = -rab[2] * gab;

            atoms[b]->addFrc(frcB);
          
            moving[a] = 1;
            moving[b] = 1;
            done = 0;
          }             
      }
    }

    for (i = 0; i < nAtoms; i++){
      moved[i] = moving[i];
      moving[i] = 0;
    }

    iteration++;
  }

  if (!done){
    cerr << "Waring: can not constraint within maxIteration" << endl;
    return -1;
  }
  else
    return 1;
}

    //calculate the value of object function
void OOPSEMinimizer::calcF(){
  calcEnergyGradient(curX, curG, curF, egEvalStatus);
}
    
void OOPSEMinimizer::calcF(vector<double>& x, double&f, int& status){
  vector<double> tempG;
  tempG.resize(x.size());
  
  calcEnergyGradient(x, tempG, f, status);
}

//calculate the gradient
void OOPSEMinimizer::calcG(){
  calcEnergyGradient(curX, curG, curF, egEvalStatus);
}

void OOPSEMinimizer::calcG(vector<double>& x, vector<double>& g, double& f, int& status){
  calcEnergyGradient(x, g, f, status);
}

void OOPSEMinimizer::calcDim(){
  DirectionalAtom* dAtom;

  ndim = 0;

  for(int i = 0; i < integrableObjects.size(); i++){
    ndim += 3;
    if (integrableObjects[i]->isDirectional())
      ndim += 3;      
  }
}

void OOPSEMinimizer::setX(vector < double > & x){

    if (x.size() != ndim && bVerbose){
      //sprintf(painCave.errMsg,
      //          "OOPSEMinimizer Error: dimesion of x and curX does not match\n");
     // painCave.isFatal = 1;
     // simError();
    }

    curX = x;
}

void OOPSEMinimizer::setG(vector < double > & g){

    if (g.size() != ndim && bVerbose){
      //sprintf(painCave.errMsg,
      //          "OOPSEMinimizer Error: dimesion of g and curG does not match\n");
     // painCave.isFatal = 1;
      //simError();
    }

    curG = g;
}

void OOPSEMinimizer::writeOut(vector<double>& x, double iter){

  setX(x);

  calcG();

  dumpOut->writeDump(iter);
  statOut->writeStat(iter);
}


void OOPSEMinimizer::printMinimizerInfo(){
  cout << "--------------------------------------------------------------------" << endl;
  cout << minimizerName << endl;
  cout << "minimization parameter set" << endl;
  cout << "function tolerance = " << paramSet->getFTol() << endl;
  cout << "gradient tolerance = " << paramSet->getGTol() << endl;
  cout << "step tolerance = "<< paramSet->getFTol() << endl;
  cout << "absolute gradient tolerance = " << endl;
  cout << "max iteration = " << paramSet->getMaxIteration() << endl;
  cout << "max line search iteration = " << paramSet->getLineSearchMaxIteration() <<endl;
  cout << "shake algorithm = " << bShake << endl;
  cout << "--------------------------------------------------------------------" << endl;

}

/**
 * In thoery, we need to find the minimum along the search direction
 * 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,
 * if the right end point has lower energy, we just take it
 *
 *
 *
 */

int OOPSEMinimizer::doLineSearch(vector<double>& direction, double stepSize){
  vector<double> xa;
  vector<double> xb;
  vector<double> xc;
  vector<double> ga;
  vector<double> gb;
  vector<double> gc;
  double fa;
  double fb;
  double fc;
  double a;
  double b;
  double c;
  int status;
  double initSlope;
  double slopeA;
  double slopeB;
  double slopeC;
  bool foundLower;
  int iter;
  int maxLSIter;
  double mu;
  double eta;
  double ftol;  
  double lsTol;

  xa.resize(ndim);
  xb.resize(ndim);
  xc.resize(ndim);

  ga.resize(ndim);
  gb.resize(ndim);
  gc.resize(ndim);

  a = 0.0;
  fa =  curF;    
  xa = curX;
  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++)
    slopeA += curG[i]*direction[i];

  initSlope = slopeA;
  
  // if  going uphill, use negative gradient as searching direction 
  if (slopeA > 0) {

    if (bVerbose){
      cout << "LineSearch Warning: initial searching direction is not a descent searching direction, "
             << " use negative gradient instead. Therefore, finding a smaller vaule of function "
             << " is guaranteed"
             << endl;
    }    
    
    for (size_t i = 0; i < ndim; i++)
      direction[i] = -curG[i];    
    
    for (size_t i = 0; i < ndim; i++)
      slopeA += curG[i]*direction[i];
    
    initSlope = slopeA;
  }
    
  // Take a trial step
  for(size_t i = 0; i < ndim; i++)
    xc[i] = curX[i] + direction[i] * c;      

  calcG(xc, gc, fc, status);

  if (status < 0){
    if (bVerbose)
      cerr << "Function Evaluation Error" << endl;
  }

  //calculate the derivative at c
  slopeC = 0;
  for (size_t i = 0; i < ndim; i++)
    slopeC += gc[i]*direction[i];

  // found a lower point
  if (fc < fa) {
    curX = xc;
    curG = gc;
    curF = fc;
    return LS_SUCCEED;
  }
  else {

    if (slopeC > 0)
    stepSize *= 0.618034;
  }    

  maxLSIter = paramSet->getLineSearchMaxIteration();
   
  iter = 0;
  
  do {
    // Select a new trial point.
    // If the derivatives at points a & c have different sign we use cubic interpolate    
    //if (slopeC > 0){     
      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 
    for(size_t i = 0; i < ndim; i++)
      xb[i] = curX[i] + direction[i] * b;      

    //function evaluation
    calcG(xb, gb, fb, status);

    if (status < 0){
      if (bVerbose)
        cerr << "Function Evaluation Error" << endl;
    }

  //calculate the derivative at c
    slopeB = 0;
    for (size_t i = 0; i < ndim; i++)
      slopeB += gb[i]*direction[i];

    //Amijo Rule to stop the line search 
    if (fb <= curF +  initSlope * ftol * b) {
      curF = fb;
      curX = xb;
      curG = gb;
      return LS_SUCCEED;
     }
 
    if (slopeB <0 &&  fb < fa) {
      //replace a by b
      fa = fb;
      a = b;
      slopeA = slopeB;

      // swap coord  a/b 
      std::swap(xa, xb);
      std::swap(ga, gb);
    }
    else {
      //replace c by b
      fc = fb;
      c = b;
      slopeC = slopeB;

      // swap coord  b/c 
      std::swap(gb, gc);
      std::swap(xb, xc);
    }
     

     iter++;
  } while((fb > fa || fb > fc) && (iter < maxLSIter));     
  
  if(fb < curF || iter >= maxLSIter) {
    //could not find a lower value, we might just go uphill.      
    return LS_ERROR;
  }

  //select the end point

  if (fa <= fc) {
    curX = xa;
    curG = ga;
    curF = fa;
  }
  else {
    curX = xc;
    curG = gc;
    curF = fc;    
  }

  return LS_SUCCEED;
    
}

void OOPSEMinimizer::minimize(){

  int convgStatus;
  int stepStatus;
  int maxIter;
  //int resetFrq;
  //int nextResetIter;
  int writeFrq;
  int nextWriteIter;
  
  if (bVerbose)
    printMinimizerInfo();

  dumpOut = new DumpWriter(info);
  statOut = new StatWriter(info);
  
  init();

  //resetFrq = paramSet->getResetFrq();
  //nextResetIter = resetFrq;

  writeFrq = paramSet->getWriteFrq();
  nextWriteIter = writeFrq;
  
  maxIter = paramSet->getMaxIteration();
  
  for (curIter = 1; curIter <= maxIter; curIter++){

    stepStatus = step();

    if (stepStatus < 0){
      saveResult();
      minStatus = MIN_LSERROR;
      cerr << "OOPSEMinimizer Error: line search error, please try a small stepsize" << endl;
      return;
    }

    if (curIter == nextWriteIter){
      nextWriteIter += writeFrq;
      writeOut(curX, curIter); 
    }

    //if (curIter == nextResetIter){
    //  reset();
    //  nextResetIter += resetFrq;      
    //}

    convgStatus = checkConvg();

    if (convgStatus > 0){
      saveResult();
      minStatus = MIN_CONVERGE;
      return;
    }
    
    prepareStep();

  }


  if (bVerbose) {
    cout << "OOPSEMinimizer Warning: "
           << minimizerName << " algorithm did not converge within "
           << maxIter << " iteration" << endl;
  }
  minStatus = MIN_MAXITER;
  saveResult();
  
}
Revision:	1144
Committed:	Sat May 1 18:52:38 2004 UTC (20 years, 2 months ago) by tim
File size:	16191 byte(s)
Log Message:	C++ pass groupList to fortran
#	User	Rev	Content
1	tim	1064	#include <math.h>
2			#include "OOPSEMinimizer.hpp"
3
4			OOPSEMinimizer::OOPSEMinimizer( SimInfo theInfo, ForceFields the_ff ,
5			MinimizerParameterSet * param)
6			:RealIntegrator(theInfo, the_ff), bVerbose(false), bShake(true){
7	tim	1144	dumpOut = NULL;
8			statOut = NULL;
9	tim	1064
10			tStats = new Thermo(info);
11	tim	1144
12	tim	1064
13			paramSet = param;
14
15			calcDim();
16
17			curX = getCoor();
18			curG.resize(ndim);
19
20	tim	1074	preMove();
21	tim	1064	}
22
23			OOPSEMinimizer::~OOPSEMinimizer(){
24			delete tStats;
25	tim	1144	if(dumpOut)
26			delete dumpOut;
27			if(statOut)
28			delete statOut;
29	tim	1064	delete paramSet;
30			}
31
32			void OOPSEMinimizer::calcEnergyGradient(vector<double>& x, vector<double>& grad,
33			double& energy, int& status){
34
35			DirectionalAtom* dAtom;
36			int index;
37			double force[3];
38			double dAtomGrad[6];
39			int shakeStatus;
40
41			setCoor(x);
42
43			if (nConstrained && bShake){
44			shakeStatus = shakeR();
45			}
46
47			calcForce(1, 1);
48
49			if (nConstrained && bShake){
50			shakeStatus \|= shakeF();
51			}
52
53			x = getCoor();
54
55
56			index = 0;
57
58	gezelter	1097	for(int i = 0; i < integrableObjects.size(); i++){
59	tim	1064
60	gezelter	1097	if (integrableObjects[i]->isDirectional()) {
61	tim	1064
62	gezelter	1097	integrableObjects[i]->getGrad(dAtomGrad);
63
64	tim	1064	//gradient is equal to -f
65			grad[index++] = -dAtomGrad[0];
66			grad[index++] = -dAtomGrad[1];
67			grad[index++] = -dAtomGrad[2];
68			grad[index++] = -dAtomGrad[3];
69			grad[index++] = -dAtomGrad[4];
70			grad[index++] = -dAtomGrad[5];
71
72			}
73			else{
74	gezelter	1097	integrableObjects[i]->getFrc(force);
75	tim	1064
76			grad[index++] = -force[0];
77			grad[index++] = -force[1];
78			grad[index++] = -force[2];
79
80			}
81
82			}
83
84			energy = tStats->getPotential();
85
86			status = shakeStatus;
87			}
88
89			/**
90			*
91			*/
92
93			void OOPSEMinimizer::setCoor(vector<double>& x){
94
95			DirectionalAtom* dAtom;
96			int index;
97			double position[3];
98			double eulerAngle[3];
99
100
101			index = 0;
102
103	gezelter	1097	for(int i = 0; i < integrableObjects.size(); i++){
104	tim	1064
105			position[0] = x[index++];
106			position[1] = x[index++];
107			position[2] = x[index++];
108
109	gezelter	1097	integrableObjects[i]->setPos(position);
110	tim	1064
111	gezelter	1097	if (integrableObjects[i]->isDirectional()){
112	tim	1064
113			eulerAngle[0] = x[index++];
114			eulerAngle[1] = x[index++];
115			eulerAngle[2] = x[index++];
116
117	gezelter	1097	integrableObjects[i]->setEuler(eulerAngle[0],
118			eulerAngle[1],
119			eulerAngle[2]);
120	tim	1064
121			}
122
123			}
124
125			}
126
127			/**
128			*
129			*/
130			vector<double> OOPSEMinimizer::getCoor(){
131
132			DirectionalAtom* dAtom;
133			int index;
134			double position[3];
135			double eulerAngle[3];
136			vector<double> x;
137
138			x.resize(getDim());
139
140			index = 0;
141
142	gezelter	1097	for(int i = 0; i < integrableObjects.size(); i++){
143			integrableObjects[i]->getPos(position);
144	tim	1064
145			x[index++] = position[0];
146			x[index++] = position[1];
147			x[index++] = position[2];
148
149	gezelter	1097	if (integrableObjects[i]->isDirectional()){
150	tim	1064
151	gezelter	1097	integrableObjects[i]->getEulerAngles(eulerAngle);
152
153	tim	1064	x[index++] = eulerAngle[0];
154			x[index++] = eulerAngle[1];
155			x[index++] = eulerAngle[2];
156
157			}
158
159			}
160
161			return x;
162
163			}
164
165			int OOPSEMinimizer::shakeR(){
166			int i, j;
167			int done;
168			double posA[3], posB[3];
169			double velA[3], velB[3];
170			double pab[3];
171			double rab[3];
172			int a, b, ax, ay, az, bx, by, bz;
173			double rma, rmb;
174			double dx, dy, dz;
175			double rpab;
176			double rabsq, pabsq, rpabsq;
177			double diffsq;
178			double gab;
179			int iteration;
180
181			for (i = 0; i < nAtoms; i++){
182			moving[i] = 0;
183			moved[i] = 1;
184			}
185
186			iteration = 0;
187			done = 0;
188			while (!done && (iteration < maxIteration)){
189			done = 1;
190			for (i = 0; i < nConstrained; i++){
191			a = constrainedA[i];
192			b = constrainedB[i];
193
194			ax = (a * 3) + 0;
195			ay = (a * 3) + 1;
196			az = (a * 3) + 2;
197
198			bx = (b * 3) + 0;
199			by = (b * 3) + 1;
200			bz = (b * 3) + 2;
201
202			if (moved[a] \|\| moved[b]){
203			atoms[a]->getPos(posA);
204			atoms[b]->getPos(posB);
205
206			for (j = 0; j < 3; j++)
207			pab[j] = posA[j] - posB[j];
208
209			//periodic boundary condition
210
211			info->wrapVector(pab);
212
213			pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
214
215			rabsq = constrainedDsqr[i];
216			diffsq = rabsq - pabsq;
217
218			// the original rattle code from alan tidesley
219			if (fabs(diffsq) > (tol * rabsq * 2)){
220			rab[0] = oldPos[ax] - oldPos[bx];
221			rab[1] = oldPos[ay] - oldPos[by];
222			rab[2] = oldPos[az] - oldPos[bz];
223
224			info->wrapVector(rab);
225
226			rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
227
228			rpabsq = rpab * rpab;
229
230
231			if (rpabsq < (rabsq * -diffsq)){
232			#ifdef IS_MPI
233			a = atoms[a]->getGlobalIndex();
234			b = atoms[b]->getGlobalIndex();
235			#endif //is_mpi
236			//cerr << "Waring: constraint failure" << endl;
237			gab = sqrt(rabsq/pabsq);
238			rab[0] = (posA[0] - posB[0])*gab;
239			rab[1]= (posA[1] - posB[1])*gab;
240			rab[2] = (posA[2] - posB[2])*gab;
241
242			info->wrapVector(rab);
243
244			rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
245
246			}
247
248			//rma = 1.0 / atoms[a]->getMass();
249			//rmb = 1.0 / atoms[b]->getMass();
250			rma = 1.0;
251			rmb =1.0;
252
253			gab = diffsq / (2.0 * (rma + rmb) * rpab);
254
255			dx = rab[0] * gab;
256			dy = rab[1] * gab;
257			dz = rab[2] * gab;
258
259			posA[0] += rma * dx;
260			posA[1] += rma * dy;
261			posA[2] += rma * dz;
262
263			atoms[a]->setPos(posA);
264
265			posB[0] -= rmb * dx;
266			posB[1] -= rmb * dy;
267			posB[2] -= rmb * dz;
268
269			atoms[b]->setPos(posB);
270
271			moving[a] = 1;
272			moving[b] = 1;
273			done = 0;
274			}
275			}
276			}
277
278			for (i = 0; i < nAtoms; i++){
279			moved[i] = moving[i];
280			moving[i] = 0;
281			}
282
283			iteration++;
284			}
285
286			if (!done){
287			cerr << "Waring: can not constraint within maxIteration" << endl;
288			return -1;
289			}
290			else
291			return 1;
292			}
293
294
295			//remove constraint force along the bond direction
296			int OOPSEMinimizer::shakeF(){
297			int i, j;
298			int done;
299			double posA[3], posB[3];
300			double frcA[3], frcB[3];
301			double rab[3], fpab[3];
302			int a, b, ax, ay, az, bx, by, bz;
303			double rma, rmb;
304			double rvab;
305			double gab;
306			double rabsq;
307			double rfab;
308			int iteration;
309
310			for (i = 0; i < nAtoms; i++){
311			moving[i] = 0;
312			moved[i] = 1;
313			}
314
315			done = 0;
316			iteration = 0;
317			while (!done && (iteration < maxIteration)){
318			done = 1;
319
320			for (i = 0; i < nConstrained; i++){
321			a = constrainedA[i];
322			b = constrainedB[i];
323
324			ax = (a * 3) + 0;
325			ay = (a * 3) + 1;
326			az = (a * 3) + 2;
327
328			bx = (b * 3) + 0;
329			by = (b * 3) + 1;
330			bz = (b * 3) + 2;
331
332			if (moved[a] \|\| moved[b]){
333
334			atoms[a]->getPos(posA);
335			atoms[b]->getPos(posB);
336
337			for (j = 0; j < 3; j++)
338			rab[j] = posA[j] - posB[j];
339
340			info->wrapVector(rab);
341
342			atoms[a]->getFrc(frcA);
343			atoms[b]->getFrc(frcB);
344
345			//rma = 1.0 / atoms[a]->getMass();
346			//rmb = 1.0 / atoms[b]->getMass();
347			rma = 1.0;
348			rmb = 1.0;
349
350
351			fpab[0] = frcA[0] * rma - frcB[0] * rmb;
352			fpab[1] = frcA[1] * rma - frcB[1] * rmb;
353			fpab[2] = frcA[2] * rma - frcB[2] * rmb;
354
355
356			gab=fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
357
358			if (gab < 1.0)
359			gab = 1.0;
360
361			rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
362			rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];
363
364			if (fabs(rfab) > sqrt(rabsqgab) 0.00001){
365
366			gab = -rfab /(rabsq*(rma + rmb));
367
368			frcA[0] = rab[0] * gab;
369			frcA[1] = rab[1] * gab;
370			frcA[2] = rab[2] * gab;
371
372			atoms[a]->addFrc(frcA);
373
374
375			frcB[0] = -rab[0] * gab;
376			frcB[1] = -rab[1] * gab;
377			frcB[2] = -rab[2] * gab;
378
379			atoms[b]->addFrc(frcB);
380
381			moving[a] = 1;
382			moving[b] = 1;
383			done = 0;
384			}
385			}
386			}
387
388			for (i = 0; i < nAtoms; i++){
389			moved[i] = moving[i];
390			moving[i] = 0;
391			}
392
393			iteration++;
394			}
395
396			if (!done){
397			cerr << "Waring: can not constraint within maxIteration" << endl;
398			return -1;
399			}
400			else
401			return 1;
402			}
403
404			//calculate the value of object function
405			void OOPSEMinimizer::calcF(){
406			calcEnergyGradient(curX, curG, curF, egEvalStatus);
407			}
408
409			void OOPSEMinimizer::calcF(vector<double>& x, double&f, int& status){
410			vector<double> tempG;
411			tempG.resize(x.size());
412
413			calcEnergyGradient(x, tempG, f, status);
414			}
415
416			//calculate the gradient
417			void OOPSEMinimizer::calcG(){
418			calcEnergyGradient(curX, curG, curF, egEvalStatus);
419			}
420
421			void OOPSEMinimizer::calcG(vector<double>& x, vector<double>& g, double& f, int& status){
422			calcEnergyGradient(x, g, f, status);
423			}
424
425			void OOPSEMinimizer::calcDim(){
426			DirectionalAtom* dAtom;
427
428			ndim = 0;
429
430	gezelter	1097	for(int i = 0; i < integrableObjects.size(); i++){
431	tim	1064	ndim += 3;
432	gezelter	1097	if (integrableObjects[i]->isDirectional())
433	tim	1064	ndim += 3;
434			}
435			}
436
437			void OOPSEMinimizer::setX(vector < double > & x){
438
439			if (x.size() != ndim && bVerbose){
440			//sprintf(painCave.errMsg,
441			// "OOPSEMinimizer Error: dimesion of x and curX does not match\n");
442			// painCave.isFatal = 1;
443			// simError();
444			}
445
446			curX = x;
447			}
448
449			void OOPSEMinimizer::setG(vector < double > & g){
450
451			if (g.size() != ndim && bVerbose){
452			//sprintf(painCave.errMsg,
453			// "OOPSEMinimizer Error: dimesion of g and curG does not match\n");
454			// painCave.isFatal = 1;
455			//simError();
456			}
457
458			curG = g;
459			}
460
461			void OOPSEMinimizer::writeOut(vector<double>& x, double iter){
462
463			setX(x);
464
465			calcG();
466
467			dumpOut->writeDump(iter);
468			statOut->writeStat(iter);
469			}
470
471
472			void OOPSEMinimizer::printMinimizerInfo(){
473			cout << "--------------------------------------------------------------------" << endl;
474			cout << minimizerName << endl;
475			cout << "minimization parameter set" << endl;
476			cout << "function tolerance = " << paramSet->getFTol() << endl;
477			cout << "gradient tolerance = " << paramSet->getGTol() << endl;
478			cout << "step tolerance = "<< paramSet->getFTol() << endl;
479			cout << "absolute gradient tolerance = " << endl;
480			cout << "max iteration = " << paramSet->getMaxIteration() << endl;
481			cout << "max line search iteration = " << paramSet->getLineSearchMaxIteration() <<endl;
482			cout << "shake algorithm = " << bShake << endl;
483			cout << "--------------------------------------------------------------------" << endl;
484
485			}
486
487			/**
488			* In thoery, we need to find the minimum along the search direction
489	tim	1108	* However, function evaluation is too expensive.
490	tim	1064	* At the very begining of the problem, we check the search direction and make sure
491			* it is a descent direction
492			* we will compare the energy of two end points,
493			* if the right end point has lower energy, we just take it
494			*
495			*
496			*
497			*/
498
499			int OOPSEMinimizer::doLineSearch(vector<double>& direction, double stepSize){
500			vector<double> xa;
501			vector<double> xb;
502			vector<double> xc;
503			vector<double> ga;
504			vector<double> gb;
505			vector<double> gc;
506			double fa;
507			double fb;
508			double fc;
509			double a;
510			double b;
511			double c;
512			int status;
513			double initSlope;
514			double slopeA;
515			double slopeB;
516			double slopeC;
517			bool foundLower;
518			int iter;
519			int maxLSIter;
520			double mu;
521			double eta;
522			double ftol;
523	tim	1074	double lsTol;
524	tim	1064
525			xa.resize(ndim);
526			xb.resize(ndim);
527			xc.resize(ndim);
528
529			ga.resize(ndim);
530			gb.resize(ndim);
531			gc.resize(ndim);
532
533			a = 0.0;
534			fa = curF;
535			xa = curX;
536			ga = curG;
537			c = a + stepSize;
538			ftol = paramSet->getFTol();
539	tim	1074	lsTol = paramSet->getLineSearchTol();
540	tim	1064
541			//calculate the derivative at a = 0
542			for (size_t i = 0; i < ndim; i++)
543			slopeA += curG[i]*direction[i];
544
545			initSlope = slopeA;
546
547			// if going uphill, use negative gradient as searching direction
548			if (slopeA > 0) {
549
550			if (bVerbose){
551			cout << "LineSearch Warning: initial searching direction is not a descent searching direction, "
552			<< " use negative gradient instead. Therefore, finding a smaller vaule of function "
553			<< " is guaranteed"
554			<< endl;
555			}
556
557			for (size_t i = 0; i < ndim; i++)
558			direction[i] = -curG[i];
559
560			for (size_t i = 0; i < ndim; i++)
561			slopeA += curG[i]*direction[i];
562
563			initSlope = slopeA;
564			}
565
566			// Take a trial step
567			for(size_t i = 0; i < ndim; i++)
568			xc[i] = curX[i] + direction[i] * c;
569
570			calcG(xc, gc, fc, status);
571
572			if (status < 0){
573			if (bVerbose)
574			cerr << "Function Evaluation Error" << endl;
575			}
576
577			//calculate the derivative at c
578			slopeC = 0;
579			for (size_t i = 0; i < ndim; i++)
580			slopeC += gc[i]*direction[i];
581
582			// found a lower point
583			if (fc < fa) {
584			curX = xc;
585			curG = gc;
586			curF = fc;
587			return LS_SUCCEED;
588			}
589			else {
590
591			if (slopeC > 0)
592			stepSize *= 0.618034;
593			}
594
595			maxLSIter = paramSet->getLineSearchMaxIteration();
596
597			iter = 0;
598
599			do {
600			// Select a new trial point.
601			// If the derivatives at points a & c have different sign we use cubic interpolate
602			//if (slopeC > 0){
603			eta = 3 *(fa -fc) /(c - a) + slopeA + slopeC;
604			mu = sqrt(eta * eta - slopeA * slopeC);
605			b = a + (c - a) * (1 - (slopeC + mu - eta) /(slopeC - slopeA + 2 * mu));
606	tim	1074
607			if (b < lsTol){
608			if (bVerbose)
609			cout << "stepSize is less than line search tolerance" << endl;
610			break;
611			}
612	tim	1064	//}
613
614			// Take a trial step to this new point - new coords in xb
615			for(size_t i = 0; i < ndim; i++)
616			xb[i] = curX[i] + direction[i] * b;
617
618			//function evaluation
619			calcG(xb, gb, fb, status);
620
621			if (status < 0){
622			if (bVerbose)
623			cerr << "Function Evaluation Error" << endl;
624			}
625
626			//calculate the derivative at c
627			slopeB = 0;
628			for (size_t i = 0; i < ndim; i++)
629			slopeB += gb[i]*direction[i];
630
631			//Amijo Rule to stop the line search
632			if (fb <= curF + initSlope * ftol * b) {
633			curF = fb;
634			curX = xb;
635			curG = gb;
636			return LS_SUCCEED;
637			}
638
639			if (slopeB <0 && fb < fa) {
640			//replace a by b
641			fa = fb;
642			a = b;
643			slopeA = slopeB;
644
645			// swap coord a/b
646			std::swap(xa, xb);
647			std::swap(ga, gb);
648			}
649			else {
650			//replace c by b
651			fc = fb;
652			c = b;
653			slopeC = slopeB;
654
655			// swap coord b/c
656			std::swap(gb, gc);
657			std::swap(xb, xc);
658			}
659
660
661			iter++;
662			} while((fb > fa \|\| fb > fc) && (iter < maxLSIter));
663
664			if(fb < curF \|\| iter >= maxLSIter) {
665			//could not find a lower value, we might just go uphill.
666			return LS_ERROR;
667			}
668
669			//select the end point
670
671			if (fa <= fc) {
672			curX = xa;
673			curG = ga;
674			curF = fa;
675			}
676			else {
677			curX = xc;
678			curG = gc;
679			curF = fc;
680			}
681
682			return LS_SUCCEED;
683
684			}
685
686			void OOPSEMinimizer::minimize(){
687
688			int convgStatus;
689			int stepStatus;
690			int maxIter;
691			//int resetFrq;
692			//int nextResetIter;
693			int writeFrq;
694			int nextWriteIter;
695
696			if (bVerbose)
697			printMinimizerInfo();
698	tim	1144
699			dumpOut = new DumpWriter(info);
700			statOut = new StatWriter(info);
701	tim	1064
702			init();
703
704			//resetFrq = paramSet->getResetFrq();
705			//nextResetIter = resetFrq;
706
707			writeFrq = paramSet->getWriteFrq();
708			nextWriteIter = writeFrq;
709
710			maxIter = paramSet->getMaxIteration();
711
712			for (curIter = 1; curIter <= maxIter; curIter++){
713
714			stepStatus = step();
715
716			if (stepStatus < 0){
717			saveResult();
718			minStatus = MIN_LSERROR;
719			cerr << "OOPSEMinimizer Error: line search error, please try a small stepsize" << endl;
720			return;
721			}
722
723			if (curIter == nextWriteIter){
724			nextWriteIter += writeFrq;
725			writeOut(curX, curIter);
726			}
727
728			//if (curIter == nextResetIter){
729			// reset();
730			// nextResetIter += resetFrq;
731			//}
732
733			convgStatus = checkConvg();
734
735			if (convgStatus > 0){
736			saveResult();
737			minStatus = MIN_CONVERGE;
738			return;
739			}
740
741			prepareStep();
742
743			}
744
745
746			if (bVerbose) {
747			cout << "OOPSEMinimizer Warning: "
748			<< minimizerName << " algorithm did not converge within "
749			<< maxIter << " iteration" << endl;
750			}
751			minStatus = MIN_MAXITER;
752			saveResult();
753
754			}