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

  tStats = new Thermo(info);
  dumpOut = new DumpWriter(info);
  statOut = new StatWriter(info);
  
  paramSet = param;

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

  preMove();
}

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