src/minimizers/OOPSEMinimizer.cpp

#include <math.h>

#include "minimizers/OOPSEMinimizer.hpp"

#include "integrators/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);


  preMove();

}


OOPSEMinimizer::~OOPSEMinimizer(){

  delete tStats;

  if(dumpOut)

    delete dumpOut;

  if(statOut)

    delete statOut;

  delete paramSet;

}


void OOPSEMinimizer::calcEnergyGradient(vector<double>& x,  std::vector<double>& grad,

                                                                    double& energy, int& status){

 
  DirectionalAtom* dAtom;

  int index;

  double force[3];

  double dAtomGrad[6];

  int shakeStatus;


  status = 1;

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


}


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


    }

    
  }

  
}


 std::vector<double> OOPSEMinimizer::getCoor(){

 
  DirectionalAtom* dAtom;

  int index;

  double position[3];

  double eulerAngle[3];

   std::vector<double> x;


  x.resize(getDim());


  index = 0;

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

    position = integrableObjects[i]->getPos();


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

        posA = atoms[a]->getPos();

        posB = atoms[b]->getPos();


        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

            //std::cerr << "Waring: constraint failure" << std::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){

    std::cerr << "Waring: can not constraint within maxIteration" << std::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]){


        posA = atoms[a]->getPos();

        posB = atoms[b]->getPos();


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

    std::cerr << "Waring: can not constraint within maxIteration" << std::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){

   std::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,  std::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 << "--------------------------------------------------------------------" << std::endl;

  cout << minimizerName << std::endl;

  cout << "minimization parameter set" << std::endl;

  cout << "function tolerance = " << paramSet->getFTol() << std::endl;

  cout << "gradient tolerance = " << paramSet->getGTol() << std::endl;

  cout << "step tolerance = "<< paramSet->getFTol() << std::endl;

  cout << "absolute gradient tolerance = " << std::endl;

  cout << "max iteration = " << paramSet->getMaxIteration() << std::endl;

  cout << "max line search iteration = " << paramSet->getLineSearchMaxIteration() <<std::endl;

  cout << "shake algorithm = " << bShake << std::endl;

  cout << "--------------------------------------------------------------------" << std::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){

   std::vector<double> xa;

   std::vector<double> xb;

   std::vector<double> xc;

   std::vector<double> ga;

   std::vector<double> gb;

   std::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"

             << std::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)

      std::cerr << "Function Evaluation Error" << std::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" << std::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)

        std::cerr << "Function Evaluation Error" << std::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 writeFrq;

  int nextWriteIter;

  
  if (bVerbose)

    printMinimizerInfo();


  dumpOut = new DumpWriter(info);

  statOut = new StatWriter(info);

  
  init();


  writeFrq = paramSet->getWriteFrq();

  nextWriteIter = writeFrq;

  
  maxIter = paramSet->getMaxIteration();

  
  for (curIter = 1; curIter <= maxIter; curIter++){


    stepStatus = step();


    if (nConstrained && bShake)

      preMove();

    
    if (stepStatus < 0){

      saveResult();

      minStatus = MIN_LSERROR;

      std::cerr << "OOPSEMinimizer Error: line search error, please try a small stepsize" << std::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" << std::endl;

  }

  minStatus = MIN_MAXITER;

  saveResult();

  
}

