ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/OOPSEMinimizer.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/OOPSEMinimizer.cpp (file contents):
Revision 1064 by tim, Tue Feb 24 15:44:45 2004 UTC vs.
Revision 1330 by gezelter, Fri Jul 16 16:29:44 2004 UTC

# Line 1 | Line 1
1   #include <math.h>
2   #include "OOPSEMinimizer.hpp"
3 + #include "ShakeMin.hpp"
4 + #include "Integrator.cpp"
5  
6   OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, ForceFields* the_ff ,
7 <                                              MinimizerParameterSet * param)
8 <                     :RealIntegrator(theInfo, the_ff), bVerbose(false), bShake(true){
7 >                                              MinimizerParameterSet * param) :
8 >              RealIntegrator(theInfo, the_ff), bShake(true), bVerbose(false) {
9 >  dumpOut = NULL;
10 >  statOut = NULL;
11  
8 –  atoms = info->atoms;
9 –
12    tStats = new Thermo(info);
13 <  dumpOut = new DumpWriter(info);
12 <  statOut = new StatWriter(info);
13 >
14    
15    paramSet = param;
16  
# Line 18 | Line 19 | OOPSEMinimizer::OOPSEMinimizer( SimInfo *theInfo, Forc
19    curX = getCoor();
20    curG.resize(ndim);
21  
22 +  shakeAlgo = new ShakeMinFramework(theInfo);
23 +  shakeAlgo->doPreConstraint();
24   }
25  
26   OOPSEMinimizer::~OOPSEMinimizer(){
27    delete tStats;
28 <  delete dumpOut;
29 <  delete statOut;
28 >  if(dumpOut)
29 >    delete dumpOut;
30 >  if(statOut)
31 >    delete statOut;
32    delete paramSet;
33   }
34  
# Line 35 | Line 40 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
40    double force[3];
41    double dAtomGrad[6];
42    int shakeStatus;
43 +
44 +  status = 1;
45    
46    setCoor(x);
47  
48 <  if (nConstrained && bShake){
49 <    shakeStatus = shakeR();
48 >  if (bShake){
49 >    shakeStatus = shakeAlgo->doShakeR();
50 >    if(shakeStatus < 0)
51 >      status = -1;
52    }
53 <      
53 >
54    calcForce(1, 1);
55 <  
56 <  if (nConstrained && bShake){
57 <    shakeStatus |= shakeF();
55 >
56 >  if (bShake){
57 >    shakeStatus = shakeAlgo->doShakeF();
58 >    if(shakeStatus < 0)
59 >      status = -1;
60    }
61 <  
61 >
62    x = getCoor();
63    
64  
65    index = 0;
66  
67 <  for(int i = 0; i < nAtoms; i++){
67 >  for(int i = 0; i < integrableObjects.size(); i++){
68  
69 <    if(atoms[i]->isDirectional()){
59 <      dAtom = (DirectionalAtom*) atoms[i];
60 <      dAtom->getGrad(dAtomGrad);
69 >    if (integrableObjects[i]->isDirectional()) {
70  
71 +      integrableObjects[i]->getGrad(dAtomGrad);
72 +
73        //gradient is equal to -f
74        grad[index++] = -dAtomGrad[0];
75        grad[index++] = -dAtomGrad[1];
# Line 69 | Line 80 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
80  
81      }
82      else{
83 <      atoms[i]->getFrc(force);
83 >      integrableObjects[i]->getFrc(force);
84  
85        grad[index++] = -force[0];
86        grad[index++] = -force[1];
# Line 81 | Line 92 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
92  
93    energy = tStats->getPotential();
94  
84 –  status = shakeStatus;
95   }
96  
97   /**
# Line 98 | Line 108 | void OOPSEMinimizer::setCoor(vector<double>& x){
108  
109    index = 0;
110    
111 <  for(int i = 0; i < nAtoms; i++){
111 >  for(int i = 0; i < integrableObjects.size(); i++){
112      
113      position[0] = x[index++];
114      position[1] = x[index++];
115      position[2] = x[index++];
116  
117 <    atoms[i]->setPos(position);
117 >    integrableObjects[i]->setPos(position);
118  
119 <    if (atoms[i]->isDirectional()){
110 <      dAtom = (DirectionalAtom*) atoms[i];
119 >    if (integrableObjects[i]->isDirectional()){
120  
121        eulerAngle[0] = x[index++];
122        eulerAngle[1] = x[index++];
123        eulerAngle[2] = x[index++];
124  
125 <       dAtom->setEuler(eulerAngle[0], eulerAngle[1], eulerAngle[2]);
125 >      integrableObjects[i]->setEuler(eulerAngle[0],
126 >                                     eulerAngle[1],
127 >                                     eulerAngle[2]);
128  
129      }
130      
# Line 136 | Line 147 | vector<double> OOPSEMinimizer::getCoor(){
147  
148    index = 0;
149    
150 <  for(int i = 0; i < nAtoms; i++){
151 <    atoms[i]->getPos(position);
150 >  for(int i = 0; i < integrableObjects.size(); i++){
151 >    integrableObjects[i]->getPos(position);
152  
153      x[index++] = position[0];
154      x[index++] = position[1];
155      x[index++] = position[2];
156  
157 <    if (atoms[i]->isDirectional()){
147 <      dAtom = (DirectionalAtom*) atoms[i];
148 <      dAtom->getEulerAngles(eulerAngle);
157 >    if (integrableObjects[i]->isDirectional()){
158  
159 +      integrableObjects[i]->getEulerAngles(eulerAngle);
160 +      
161        x[index++] = eulerAngle[0];
162        x[index++] = eulerAngle[1];
163        x[index++] = eulerAngle[2];
# Line 159 | Line 170 | vector<double> OOPSEMinimizer::getCoor(){
170  
171   }
172  
173 + /*
174   int OOPSEMinimizer::shakeR(){
175    int i, j;
176    int done;
# Line 398 | Line 410 | int OOPSEMinimizer::shakeF(){
410      return 1;
411   }
412  
413 <    //calculate the value of object function
413 > */
414 >
415 > //calculate the value of object function
416   void OOPSEMinimizer::calcF(){
417    calcEnergyGradient(curX, curG, curF, egEvalStatus);
418   }
# Line 424 | Line 438 | void OOPSEMinimizer::calcDim(){
438  
439    ndim = 0;
440  
441 <  for(int i = 0; i < nAtoms; i++){
441 >  for(int i = 0; i < integrableObjects.size(); i++){
442      ndim += 3;
443 <    if (atoms[i]->isDirectional())
443 >    if (integrableObjects[i]->isDirectional())
444        ndim += 3;      
445    }
446   }
# Line 483 | Line 497 | void OOPSEMinimizer::printMinimizerInfo(){
497  
498   /**
499   * In thoery, we need to find the minimum along the search direction
500 < * However, function evaluation is too expensive. I
500 > * However, function evaluation is too expensive.
501   * 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,
# Line 517 | Line 531 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
531    double mu;
532    double eta;
533    double ftol;  
534 +  double lsTol;
535  
536    xa.resize(ndim);
537    xb.resize(ndim);
# Line 532 | Line 547 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
547    ga = curG;
548    c = a + stepSize;
549    ftol = paramSet->getFTol();
550 +  lsTol = paramSet->getLineSearchTol();
551          
552    //calculate the derivative at a = 0
553 +  slopeA = 0;
554    for (size_t i = 0; i < ndim; i++)
555      slopeA += curG[i]*direction[i];
556  
# Line 598 | Line 615 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
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 +
619 +      if (b < lsTol){
620 +        if (bVerbose)
621 +          cout << "stepSize is less than line search tolerance" << endl;
622 +        break;        
623 +      }
624      //}
625  
626      // Take a trial step to this new point - new coords in xb
# Line 684 | Line 707 | void OOPSEMinimizer::minimize(){
707    
708    if (bVerbose)
709      printMinimizerInfo();
710 +
711 +  dumpOut = new DumpWriter(info);
712 +  statOut = new StatWriter(info);
713    
714    init();
715  
# Line 699 | Line 725 | void OOPSEMinimizer::minimize(){
725  
726      stepStatus = step();
727  
728 +    if (bShake)
729 +      shakeAlgo->doPreConstraint();
730 +    
731      if (stepStatus < 0){
732        saveResult();
733        minStatus = MIN_LSERROR;
# Line 711 | Line 740 | void OOPSEMinimizer::minimize(){
740        writeOut(curX, curIter);
741      }
742  
714 –    //if (curIter == nextResetIter){
715 –    //  reset();
716 –    //  nextResetIter += resetFrq;      
717 –    //}
718 –
743      convgStatus = checkConvg();
744  
745      if (convgStatus > 0){

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines