ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/minimizers/OOPSEMinimizer.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-4/src/minimizers/OOPSEMinimizer.cpp (file contents):
Revision 1899 by tim, Tue Dec 14 19:08:44 2004 UTC vs.
Revision 1900 by tim, Tue Jan 4 22:18:06 2005 UTC

# Line 1 | Line 1
1   #include <cmath>
2  
3 – #include "minimizers/OOPSEMinimizer.hpp"
4 – #include "integrators/Integrator.cpp"
3  
4 + #include "integrators/Integrator.cpp"
5 + #include "io/StatWriter.hpp"
6 + #include "minimizers/OOPSEMinimizer.hpp"
7 + #include "primitives/Molecule.hpp"
8   namespace oopse {
9  
10 < OOPSEMinimizer::OOPSEMinimizer(SimInfo* theInfo, MinimizerParameterSet *param) :
11 <    info(theInfo), paramSet(param), bShake(true), bVerbose(false) {
10 > OOPSEMinimizer::OOPSEMinimizer(SimInfo* rhs) :
11 >    info(rhs), usingShake(false) {
12  
13 +    forceMan = new ForceManager(info);
14 +    paramSet= new MinimizerParameterSet(info),
15      calcDim();
16      curX = getCoor();
17      curG.resize(ndim);
# Line 15 | Line 19 | OOPSEMinimizer::~OOPSEMinimizer() {
19   }
20  
21   OOPSEMinimizer::~OOPSEMinimizer() {
22 +    delete forceMan;
23      delete paramSet;
24   }
25  
26 < void OOPSEMinimizer::calcEnergyGradient(std::vector < double > &x,
27 <    std::vector < double > &grad, double&energy, int&status) {
26 > void OOPSEMinimizer::calcEnergyGradient(std::vector<double> &x,
27 >    std::vector<double> &grad, double&energy, int&status) {
28  
29 <    DirectionalAtom *dAtom;
30 <    int    index;
31 <    double force[3];
32 <    double dAtomGrad[6];
29 >    SimInfo::MoleculeIterator i;
30 >    Molecule::IntegrableObjectIterator  j;
31 >    Molecule* mol;
32 >    StuntDouble* integrableObject;    
33 >    std::vector<double> myGrad;    
34 >    int shakeStatus;
35  
29 –    int    shakeStatus;
30 –
36      status = 1;
37  
38      setCoor(x);
39  
40 <    if (nConstrained && bShake) {
40 >    if (usingShake) {
41          shakeStatus = shakeR();
42      }
43  
44 <    calcForce(1, 1);
44 >    energy = calcPotential();
45  
46 <    if (nConstrained && bShake) {
46 >    if (usingShake) {
47          shakeStatus = shakeF();
48      }
49  
50      x = getCoor();
51  
52 <    index = 0;
52 >    int index = 0;
53  
54 <    for(int i = 0; i < integrableObjects.size(); i++) {
55 <        if (integrableObjects[i]->isDirectional()) {
56 <            integrableObjects[i]->getGrad(dAtomGrad);
52 <
53 <            //gradient is equal to -f
54 <
55 <            grad[index++] = -dAtomGrad[0];
56 <
57 <            grad[index++] = -dAtomGrad[1];
58 <
59 <            grad[index++] = -dAtomGrad[2];
54 >    for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
55 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
56 >               integrableObject = mol->nextIntegrableObject(j)) {
57  
58 <            grad[index++] = -dAtomGrad[3];
59 <
60 <            grad[index++] = -dAtomGrad[4];
61 <
62 <            grad[index++] = -dAtomGrad[5];
63 <        } else {
67 <            integrableObjects[i]->getFrc(force);
68 <
69 <            grad[index++] = -force[0];
70 <
71 <            grad[index++] = -force[1];
72 <
73 <            grad[index++] = -force[2];
74 <        }
58 >            myGrad = integrableObject->getGrad();
59 >            for (unsigned int k = 0; k < myGrad.size(); ++k) {
60 >                //gradient is equal to -f
61 >                grad[index++] = -myGrad[k];
62 >            }
63 >        }            
64      }
65  
77 –    energy = tStats->getPotential();
66   }
67  
68   void OOPSEMinimizer::setCoor(std::vector<double> &x) {
# Line 91 | Line 79 | void OOPSEMinimizer::setCoor(std::vector<double> &x) {
79                 integrableObject = mol->nextIntegrableObject(j)) {
80  
81              position[0] = x[index++];
94 –
82              position[1] = x[index++];
96 –
83              position[2] = x[index++];
84  
85              integrableObject->setPos(position);
86  
87              if (integrableObject->isDirectional()) {
88                  eulerAngle[0] = x[index++];
103 –
89                  eulerAngle[1] = x[index++];
105 –
90                  eulerAngle[2] = x[index++];
91  
92                  integrableObject->setEuler(eulerAngle);
# Line 517 | Line 501 | void OOPSEMinimizer::calcDim() {
501          for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
502                 integrableObject = mol->nextIntegrableObject(j)) {
503  
504 <        ndim += 3;
504 >            ndim += 3;
505  
506 <        if (integrableObject->isDirectional()) {
507 <            ndim += 3;
506 >            if (integrableObject->isDirectional()) {
507 >                ndim += 3;
508 >            }
509          }
525 –    }
510  
511 <
511 >    }
512   }
513  
514   void OOPSEMinimizer::setX(std::vector < double > &x) {
# Line 546 | Line 530 | void OOPSEMinimizer::setG(std::vector < double > &g) {
530  
531      curG = g;
532   }
549 –
550 – void OOPSEMinimizer::writeOut(std::vector < double > &x, double iter) {
551 –    setX(x);
552 –
553 –    calcG();
554 –
555 –    dumpOut->writeDump(iter);
556 –
557 –    statOut->writeStat(iter);
558 – }
559 –
560 – void OOPSEMinimizer::printMinimizerInfo() {
561 –    cout
562 –        << "--------------------------------------------------------------------"
563 –        << std::endl;
533  
565 –    cout << minimizerName << std::endl;
566 –
567 –    cout << "minimization parameter set" << std::endl;
568 –
569 –    cout << "function tolerance = " << paramSet->getFTol() << std::endl;
570 –
571 –    cout << "gradient tolerance = " << paramSet->getGTol() << std::endl;
572 –
573 –    cout << "step tolerance = " << paramSet->getFTol() << std::endl;
574 –
575 –    cout << "absolute gradient tolerance = " << std::endl;
576 –
577 –    cout << "max iteration = " << paramSet->getMaxIteration() << std::endl;
578 –
579 –    cout << "max line search iteration = "
580 –        << paramSet->getLineSearchMaxIteration() << std::endl;
581 –
582 –    cout << "shake algorithm = " << bShake << std::endl;
583 –
584 –    cout
585 –        << "--------------------------------------------------------------------"
586 –        << std::endl;
587 – }
534  
535   /**
536  
# Line 594 | Line 540 | void OOPSEMinimizer::printMinimizerInfo() {
540   * it is a descent direction
541   * we will compare the energy of two end points,
542   * if the right end point has lower energy, we just take it
543 + * @todo optimize this line search algorithm
544   */
545  
546   int OOPSEMinimizer::doLineSearch(std::vector<double> &direction,
# Line 658 | Line 605 | int OOPSEMinimizer::doLineSearch(std::vector<double> &
605      // if  going uphill, use negative gradient as searching direction
606  
607      if (slopeA > 0) {
661 –        if (bVerbose) {
662 –            cout
663 –                << "LineSearch Warning: initial searching direction is not a descent searching direction, " << " use negative gradient instead. Therefore, finding a smaller vaule of function " << " is guaranteed"
664 –                << std::endl;
665 –        }
608  
609          for(size_t i = 0; i < ndim; i++) {
610              direction[i] = -curG[i];
# Line 730 | Line 672 | int OOPSEMinimizer::doLineSearch(std::vector<double> &
672                  * (1 - (slopeC + mu - eta) / (slopeC - slopeA + 2 * mu));
673  
674          if (b < lsTol) {
733 –            if (bVerbose)
734 –                cout << "stepSize is less than line search tolerance"
735 –                    << std::endl;
736 –
675              break;
676          }
677  
# Line 837 | Line 775 | void OOPSEMinimizer::minimize() {
775  
776   void OOPSEMinimizer::minimize() {
777      int convgStatus;
840 –
778      int stepStatus;
842 –
779      int maxIter;
844 –
780      int writeFrq;
846 –
781      int nextWriteIter;
782 <
783 <    if (bVerbose)
784 <        printMinimizerInfo();
785 <
786 <    dumpOut = new DumpWriter(info);
787 <
854 <    statOut = new StatWriter(info);
782 >    Snapshot* curSnapshot =info->getSnapshotManager()->getCurrentSnapshot();
783 >    DumpWriter dumpWriter(info, info->getDumpFileName());    
784 >    StatsBitSet mask;
785 >    mask.set(Stats::TIME);
786 >    mask.set(Stats::POTENTIAL_ENERGY);
787 >    StatWriter statWriter(info->getStatFileName(), mask);
788  
789      init();
790  
# Line 864 | Line 797 | void OOPSEMinimizer::minimize() {
797      for(curIter = 1; curIter <= maxIter; curIter++) {
798          stepStatus = step();
799  
800 <        if (nConstrained && bShake)
801 <            preMove();
800 >        //if (usingShake)
801 >        //    preMove();
802  
803          if (stepStatus < 0) {
804              saveResult();
# Line 879 | Line 812 | void OOPSEMinimizer::minimize() {
812              return;
813          }
814  
815 +        //save snapshot
816 +        info->getSnapshotManager()->advance();
817 +        //increase time
818 +        curSnapshot->increaseTime(1);    
819 +        
820          if (curIter == nextWriteIter) {
821              nextWriteIter += writeFrq;
822 <
823 <            writeOut(curX, curIter);
822 >            calcF();
823 >            dumpWriter.writeDump();
824 >            statWriter.writeStat(curSnapshot->statData);
825          }
826  
827          convgStatus = checkConvg();
# Line 899 | Line 838 | void OOPSEMinimizer::minimize() {
838      }
839  
840      if (bVerbose) {
841 <        cout << "OOPSEMinimizer Warning: " << minimizerName
841 >        std::cout << "OOPSEMinimizer Warning: " << minimizerName
842              << " algorithm did not converge within " << maxIter << " iteration"
843              << std::endl;
844      }
# Line 909 | Line 848 | void OOPSEMinimizer::minimize() {
848      saveResult();
849   }
850  
851 +
852 + double OOPSEMinimizer::calcPotential() {
853 +    forceMan->calcForces(true, false);
854 +
855 +    Snapshot* curSnapshot = info->getSnapshotManager()->getCurrentSnapshot();
856 +    double potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] +
857 +                                             curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;    
858 +    double potential;
859 +
860 + #ifdef IS_MPI
861 +    MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM,
862 +                  MPI_COMM_WORLD);
863 + #else
864 +    potential = potential_local;
865 + #endif
866 +
867 +    //save total potential
868 +    curSnapshot->statData[Stats::POTENTIAL_ENERGY] = potential;
869 +    return potential;
870   }
871 +
872 + }

Diff Legend

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