OOPSE/libmdtools/OOPSEMinimizer.cpp

#include <math.h>
#include "OOPSEMinimizer.hpp"
#include "ShakeMin.hpp"
#include "Integrator.cpp"

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