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 1250 by gezelter, Fri Jun 4 21:00:20 2004 UTC

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

Diff Legend

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