Revision:	1830
Committed:	Thu Dec 2 05:17:10 2004 UTC (19 years, 8 months ago) by tim
File size:	17000 byte(s)
Log Message:	change endl to std::endl
#	User	Rev	Content
1	tim	1826	#include <math.h>
2
3			#include "minimizers/OOPSEMinimizer.hpp"
4
5			#include "integrators/Integrator.cpp"
6
7
8
9			OOPSEMinimizer::OOPSEMinimizer( SimInfo theInfo, ForceFields the_ff ,
10
11			MinimizerParameterSet * param) :
12
13			RealIntegrator(theInfo, the_ff), bShake(true), bVerbose(false) {
14
15			dumpOut = NULL;
16
17			statOut = NULL;
18
19
20
21			tStats = new Thermo(info);
22
23
24
25
26
27			paramSet = param;
28
29
30
31			calcDim();
32
33
34
35			curX = getCoor();
36
37			curG.resize(ndim);
38
39
40
41			preMove();
42
43			}
44
45
46
47			OOPSEMinimizer::~OOPSEMinimizer(){
48
49			delete tStats;
50
51			if(dumpOut)
52
53			delete dumpOut;
54
55			if(statOut)
56
57			delete statOut;
58
59			delete paramSet;
60
61			}
62
63
64
65			void OOPSEMinimizer::calcEnergyGradient(vector<double>& x, std::vector<double>& grad,
66
67			double& energy, int& status){
68
69
70
71			DirectionalAtom* dAtom;
72
73			int index;
74
75			double force[3];
76
77			double dAtomGrad[6];
78
79			int shakeStatus;
80
81
82
83			status = 1;
84
85
86
87			setCoor(x);
88
89
90
91			if (nConstrained && bShake){
92
93			shakeStatus = shakeR();
94
95			}
96
97
98
99			calcForce(1, 1);
100
101
102
103			if (nConstrained && bShake){
104
105			shakeStatus = shakeF();
106
107			}
108
109
110
111			x = getCoor();
112
113
114
115
116
117			index = 0;
118
119
120
121			for(int i = 0; i < integrableObjects.size(); i++){
122
123
124
125			if (integrableObjects[i]->isDirectional()) {
126
127
128
129			integrableObjects[i]->getGrad(dAtomGrad);
130
131
132
133			//gradient is equal to -f
134
135			grad[index++] = -dAtomGrad[0];
136
137			grad[index++] = -dAtomGrad[1];
138
139			grad[index++] = -dAtomGrad[2];
140
141			grad[index++] = -dAtomGrad[3];
142
143			grad[index++] = -dAtomGrad[4];
144
145			grad[index++] = -dAtomGrad[5];
146
147
148
149			}
150
151			else{
152
153			integrableObjects[i]->getFrc(force);
154
155
156
157			grad[index++] = -force[0];
158
159			grad[index++] = -force[1];
160
161			grad[index++] = -force[2];
162
163
164
165			}
166
167
168
169			}
170
171
172
173			energy = tStats->getPotential();
174
175
176
177			}
178
179
180
181			void OOPSEMinimizer::setCoor(vector<double>& x){
182
183
184
185			DirectionalAtom* dAtom;
186
187			int index;
188
189			double position[3];
190
191			double eulerAngle[3];
192
193
194
195
196
197			index = 0;
198
199
200
201			for(int i = 0; i < integrableObjects.size(); i++){
202
203
204
205			position[0] = x[index++];
206
207			position[1] = x[index++];
208
209			position[2] = x[index++];
210
211
212
213			integrableObjects[i]->setPos(position);
214
215
216
217			if (integrableObjects[i]->isDirectional()){
218
219
220
221			eulerAngle[0] = x[index++];
222
223			eulerAngle[1] = x[index++];
224
225			eulerAngle[2] = x[index++];
226
227
228
229			integrableObjects[i]->setEuler(eulerAngle[0],
230
231			eulerAngle[1],
232
233			eulerAngle[2]);
234
235
236
237			}
238
239
240
241			}
242
243
244
245			}
246
247
248
249			std::vector<double> OOPSEMinimizer::getCoor(){
250
251
252
253			DirectionalAtom* dAtom;
254
255			int index;
256
257			double position[3];
258
259			double eulerAngle[3];
260
261			std::vector<double> x;
262
263
264
265			x.resize(getDim());
266
267
268
269			index = 0;
270
271
272
273			for(int i = 0; i < integrableObjects.size(); i++){
274
275			position = integrableObjects[i]->getPos();
276
277
278
279			x[index++] = position[0];
280
281			x[index++] = position[1];
282
283			x[index++] = position[2];
284
285
286
287			if (integrableObjects[i]->isDirectional()){
288
289
290
291			integrableObjects[i]->getEulerAngles(eulerAngle);
292
293
294
295			x[index++] = eulerAngle[0];
296
297			x[index++] = eulerAngle[1];
298
299			x[index++] = eulerAngle[2];
300
301
302
303			}
304
305
306
307			}
308
309
310
311			return x;
312
313
314
315			}
316
317
318
319
320
321			int OOPSEMinimizer::shakeR(){
322
323			int i, j;
324
325			int done;
326
327			double posA[3], posB[3];
328
329			double velA[3], velB[3];
330
331			double pab[3];
332
333			double rab[3];
334
335			int a, b, ax, ay, az, bx, by, bz;
336
337			double rma, rmb;
338
339			double dx, dy, dz;
340
341			double rpab;
342
343			double rabsq, pabsq, rpabsq;
344
345			double diffsq;
346
347			double gab;
348
349			int iteration;
350
351
352
353			for (i = 0; i < nAtoms; i++){
354
355			moving[i] = 0;
356
357			moved[i] = 1;
358
359			}
360
361
362
363			iteration = 0;
364
365			done = 0;
366
367			while (!done && (iteration < maxIteration)){
368
369			done = 1;
370
371			for (i = 0; i < nConstrained; i++){
372
373			a = constrainedA[i];
374
375			b = constrainedB[i];
376
377
378
379			ax = (a * 3) + 0;
380
381			ay = (a * 3) + 1;
382
383			az = (a * 3) + 2;
384
385
386
387			bx = (b * 3) + 0;
388
389			by = (b * 3) + 1;
390
391			bz = (b * 3) + 2;
392
393
394
395			if (moved[a] \|\| moved[b]){
396
397			posA = atoms[a]->getPos();
398
399			posB = atoms[b]->getPos();
400
401
402
403			for (j = 0; j < 3; j++)
404
405			pab[j] = posA[j] - posB[j];
406
407
408
409			//periodic boundary condition
410
411
412
413			info->wrapVector(pab);
414
415
416
417			pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
418
419
420
421			rabsq = constrainedDsqr[i];
422
423			diffsq = rabsq - pabsq;
424
425
426
427			// the original rattle code from alan tidesley
428
429			if (fabs(diffsq) > (tol * rabsq * 2)){
430
431			rab[0] = oldPos[ax] - oldPos[bx];
432
433			rab[1] = oldPos[ay] - oldPos[by];
434
435			rab[2] = oldPos[az] - oldPos[bz];
436
437
438
439			info->wrapVector(rab);
440
441
442
443			rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
444
445
446
447			rpabsq = rpab * rpab;
448
449
450
451
452
453			if (rpabsq < (rabsq * -diffsq)){
454
455			#ifdef IS_MPI
456
457			a = atoms[a]->getGlobalIndex();
458
459			b = atoms[b]->getGlobalIndex();
460
461			#endif //is_mpi
462
463	tim	1830	//std::cerr << "Waring: constraint failure" << std::endl;
464	tim	1826
465			gab = sqrt(rabsq/pabsq);
466
467			rab[0] = (posA[0] - posB[0])*gab;
468
469			rab[1]= (posA[1] - posB[1])*gab;
470
471			rab[2] = (posA[2] - posB[2])*gab;
472
473
474
475			info->wrapVector(rab);
476
477
478
479			rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
480
481
482
483			}
484
485
486
487			//rma = 1.0 / atoms[a]->getMass();
488
489			//rmb = 1.0 / atoms[b]->getMass();
490
491			rma = 1.0;
492
493			rmb =1.0;
494
495
496
497			gab = diffsq / (2.0 * (rma + rmb) * rpab);
498
499
500
501			dx = rab[0] * gab;
502
503			dy = rab[1] * gab;
504
505			dz = rab[2] * gab;
506
507
508
509			posA[0] += rma * dx;
510
511			posA[1] += rma * dy;
512
513			posA[2] += rma * dz;
514
515
516
517			atoms[a]->setPos(posA);
518
519
520
521			posB[0] -= rmb * dx;
522
523			posB[1] -= rmb * dy;
524
525			posB[2] -= rmb * dz;
526
527
528
529			atoms[b]->setPos(posB);
530
531
532
533			moving[a] = 1;
534
535			moving[b] = 1;
536
537			done = 0;
538
539			}
540
541			}
542
543			}
544
545
546
547			for (i = 0; i < nAtoms; i++){
548
549			moved[i] = moving[i];
550
551			moving[i] = 0;
552
553			}
554
555
556
557			iteration++;
558
559			}
560
561
562
563			if (!done){
564
565	tim	1830	std::cerr << "Waring: can not constraint within maxIteration" << std::endl;
566	tim	1826
567			return -1;
568
569			}
570
571			else
572
573			return 1;
574
575			}
576
577
578
579
580
581			//remove constraint force along the bond direction
582
583			int OOPSEMinimizer::shakeF(){
584
585			int i, j;
586
587			int done;
588
589			double posA[3], posB[3];
590
591			double frcA[3], frcB[3];
592
593			double rab[3], fpab[3];
594
595			int a, b, ax, ay, az, bx, by, bz;
596
597			double rma, rmb;
598
599			double rvab;
600
601			double gab;
602
603			double rabsq;
604
605			double rfab;
606
607			int iteration;
608
609
610
611			for (i = 0; i < nAtoms; i++){
612
613			moving[i] = 0;
614
615			moved[i] = 1;
616
617			}
618
619
620
621			done = 0;
622
623			iteration = 0;
624
625			while (!done && (iteration < maxIteration)){
626
627			done = 1;
628
629
630
631			for (i = 0; i < nConstrained; i++){
632
633			a = constrainedA[i];
634
635			b = constrainedB[i];
636
637
638
639			ax = (a * 3) + 0;
640
641			ay = (a * 3) + 1;
642
643			az = (a * 3) + 2;
644
645
646
647			bx = (b * 3) + 0;
648
649			by = (b * 3) + 1;
650
651			bz = (b * 3) + 2;
652
653
654
655			if (moved[a] \|\| moved[b]){
656
657
658
659			posA = atoms[a]->getPos();
660
661			posB = atoms[b]->getPos();
662
663
664
665			for (j = 0; j < 3; j++)
666
667			rab[j] = posA[j] - posB[j];
668
669
670
671			info->wrapVector(rab);
672
673
674
675			atoms[a]->getFrc(frcA);
676
677			atoms[b]->getFrc(frcB);
678
679
680
681			//rma = 1.0 / atoms[a]->getMass();
682
683			//rmb = 1.0 / atoms[b]->getMass();
684
685			rma = 1.0;
686
687			rmb = 1.0;
688
689
690
691
692
693			fpab[0] = frcA[0] * rma - frcB[0] * rmb;
694
695			fpab[1] = frcA[1] * rma - frcB[1] * rmb;
696
697			fpab[2] = frcA[2] * rma - frcB[2] * rmb;
698
699
700
701
702
703			gab=fpab[0] * fpab[0] + fpab[1] * fpab[1] + fpab[2] * fpab[2];
704
705
706
707			if (gab < 1.0)
708
709			gab = 1.0;
710
711
712
713			rabsq = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
714
715			rfab = rab[0] * fpab[0] + rab[1] * fpab[1] + rab[2] * fpab[2];
716
717
718
719			if (fabs(rfab) > sqrt(rabsqgab) 0.00001){
720
721
722
723			gab = -rfab /(rabsq*(rma + rmb));
724
725
726
727			frcA[0] = rab[0] * gab;
728
729			frcA[1] = rab[1] * gab;
730
731			frcA[2] = rab[2] * gab;
732
733
734
735			atoms[a]->addFrc(frcA);
736
737
738
739
740
741			frcB[0] = -rab[0] * gab;
742
743			frcB[1] = -rab[1] * gab;
744
745			frcB[2] = -rab[2] * gab;
746
747
748
749			atoms[b]->addFrc(frcB);
750
751
752
753			moving[a] = 1;
754
755			moving[b] = 1;
756
757			done = 0;
758
759			}
760
761			}
762
763			}
764
765
766
767			for (i = 0; i < nAtoms; i++){
768
769			moved[i] = moving[i];
770
771			moving[i] = 0;
772
773			}
774
775
776
777			iteration++;
778
779			}
780
781
782
783			if (!done){
784
785	tim	1830	std::cerr << "Waring: can not constraint within maxIteration" << std::endl;
786	tim	1826
787			return -1;
788
789			}
790
791			else
792
793			return 1;
794
795			}
796
797
798
799
800
801
802
803			//calculate the value of object function
804
805			void OOPSEMinimizer::calcF(){
806
807			calcEnergyGradient(curX, curG, curF, egEvalStatus);
808
809			}
810
811
812
813			void OOPSEMinimizer::calcF(vector<double>& x, double&f, int& status){
814
815			std::vector<double> tempG;
816
817			tempG.resize(x.size());
818
819
820
821			calcEnergyGradient(x, tempG, f, status);
822
823			}
824
825
826
827			//calculate the gradient
828
829			void OOPSEMinimizer::calcG(){
830
831			calcEnergyGradient(curX, curG, curF, egEvalStatus);
832
833			}
834
835
836
837			void OOPSEMinimizer::calcG(vector<double>& x, std::vector<double>& g, double& f, int& status){
838
839			calcEnergyGradient(x, g, f, status);
840
841			}
842
843
844
845			void OOPSEMinimizer::calcDim(){
846
847			DirectionalAtom* dAtom;
848
849
850
851			ndim = 0;
852
853
854
855			for(int i = 0; i < integrableObjects.size(); i++){
856
857			ndim += 3;
858
859			if (integrableObjects[i]->isDirectional())
860
861			ndim += 3;
862
863			}
864
865			}
866
867
868
869			void OOPSEMinimizer::setX(vector < double > & x){
870
871
872
873			if (x.size() != ndim && bVerbose){
874
875			//sprintf(painCave.errMsg,
876
877			// "OOPSEMinimizer Error: dimesion of x and curX does not match\n");
878
879			// painCave.isFatal = 1;
880
881			// simError();
882
883			}
884
885
886
887			curX = x;
888
889			}
890
891
892
893			void OOPSEMinimizer::setG(vector < double > & g){
894
895
896
897			if (g.size() != ndim && bVerbose){
898
899			//sprintf(painCave.errMsg,
900
901			// "OOPSEMinimizer Error: dimesion of g and curG does not match\n");
902
903			// painCave.isFatal = 1;
904
905			//simError();
906
907			}
908
909
910
911			curG = g;
912
913			}
914
915
916
917			void OOPSEMinimizer::writeOut(vector<double>& x, double iter){
918
919
920
921			setX(x);
922
923
924
925			calcG();
926
927
928
929			dumpOut->writeDump(iter);
930
931			statOut->writeStat(iter);
932
933			}
934
935
936
937
938
939			void OOPSEMinimizer::printMinimizerInfo(){
940
941	tim	1830	cout << "--------------------------------------------------------------------" << std::endl;
942	tim	1826
943	tim	1830	cout << minimizerName << std::endl;
944	tim	1826
945	tim	1830	cout << "minimization parameter set" << std::endl;
946	tim	1826
947	tim	1830	cout << "function tolerance = " << paramSet->getFTol() << std::endl;
948	tim	1826
949	tim	1830	cout << "gradient tolerance = " << paramSet->getGTol() << std::endl;
950	tim	1826
951	tim	1830	cout << "step tolerance = "<< paramSet->getFTol() << std::endl;
952	tim	1826
953	tim	1830	cout << "absolute gradient tolerance = " << std::endl;
954	tim	1826
955	tim	1830	cout << "max iteration = " << paramSet->getMaxIteration() << std::endl;
956	tim	1826
957	tim	1830	cout << "max line search iteration = " << paramSet->getLineSearchMaxIteration() <<std::endl;
958	tim	1826
959	tim	1830	cout << "shake algorithm = " << bShake << std::endl;
960	tim	1826
961	tim	1830	cout << "--------------------------------------------------------------------" << std::endl;
962	tim	1826
963
964
965			}
966
967
968
969			/**
970
971			* In thoery, we need to find the minimum along the search direction
972
973			* However, function evaluation is too expensive.
974
975			* At the very begining of the problem, we check the search direction and make sure
976
977			* it is a descent direction
978
979			* we will compare the energy of two end points,
980
981			* if the right end point has lower energy, we just take it
982
983			*
984
985			*
986
987			*
988
989			*/
990
991
992
993			int OOPSEMinimizer::doLineSearch(vector<double>& direction, double stepSize){
994
995			std::vector<double> xa;
996
997			std::vector<double> xb;
998
999			std::vector<double> xc;
1000
1001			std::vector<double> ga;
1002
1003			std::vector<double> gb;
1004
1005			std::vector<double> gc;
1006
1007			double fa;
1008
1009			double fb;
1010
1011			double fc;
1012
1013			double a;
1014
1015			double b;
1016
1017			double c;
1018
1019			int status;
1020
1021			double initSlope;
1022
1023			double slopeA;
1024
1025			double slopeB;
1026
1027			double slopeC;
1028
1029			bool foundLower;
1030
1031			int iter;
1032
1033			int maxLSIter;
1034
1035			double mu;
1036
1037			double eta;
1038
1039			double ftol;
1040
1041			double lsTol;
1042
1043
1044
1045			xa.resize(ndim);
1046
1047			xb.resize(ndim);
1048
1049			xc.resize(ndim);
1050
1051
1052
1053			ga.resize(ndim);
1054
1055			gb.resize(ndim);
1056
1057			gc.resize(ndim);
1058
1059
1060
1061			a = 0.0;
1062
1063			fa = curF;
1064
1065			xa = curX;
1066
1067			ga = curG;
1068
1069			c = a + stepSize;
1070
1071			ftol = paramSet->getFTol();
1072
1073			lsTol = paramSet->getLineSearchTol();
1074
1075
1076
1077			//calculate the derivative at a = 0
1078
1079			slopeA = 0;
1080
1081			for (size_t i = 0; i < ndim; i++)
1082
1083			slopeA += curG[i]*direction[i];
1084
1085
1086
1087			initSlope = slopeA;
1088
1089
1090
1091			// if going uphill, use negative gradient as searching direction
1092
1093			if (slopeA > 0) {
1094
1095
1096
1097			if (bVerbose){
1098
1099			cout << "LineSearch Warning: initial searching direction is not a descent searching direction, "
1100
1101			<< " use negative gradient instead. Therefore, finding a smaller vaule of function "
1102
1103			<< " is guaranteed"
1104
1105	tim	1830	<< std::endl;
1106	tim	1826
1107			}
1108
1109
1110
1111			for (size_t i = 0; i < ndim; i++)
1112
1113			direction[i] = -curG[i];
1114
1115
1116
1117			for (size_t i = 0; i < ndim; i++)
1118
1119			slopeA += curG[i]*direction[i];
1120
1121
1122
1123			initSlope = slopeA;
1124
1125			}
1126
1127
1128
1129			// Take a trial step
1130
1131			for(size_t i = 0; i < ndim; i++)
1132
1133			xc[i] = curX[i] + direction[i] * c;
1134
1135
1136
1137			calcG(xc, gc, fc, status);
1138
1139
1140
1141			if (status < 0){
1142
1143			if (bVerbose)
1144
1145	tim	1830	std::cerr << "Function Evaluation Error" << std::endl;
1146	tim	1826
1147			}
1148
1149
1150
1151			//calculate the derivative at c
1152
1153			slopeC = 0;
1154
1155			for (size_t i = 0; i < ndim; i++)
1156
1157			slopeC += gc[i]*direction[i];
1158
1159
1160
1161			// found a lower point
1162
1163			if (fc < fa) {
1164
1165			curX = xc;
1166
1167			curG = gc;
1168
1169			curF = fc;
1170
1171			return LS_SUCCEED;
1172
1173			}
1174
1175			else {
1176
1177
1178
1179			if (slopeC > 0)
1180
1181			stepSize *= 0.618034;
1182
1183			}
1184
1185
1186
1187			maxLSIter = paramSet->getLineSearchMaxIteration();
1188
1189
1190
1191			iter = 0;
1192
1193
1194
1195			do {
1196
1197			// Select a new trial point.
1198
1199			// If the derivatives at points a & c have different sign we use cubic interpolate
1200
1201			//if (slopeC > 0){
1202
1203			eta = 3 *(fa -fc) /(c - a) + slopeA + slopeC;
1204
1205			mu = sqrt(eta * eta - slopeA * slopeC);
1206
1207			b = a + (c - a) * (1 - (slopeC + mu - eta) /(slopeC - slopeA + 2 * mu));
1208
1209
1210
1211			if (b < lsTol){
1212
1213			if (bVerbose)
1214
1215	tim	1830	cout << "stepSize is less than line search tolerance" << std::endl;
1216	tim	1826
1217			break;
1218
1219			}
1220
1221			//}
1222
1223
1224
1225			// Take a trial step to this new point - new coords in xb
1226
1227			for(size_t i = 0; i < ndim; i++)
1228
1229			xb[i] = curX[i] + direction[i] * b;
1230
1231
1232
1233			//function evaluation
1234
1235			calcG(xb, gb, fb, status);
1236
1237
1238
1239			if (status < 0){
1240
1241			if (bVerbose)
1242
1243	tim	1830	std::cerr << "Function Evaluation Error" << std::endl;
1244	tim	1826
1245			}
1246
1247
1248
1249			//calculate the derivative at c
1250
1251			slopeB = 0;
1252
1253			for (size_t i = 0; i < ndim; i++)
1254
1255			slopeB += gb[i]*direction[i];
1256
1257
1258
1259			//Amijo Rule to stop the line search
1260
1261			if (fb <= curF + initSlope * ftol * b) {
1262
1263			curF = fb;
1264
1265			curX = xb;
1266
1267			curG = gb;
1268
1269			return LS_SUCCEED;
1270
1271			}
1272
1273
1274
1275			if (slopeB <0 && fb < fa) {
1276
1277			//replace a by b
1278
1279			fa = fb;
1280
1281			a = b;
1282
1283			slopeA = slopeB;
1284
1285
1286
1287			// swap coord a/b
1288
1289			std::swap(xa, xb);
1290
1291			std::swap(ga, gb);
1292
1293			}
1294
1295			else {
1296
1297			//replace c by b
1298
1299			fc = fb;
1300
1301			c = b;
1302
1303			slopeC = slopeB;
1304
1305
1306
1307			// swap coord b/c
1308
1309			std::swap(gb, gc);
1310
1311			std::swap(xb, xc);
1312
1313			}
1314
1315
1316
1317
1318
1319			iter++;
1320
1321			} while((fb > fa \|\| fb > fc) && (iter < maxLSIter));
1322
1323
1324
1325			if(fb < curF \|\| iter >= maxLSIter) {
1326
1327			//could not find a lower value, we might just go uphill.
1328
1329			return LS_ERROR;
1330
1331			}
1332
1333
1334
1335			//select the end point
1336
1337
1338
1339			if (fa <= fc) {
1340
1341			curX = xa;
1342
1343			curG = ga;
1344
1345			curF = fa;
1346
1347			}
1348
1349			else {
1350
1351			curX = xc;
1352
1353			curG = gc;
1354
1355			curF = fc;
1356
1357			}
1358
1359
1360
1361			return LS_SUCCEED;
1362
1363
1364
1365			}
1366
1367
1368
1369			void OOPSEMinimizer::minimize(){
1370
1371
1372
1373			int convgStatus;
1374
1375			int stepStatus;
1376
1377			int maxIter;
1378
1379			int writeFrq;
1380
1381			int nextWriteIter;
1382
1383
1384
1385			if (bVerbose)
1386
1387			printMinimizerInfo();
1388
1389
1390
1391			dumpOut = new DumpWriter(info);
1392
1393			statOut = new StatWriter(info);
1394
1395
1396
1397			init();
1398
1399
1400
1401			writeFrq = paramSet->getWriteFrq();
1402
1403			nextWriteIter = writeFrq;
1404
1405
1406
1407			maxIter = paramSet->getMaxIteration();
1408
1409
1410
1411			for (curIter = 1; curIter <= maxIter; curIter++){
1412
1413
1414
1415			stepStatus = step();
1416
1417
1418
1419			if (nConstrained && bShake)
1420
1421			preMove();
1422
1423
1424
1425			if (stepStatus < 0){
1426
1427			saveResult();
1428
1429			minStatus = MIN_LSERROR;
1430
1431	tim	1830	std::cerr << "OOPSEMinimizer Error: line search error, please try a small stepsize" << std::endl;
1432	tim	1826
1433			return;
1434
1435			}
1436
1437
1438
1439			if (curIter == nextWriteIter){
1440
1441			nextWriteIter += writeFrq;
1442
1443			writeOut(curX, curIter);
1444
1445			}
1446
1447
1448
1449			convgStatus = checkConvg();
1450
1451
1452
1453			if (convgStatus > 0){
1454
1455			saveResult();
1456
1457			minStatus = MIN_CONVERGE;
1458
1459			return;
1460
1461			}
1462
1463
1464
1465			prepareStep();
1466
1467
1468
1469			}
1470
1471
1472
1473
1474
1475			if (bVerbose) {
1476
1477			cout << "OOPSEMinimizer Warning: "
1478
1479			<< minimizerName << " algorithm did not converge within "
1480
1481	tim	1830	<< maxIter << " iteration" << std::endl;
1482	tim	1826
1483			}
1484
1485			minStatus = MIN_MAXITER;
1486
1487			saveResult();
1488
1489
1490
1491			}
1492