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

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

  shakeAlgo->doShakeR();
      
  calcForce(1, 1);
  
  //if (nConstrained && bShake){
  //  shakeStatus |= shakeF();
  //}

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