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 1234 by tim, Fri Jun 4 03:15:31 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 +  //preMove();
22   }
23  
24   OOPSEMinimizer::~OOPSEMinimizer(){
25    delete tStats;
26 <  delete dumpOut;
27 <  delete statOut;
26 >  if(dumpOut)
27 >    delete dumpOut;
28 >  if(statOut)
29 >    delete statOut;
30    delete paramSet;
31   }
32  
# Line 38 | Line 41 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
41    
42    setCoor(x);
43  
44 <  if (nConstrained && bShake){
45 <    shakeStatus = shakeR();
46 <  }
44 >  //if (nConstrained && bShake){
45 >  //  shakeStatus = shakeR();
46 >  //}
47 >
48 >  shakeAlgo->doShakeR();
49        
50    calcForce(1, 1);
51    
52 <  if (nConstrained && bShake){
53 <    shakeStatus |= shakeF();
54 <  }
52 >  //if (nConstrained && bShake){
53 >  //  shakeStatus |= shakeF();
54 >  //}
55 >
56 >  shakeAlgo->doShakeF();
57    
58    x = getCoor();
59    
60  
61    index = 0;
62  
63 <  for(int i = 0; i < nAtoms; i++){
63 >  for(int i = 0; i < integrableObjects.size(); i++){
64  
65 <    if(atoms[i]->isDirectional()){
59 <      dAtom = (DirectionalAtom*) atoms[i];
60 <      dAtom->getGrad(dAtomGrad);
65 >    if (integrableObjects[i]->isDirectional()) {
66  
67 +      integrableObjects[i]->getGrad(dAtomGrad);
68 +
69        //gradient is equal to -f
70        grad[index++] = -dAtomGrad[0];
71        grad[index++] = -dAtomGrad[1];
# Line 69 | Line 76 | void OOPSEMinimizer::calcEnergyGradient(vector<double>
76  
77      }
78      else{
79 <      atoms[i]->getFrc(force);
79 >      integrableObjects[i]->getFrc(force);
80  
81        grad[index++] = -force[0];
82        grad[index++] = -force[1];
# Line 98 | Line 105 | void OOPSEMinimizer::setCoor(vector<double>& x){
105  
106    index = 0;
107    
108 <  for(int i = 0; i < nAtoms; i++){
108 >  for(int i = 0; i < integrableObjects.size(); i++){
109      
110      position[0] = x[index++];
111      position[1] = x[index++];
112      position[2] = x[index++];
113  
114 <    atoms[i]->setPos(position);
114 >    integrableObjects[i]->setPos(position);
115  
116 <    if (atoms[i]->isDirectional()){
110 <      dAtom = (DirectionalAtom*) atoms[i];
116 >    if (integrableObjects[i]->isDirectional()){
117  
118        eulerAngle[0] = x[index++];
119        eulerAngle[1] = x[index++];
120        eulerAngle[2] = x[index++];
121  
122 <       dAtom->setEuler(eulerAngle[0], eulerAngle[1], eulerAngle[2]);
122 >      integrableObjects[i]->setEuler(eulerAngle[0],
123 >                                     eulerAngle[1],
124 >                                     eulerAngle[2]);
125  
126      }
127      
# Line 136 | Line 144 | vector<double> OOPSEMinimizer::getCoor(){
144  
145    index = 0;
146    
147 <  for(int i = 0; i < nAtoms; i++){
148 <    atoms[i]->getPos(position);
147 >  for(int i = 0; i < integrableObjects.size(); i++){
148 >    integrableObjects[i]->getPos(position);
149  
150      x[index++] = position[0];
151      x[index++] = position[1];
152      x[index++] = position[2];
153  
154 <    if (atoms[i]->isDirectional()){
147 <      dAtom = (DirectionalAtom*) atoms[i];
148 <      dAtom->getEulerAngles(eulerAngle);
154 >    if (integrableObjects[i]->isDirectional()){
155  
156 +      integrableObjects[i]->getEulerAngles(eulerAngle);
157 +      
158        x[index++] = eulerAngle[0];
159        x[index++] = eulerAngle[1];
160        x[index++] = eulerAngle[2];
# Line 159 | Line 167 | vector<double> OOPSEMinimizer::getCoor(){
167  
168   }
169  
170 + /*
171   int OOPSEMinimizer::shakeR(){
172    int i, j;
173    int done;
# Line 398 | Line 407 | int OOPSEMinimizer::shakeF(){
407      return 1;
408   }
409  
410 <    //calculate the value of object function
410 > */
411 >
412 > //calculate the value of object function
413   void OOPSEMinimizer::calcF(){
414    calcEnergyGradient(curX, curG, curF, egEvalStatus);
415   }
# Line 424 | Line 435 | void OOPSEMinimizer::calcDim(){
435  
436    ndim = 0;
437  
438 <  for(int i = 0; i < nAtoms; i++){
438 >  for(int i = 0; i < integrableObjects.size(); i++){
439      ndim += 3;
440 <    if (atoms[i]->isDirectional())
440 >    if (integrableObjects[i]->isDirectional())
441        ndim += 3;      
442    }
443   }
# Line 483 | Line 494 | void OOPSEMinimizer::printMinimizerInfo(){
494  
495   /**
496   * In thoery, we need to find the minimum along the search direction
497 < * However, function evaluation is too expensive. I
497 > * However, function evaluation is too expensive.
498   * At the very begining of the problem, we check the search direction and make sure
499   * it is a descent direction
500   * we will compare the energy of two end points,
# Line 517 | Line 528 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
528    double mu;
529    double eta;
530    double ftol;  
531 +  double lsTol;
532  
533    xa.resize(ndim);
534    xb.resize(ndim);
# Line 532 | Line 544 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
544    ga = curG;
545    c = a + stepSize;
546    ftol = paramSet->getFTol();
547 +  lsTol = paramSet->getLineSearchTol();
548          
549    //calculate the derivative at a = 0
550    for (size_t i = 0; i < ndim; i++)
# Line 598 | Line 611 | int OOPSEMinimizer::doLineSearch(vector<double>& direc
611        eta = 3 *(fa -fc) /(c - a) + slopeA + slopeC;
612        mu = sqrt(eta * eta - slopeA * slopeC);      
613        b = a + (c - a) * (1 - (slopeC + mu - eta) /(slopeC - slopeA + 2 * mu));      
614 +
615 +      if (b < lsTol){
616 +        if (bVerbose)
617 +          cout << "stepSize is less than line search tolerance" << endl;
618 +        break;        
619 +      }
620      //}
621  
622      // Take a trial step to this new point - new coords in xb
# Line 684 | Line 703 | void OOPSEMinimizer::minimize(){
703    
704    if (bVerbose)
705      printMinimizerInfo();
706 +
707 +  dumpOut = new DumpWriter(info);
708 +  statOut = new StatWriter(info);
709    
710    init();
711  

Diff Legend

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