ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/minimizers/Minimizer.cpp
(Generate patch)

Comparing branches/development/src/minimizers/Minimizer.cpp (file contents):
Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC vs.
Revision 1744 by gezelter, Tue Jun 5 18:07:08 2012 UTC

# Line 50 | Line 50 | namespace OpenMD {
50   #include <mpi.h>
51   #endif
52   namespace OpenMD {
53 –  RealType dotProduct(const std::vector<RealType>& v1, const std::vector<RealType>& v2) {
54 –    if (v1.size() != v2.size()) {
53  
56 –    }
57 –
58 –
59 –    RealType result = 0.0;    
60 –    for (unsigned int i = 0; i < v1.size(); ++i) {
61 –      result += v1[i] * v2[i];
62 –    }
63 –
64 –    return result;
65 –  }
66 –
54    Minimizer::Minimizer(SimInfo* rhs) :
55      info(rhs), usingShake(false) {
56 <
57 <      forceMan = new ForceManager(info);
58 <      paramSet= new MinimizerParameterSet(info), calcDim();
59 <      curX = getCoor();
60 <      curG.resize(ndim);
61 <
62 <    }
76 <
56 >    
57 >    forceMan = new ForceManager(info);
58 >    paramSet= new MinimizerParameterSet(info), calcDim();
59 >    curX = getCoor();
60 >    curG.resize(ndim);
61 >  }
62 >  
63    Minimizer::~Minimizer() {
64      delete forceMan;
65      delete paramSet;
66    }
67 +  
68 +  // void Minimizer::calcEnergyGradient(std::vector<RealType> &x,
69 +  //                                 std::vector<RealType> &grad,
70 +  //                                    RealType&energy, int&status) {
71  
72 <  void Minimizer::calcEnergyGradient(std::vector<RealType> &x,
73 <                                     std::vector<RealType> &grad, RealType&energy, int&status) {
72 >  //   SimInfo::MoleculeIterator i;
73 >  //   Molecule::IntegrableObjectIterator  j;
74 >  //   Molecule* mol;
75 >  //   StuntDouble* integrableObject;    
76 >  //   std::vector<RealType> myGrad;    
77 >  //   int shakeStatus;
78  
79 <    SimInfo::MoleculeIterator i;
86 <    Molecule::IntegrableObjectIterator  j;
87 <    Molecule* mol;
88 <    StuntDouble* integrableObject;    
89 <    std::vector<RealType> myGrad;    
90 <    int shakeStatus;
79 >  //   status = 1;
80  
81 <    status = 1;
81 >  //   setCoor(x);
82  
83 <    setCoor(x);
83 >  //   if (usingShake) {
84 >  //     shakeStatus = shakeR();
85 >  //   }
86  
87 <    if (usingShake) {
97 <      shakeStatus = shakeR();
98 <    }
87 >  //   energy = calcPotential();
88  
89 <    energy = calcPotential();
89 >  //   if (usingShake) {
90 >  //     shakeStatus = shakeF();
91 >  //   }
92  
93 <    if (usingShake) {
103 <      shakeStatus = shakeF();
104 <    }
93 >  //   x = getCoordinates();
94  
95 <    x = getCoor();
95 >  //   int index = 0;
96  
97 <    int index = 0;
97 >  //   for (mol = info->beginMolecule(i); mol != NULL;
98 >  //        mol = info->nextMolecule(i)) {
99 >  //     for (integrableObject = mol->beginIntegrableObject(j);
100 >  //          integrableObject != NULL;
101 >  //          integrableObject = mol->nextIntegrableObject(j)) {        
102 >  //       myGrad = integrableObject->getGrad();
103 >  //       for (unsigned int k = 0; k < myGrad.size(); ++k) {
104 >  //         grad[index++] = myGrad[k];
105 >  //       }
106 >  //     }            
107 >  //   }    
108 >  // }
109  
110 –    for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
111 –      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
112 –           integrableObject = mol->nextIntegrableObject(j)) {
113 –
114 –        myGrad = integrableObject->getGrad();
115 –        for (unsigned int k = 0; k < myGrad.size(); ++k) {
116 –
117 –          grad[index++] = myGrad[k];
118 –        }
119 –      }            
120 –    }
121 –
122 –  }
123 –
110    void Minimizer::setCoor(std::vector<RealType> &x) {
111      Vector3d position;
112      Vector3d eulerAngle;
# Line 130 | Line 116 | namespace OpenMD {
116      StuntDouble* integrableObject;    
117      int index = 0;
118  
119 <    for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
120 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
119 >    for (mol = info->beginMolecule(i); mol != NULL;
120 >         mol = info->nextMolecule(i)) {
121 >      for (integrableObject = mol->beginIntegrableObject(j);
122 >           integrableObject != NULL;
123             integrableObject = mol->nextIntegrableObject(j)) {
124 <
124 >        
125          position[0] = x[index++];
126          position[1] = x[index++];
127          position[2] = x[index++];
128 <
128 >        
129          integrableObject->setPos(position);
130 <
130 >        
131          if (integrableObject->isDirectional()) {
132            eulerAngle[0] = x[index++];
133            eulerAngle[1] = x[index++];
134            eulerAngle[2] = x[index++];
135 <
135 >          
136            integrableObject->setEuler(eulerAngle);
137          }
138        }
139 <    }
152 <    
139 >    }    
140    }
141  
142    std::vector<RealType> Minimizer::getCoor() {
# Line 162 | Line 149 | namespace OpenMD {
149      int index = 0;
150      std::vector<RealType> x(getDim());
151  
152 <    for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
153 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
152 >    for (mol = info->beginMolecule(i); mol != NULL;
153 >         mol = info->nextMolecule(i)) {
154 >      for (integrableObject = mol->beginIntegrableObject(j);
155 >           integrableObject != NULL;
156             integrableObject = mol->nextIntegrableObject(j)) {
157 <                
157 >        
158          position = integrableObject->getPos();
159          x[index++] = position[0];
160          x[index++] = position[1];
161          x[index++] = position[2];
162 <
162 >        
163          if (integrableObject->isDirectional()) {
164            eulerAngle = integrableObject->getEuler();
165            x[index++] = eulerAngle[0];
# Line 181 | Line 170 | namespace OpenMD {
170      }
171      return x;
172    }
173 <
174 <
173 >  
174 >  
175    /*
176      int Minimizer::shakeR() {
177      int    i,       j;
178 <
178 >    
179      int    done;
180  
181      RealType posA[3], posB[3];
# Line 524 | Line 513 | namespace OpenMD {
513    //calculate the value of object function
514  
515    void Minimizer::calcF() {
516 <    calcEnergyGradient(curX, curG, curF, egEvalStatus);
516 >    egEvalStatus = calcEnergyGradient(curX, curG, curF);
517    }
518 <
518 >  
519    void Minimizer::calcF(std::vector < RealType > &x, RealType&f, int&status) {
520      std::vector < RealType > tempG;
521 <
521 >    
522      tempG.resize(x.size());
523 <
524 <    calcEnergyGradient(x, tempG, f, status);
523 >    
524 >    status = calcEnergyGradient(x, tempG, f);
525    }
526 <
526 >  
527    //calculate the gradient
528 <
528 >  
529    void Minimizer::calcG() {
530 <    calcEnergyGradient(curX, curG, curF, egEvalStatus);
530 >    egEvalStatus = calcEnergyGradient(curX, curG, curF);
531    }
532 <
533 <  void Minimizer::calcG(std::vector<RealType>& x, std::vector<RealType>& g, RealType&f, int&status) {
534 <    calcEnergyGradient(x, g, f, status);
535 <  }
536 <
532 >  
533 >  void Minimizer::calcG(std::vector<RealType>& x,
534 >                        std::vector<RealType>& g, RealType&f, int&status) {
535 >    status = calcEnergyGradient(x, g, f);
536 >  }
537 >  
538    void Minimizer::calcDim() {
539 <
539 >    
540      SimInfo::MoleculeIterator i;
541      Molecule::IntegrableObjectIterator  j;
542      Molecule* mol;
543      StuntDouble* integrableObject;    
544      ndim = 0;
545 <
546 <    for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
547 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
545 >    
546 >    for (mol = info->beginMolecule(i); mol != NULL;
547 >         mol = info->nextMolecule(i)) {
548 >      for (integrableObject = mol->beginIntegrableObject(j);
549 >           integrableObject != NULL;
550             integrableObject = mol->nextIntegrableObject(j)) {
551 <
551 >        
552          ndim += 3;
553 <
553 >        
554          if (integrableObject->isDirectional()) {
555            ndim += 3;
556          }
557 <      }
566 <
557 >      }      
558      }
559    }
560  
561 <  void Minimizer::setX(std::vector < RealType > &x) {
561 >  void Minimizer::setX(std::vector<RealType> &x) {
562      if (x.size() != ndim) {
563 <      sprintf(painCave.errMsg, "Minimizer Error: dimesion of x and curX does not match\n");
563 >      sprintf(painCave.errMsg,
564 >              "Minimizer setX: dimensions of x and curX do not match\n");
565        painCave.isFatal = 1;
566        simError();
567      }
# Line 577 | Line 569 | namespace OpenMD {
569      curX = x;
570    }
571  
572 <  void Minimizer::setG(std::vector < RealType > &g) {
572 >  void Minimizer::setG(std::vector <RealType> &g) {
573      if (g.size() != ndim) {
574 <      sprintf(painCave.errMsg, "Minimizer Error: dimesion of g and curG does not match\n");
574 >      sprintf(painCave.errMsg,
575 >              "Minimizer setG: dimensions of g and curG do not match\n");
576        painCave.isFatal = 1;
577        simError();
578      }
579 <
579 >    
580      curG = g;
581    }
582  
583  
584    /**
585 <
586 <  * In thoery, we need to find the minimum along the search direction
587 <  * However, function evaluation is too expensive.
588 <  * At the very begining of the problem, we check the search direction and make sure
589 <  * it is a descent direction
590 <  * we will compare the energy of two end points,
598 <  * if the right end point has lower energy, we just take it
599 <  * @todo optimize this line search algorithm
585 >  * In theory, we need to find the minimum along the search direction
586 >  * However, function evaluation is usually too expensive.  At the
587 >  * very begining of the problem, we check the search direction and
588 >  * make sure it is a descent direction we will compare the energy of
589 >  * two end points, if the right end point has lower energy, we'll
590 >  * just take it.
591    */
592  
593    int Minimizer::doLineSearch(std::vector<RealType> &direction,
# Line 635 | Line 626 | namespace OpenMD {
626      gc.resize(ndim);
627  
628      a = 0.0;
638 –
629      fa = curF;
640 –
630      xa = curX;
642 –
631      ga = curG;
632  
633      c = a + stepSize;
634  
635      ftol = paramSet->getFTol();
648 –
636      lsTol = paramSet->getLineSearchTol();
637  
638      //calculate the derivative at a = 0
# Line 655 | Line 642 | namespace OpenMD {
642      for(size_t i = 0; i < ndim; i++) {
643        slopeA += curG[i] * direction[i];
644      }
645 <    
645 >
646 > #ifdef IS_MPI
647 >    // in parallel, we need to add up the contributions from all
648 >    // processors:
649 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &slopeA, 1, MPI::REALTYPE,
650 >                              MPI::SUM);
651 > #endif
652 >
653      initSlope = slopeA;
654  
655      // if  going uphill, use negative gradient as searching direction
# Line 670 | Line 664 | namespace OpenMD {
664          slopeA += curG[i] * direction[i];
665        }
666          
667 + #ifdef IS_MPI
668 +      // in parallel, we need to add up the contributions from all
669 +      // processors:
670 +      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &slopeA, 1, MPI::REALTYPE,
671 +                              MPI::SUM);
672 + #endif
673        initSlope = slopeA;
674      }
675  
# Line 697 | Line 697 | namespace OpenMD {
697  
698      if (fc < fa) {
699        curX = xc;
700 –
700        curG = gc;
702 –
701        curF = fc;
704 –
702        return LS_SUCCEED;
703      } else {
704        if (slopeC > 0)

Diff Legend

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