OOPSE/libmdtools/OOPSEMinimizer.cpp

#include <math.h>
#include "OOPSEMinimizer.hpp"
#include "ShakeMin.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);

  shakeAlgo = new ShakeMinFramework(theInfo);
  shakeAlgo->doPreConstraint();
}

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;

  status = 1;
  
  setCoor(x);

  if (bShake){
    shakeStatus = shakeAlgo->doShakeR();
    if(shakeStatus < 0)
      status = -1;
  }

  calcForce(1, 1);

  if (bShake){
    shakeStatus = shakeAlgo->doShakeF();
    if(shakeStatus < 0)
      status = -1;
  }

  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();

}

/**
 *
 */

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 
  slopeA = 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 (bShake)
      shakeAlgo->doPreConstraint();
    
    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); 
    }

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