Commit: d5b846f87409da25db9520322b2640389622bb03 - MolDS (git)

Commit MetaInfo

• Revision: d5b846f87409da25db9520322b2640389622bb03 (tree)
• Time: 2014-01-11 03:44:51
• Author: ktns <ktns@1136...>
• Commiter: ktns

Log Message

Refactor BFGS::SearchMinimum using BFGSState. #32881

git-svn-id: https://svn.sourceforge.jp/svnroot/molds/branches/refactor_opt@1638 1136aad2-a195-0410-b898-f5ea1d11b9d8

Change Summary

• modified: src/optimization/BFGS.cpp
• modified: src/optimization/BFGS.h
• modified: src/optimization/Optimizer.h

Incremental Difference

--- a/src/optimization/BFGS.cpp
+++ b/src/optimization/BFGS.cpp
@@ -52,11 +52,35 @@ using namespace MolDS_base;
 using namespace MolDS_base_atoms;
 
 namespace MolDS_optimization{
+
+BFGS::BFGSState::BFGSState(Molecule& molecule):
+      matrixHessian(NULL),
+      matrixOldForce(NULL),
+      matrixStep(NULL),
+      matrixOldCoordinates(NULL),
+      matrixDisplacement(NULL),
+      numAtoms(molecule.GetAtomVect().size())
+{
+   const int dimension = numAtoms * CartesianType_end;
+   MallocerFreer::GetInstance()->Malloc(&this->matrixHessian,        dimension,      dimension);
+   MallocerFreer::GetInstance()->Malloc(&this->matrixOldForce,       this->numAtoms, CartesianType_end);
+   MallocerFreer::GetInstance()->Malloc(&this->matrixStep,           this->numAtoms, CartesianType_end);
+   MallocerFreer::GetInstance()->Malloc(&this->matrixOldCoordinates, this->numAtoms, CartesianType_end);
+   MallocerFreer::GetInstance()->Malloc(&this->matrixDisplacement,   this->numAtoms, CartesianType_end);
+}
+BFGS::BFGSState::~BFGSState(){
+   const int dimension = numAtoms * CartesianType_end;
+   MallocerFreer::GetInstance()->Free(&this->matrixHessian,        dimension,      dimension);
+   MallocerFreer::GetInstance()->Free(&this->matrixOldForce,       this->numAtoms, CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&this->matrixStep,           this->numAtoms, CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&this->matrixOldCoordinates, this->numAtoms, CartesianType_end);
+   MallocerFreer::GetInstance()->Free(&this->matrixDisplacement,   this->numAtoms, CartesianType_end);
+}
+
 BFGS::BFGS(){
    this->SetMessages();
    //this->OutputLog("BFGS created\n");
 }
-
 BFGS::~BFGS(){
    //this->OutputLog("BFGS deleted\n");
 }
@@ -110,129 +134,98 @@ void BFGS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStruct
    int totalSteps = Parameters::GetInstance()->GetTotalStepsOptimization();
    double maxGradientThreshold = Parameters::GetInstance()->GetMaxGradientOptimization();
    double rmsGradientThreshold = Parameters::GetInstance()->GetRmsGradientOptimization();
-   double lineSearchCurrentEnergy   = 0.0;
-   double lineSearchInitialEnergy   = 0.0;
-   double const* const* matrixForce = NULL;
-   double const* vectorForce        = NULL;
    const int dimension = molecule.GetAtomVect().size()*CartesianType_end;
-   double** matrixHessian        = NULL;
-   double*  vectorOldForce       = NULL;
-   double*  vectorStep           = NULL;
-   double** matrixStep           = NULL;
-   double** matrixOldCoordinates = NULL;
-   double*  vectorOldCoordinates = NULL;
-   double** matrixDisplacement   = NULL;
    double       trustRadius      = Parameters::GetInstance()->GetInitialTrustRadiusOptimization();
    const double maxNormStep      = Parameters::GetInstance()->GetMaxNormStepOptimization();
+   BFGSState state(molecule);
 
-   try{
-      // initialize Hessian with unit matrix
-      MallocerFreer::GetInstance()->Malloc(&matrixHessian, dimension, dimension);
-      const double one = 1;
-      MolDS_wrappers::Blas::GetInstance()->Dcopy(dimension, &one, 0, &matrixHessian[0][0], dimension+1);
-
-      // initial calculation
-      bool requireGuess = true;
-      this->UpdateElectronicStructure(electronicStructure, molecule, requireGuess, this->CanOutputLogs());
-      lineSearchCurrentEnergy = electronicStructure->GetElectronicEnergy(elecState);
-
-      requireGuess = false;
-      matrixForce = electronicStructure->GetForce(elecState);
-      vectorForce = &matrixForce[0][0];
-
-      for(int s=0; s<totalSteps; s++){
-         this->OutputLog(boost::format("%s%d\n\n") % this->messageStartBFGSStep % (s+1));
-
-         // Store old Force data
-         MallocerFreer::GetInstance()->Malloc(&vectorOldForce, dimension);
-         MolDS_wrappers::Blas::GetInstance()->Dcopy(dimension,
-                                                    static_cast<double const*>(vectorForce),
-                                                    vectorOldForce);
-
-         this->StoreMolecularGeometry(matrixOldCoordinates, molecule);
-
-         // Level shift Hessian redundant modes
-         this->ShiftHessianRedundantMode(matrixHessian, molecule);
-
-         // Limit the trustRadius to maxNormStep
-         trustRadius=min(trustRadius,maxNormStep);
-
-         //Calculate RFO step
-         MallocerFreer::GetInstance()->Malloc(&matrixStep, molecule.GetAtomVect().size(), CartesianType_end);
-         vectorStep = &matrixStep[0][0];
-         this->CalcRFOStep(vectorStep, matrixHessian, vectorForce, trustRadius, dimension);
-
-         double approximateChange = this->ApproximateEnergyChange(dimension, matrixHessian, vectorForce, vectorStep);
-
-         // Take a RFO step
-         bool doLineSearch = false;
-         bool tempCanOutputLogs = false;
-         lineSearchInitialEnergy = lineSearchCurrentEnergy;
-         if(doLineSearch){
-            this->LineSearch(electronicStructure, molecule, lineSearchCurrentEnergy, matrixStep, elecState, dt);
-         }
-         else{
-            this->UpdateMolecularCoordinates(molecule, matrixStep);
+   // initialize Hessian with unit matrix
+   const double one = 1;
+   MolDS_wrappers::Blas::GetInstance()->Dcopy(dimension, &one, 0, &state.GetMatrixHessian()[0][0], dimension+1);
 
-            // Broadcast to all processes
-            int root = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
-            molecule.BroadcastConfigurationToAllProcesses(root);
+   // initial calculation
+   bool requireGuess = true;
+   this->UpdateElectronicStructure(electronicStructure, molecule, requireGuess, this->CanOutputLogs());
+   state.SetCurrentEnergy(electronicStructure->GetElectronicEnergy(elecState));
 
-            this->UpdateElectronicStructure(electronicStructure, molecule, requireGuess, tempCanOutputLogs);
-            lineSearchCurrentEnergy = electronicStructure->GetElectronicEnergy(elecState);
-         }
-         this->OutputMoleculeElectronicStructure(electronicStructure, molecule, this->CanOutputLogs());
-
-         this->UpdateTrustRadius(trustRadius, approximateChange, lineSearchInitialEnergy, lineSearchCurrentEnergy);
-
-         // check convergence
-         if(this->SatisfiesConvergenceCriterion(matrixForce,
-                  molecule,
-                  lineSearchInitialEnergy,
-                  lineSearchCurrentEnergy,
-                  maxGradientThreshold,
-                  rmsGradientThreshold)){
-            *obtainesOptimizedStructure = true;
-            break;
-         }
+   requireGuess = false;
+   state.SetMatrixForce(electronicStructure->GetForce(elecState));
 
-         //Calculate displacement (K_k at Eq. (15) in [SJTO_1983])
-         this->CalcDisplacement(matrixDisplacement, matrixOldCoordinates, molecule);
+   for(int s=0; s<totalSteps; s++){
+      this->OutputLog(boost::format("%s%d\n\n") % this->messageStartBFGSStep % (s+1));
 
-         //Rollback geometry and energy if energy increases
-         bool isHillClimbing = lineSearchCurrentEnergy > lineSearchInitialEnergy;
-         if(isHillClimbing){
-            this->OutputLog(this->messageHillClimbing);
-            this->RollbackMolecularGeometry(molecule, matrixOldCoordinates);
-            lineSearchCurrentEnergy = lineSearchInitialEnergy;
-         }
+      // Store old Force data
+      MolDS_wrappers::Blas::GetInstance()->Dcopy(dimension,
+                                                 static_cast<double const*>(state.GetVectorForce()),
+                                                 state.GetVectorOldForce());
 
-         matrixForce = electronicStructure->GetForce(elecState);
-         vectorForce = &matrixForce[0][0];
+      this->StoreMolecularGeometry(state.GetMatrixOldCoordinatesRef(), molecule);
 
-         // Update Hessian
-         this->UpdateHessian(matrixHessian, dimension, vectorForce, vectorOldForce, &matrixDisplacement[0][0]);
+      // Level shift Hessian redundant modes
+      this->ShiftHessianRedundantMode(state.GetMatrixHessian(), molecule);
 
-         //Rollback gradient if energy increases
-         if(isHillClimbing){
-            vectorForce = vectorOldForce;
-         }
+      // Limit the trustRadius to maxNormStep
+      trustRadius=min(trustRadius,maxNormStep);
+
+      //Calculate RFO step
+      this->CalcRFOStep(state.GetVectorStep(), state.GetMatrixHessian(), state.GetVectorForce(), trustRadius, dimension);
+
+      double approximateChange = this->ApproximateEnergyChange(dimension, state.GetMatrixHessian(), state.GetVectorForce(), state.GetVectorStep());
+
+      // Take a RFO step
+      bool doLineSearch = false;
+      bool tempCanOutputLogs = false;
+      state.SetInitialEnergy(state.GetCurrentEnergy());
+      if(doLineSearch){
+         this->LineSearch(electronicStructure, molecule, state.GetCurrentEnergyRef(), state.GetMatrixStep(), elecState, dt);
+      }
+      else{
+         this->UpdateMolecularCoordinates(molecule, state.GetMatrixStep());
+
+         // Broadcast to all processes
+         int root = MolDS_mpi::MpiProcess::GetInstance()->GetHeadRank();
+         molecule.BroadcastConfigurationToAllProcesses(root);
+
+         this->UpdateElectronicStructure(electronicStructure, molecule, requireGuess, tempCanOutputLogs);
+         state.SetCurrentEnergy(electronicStructure->GetElectronicEnergy(elecState));
+      }
+      this->OutputMoleculeElectronicStructure(electronicStructure, molecule, this->CanOutputLogs());
+
+      this->UpdateTrustRadius(trustRadius, approximateChange, state.GetInitialEnergy(), state.GetCurrentEnergy());
+
+      // check convergence
+      if(this->SatisfiesConvergenceCriterion(state.GetMatrixForce(),
+               molecule,
+               state.GetInitialEnergy(),
+               state.GetCurrentEnergy(),
+               maxGradientThreshold,
+               rmsGradientThreshold)){
+         *obtainesOptimizedStructure = true;
+         break;
+      }
+
+      //Calculate displacement (K_k at Eq. (15) in [SJTO_1983])
+      this->CalcDisplacement(state.GetMatrixDisplacement(), state.GetMatrixOldCoordinates(), molecule);
+
+      //Rollback geometry and energy if energy increases
+      bool isHillClimbing = state.GetCurrentEnergy() > state.GetInitialEnergy();
+      if(isHillClimbing){
+         this->OutputLog(this->messageHillClimbing);
+         this->RollbackMolecularGeometry(molecule, state.GetMatrixOldCoordinates());
+         state.SetCurrentEnergy(state.GetInitialEnergy());
+      }
+
+      state.SetMatrixForce(electronicStructure->GetForce(elecState));
+
+      // Update Hessian
+      this->UpdateHessian(state.GetMatrixHessian(), dimension, state.GetVectorForce(), state.GetVectorOldForce(), &state.GetMatrixDisplacement()[0][0]);
+
+      //Rollback gradient if energy increases
+      if(isHillClimbing){
+         state.SetMatrixForce(state.GetMatrixOldForce());
       }
-      *lineSearchedEnergy = lineSearchCurrentEnergy;
-   }
-   catch(MolDSException ex){
-      MallocerFreer::GetInstance()->Free(&matrixHessian, dimension, dimension);
-      MallocerFreer::GetInstance()->Free(&vectorOldForce, dimension);
-      MallocerFreer::GetInstance()->Free(&matrixStep, molecule.GetAtomVect().size(), CartesianType_end);
-      MallocerFreer::GetInstance()->Free(&matrixDisplacement, molecule.GetAtomVect().size(), CartesianType_end);
-      MallocerFreer::GetInstance()->Free(&matrixOldCoordinates, molecule.GetAtomVect().size(), CartesianType_end);
-      throw ex;
    }
-   MallocerFreer::GetInstance()->Free(&matrixHessian, dimension, dimension);
-   MallocerFreer::GetInstance()->Free(&vectorOldForce, dimension);
-   MallocerFreer::GetInstance()->Free(&matrixStep, molecule.GetAtomVect().size(), CartesianType_end);
-   MallocerFreer::GetInstance()->Free(&matrixDisplacement, molecule.GetAtomVect().size(), CartesianType_end);
-   MallocerFreer::GetInstance()->Free(&matrixOldCoordinates, molecule.GetAtomVect().size(), CartesianType_end);
+   *lineSearchedEnergy = state.GetCurrentEnergy();
 }
 
 void BFGS::CalcRFOStep(double* vectorStep,
@@ -566,11 +559,10 @@ void BFGS::RollbackMolecularGeometry(MolDS_base::Molecule& molecule,
    molecule.SetCanOutputLogs(tempCanOutputLogs);
 }
 
-void BFGS::CalcDisplacement(double      *      *& matrixDisplacement,
-                            double const* const*  matrixOldCoordinates,
+void BFGS::CalcDisplacement(double      *      * matrixDisplacement,
+                            double const* const* matrixOldCoordinates,
                             const MolDS_base::Molecule& molecule)const{
    //Calculate displacement (K_k at Eq. (15) in [SJTO_1983])
-   MallocerFreer::GetInstance()->Malloc(&matrixDisplacement, molecule.GetAtomVect().size(), CartesianType_end);
    for(int i=0;i<molecule.GetAtomVect().size();i++){
       const Atom*   atom = molecule.GetAtomVect()[i];
       const double* xyz  = atom->GetXyz();

--- a/src/optimization/BFGS.h
+++ b/src/optimization/BFGS.h
@@ -22,6 +22,33 @@
 namespace MolDS_optimization{
 
 class BFGS : public MolDS_optimization::Optimizer{
+private:
+   class BFGSState: public OptimizerState{
+   protected:
+      double** matrixHessian;
+      double** matrixOldForce;
+      double** matrixStep;
+      double** matrixOldCoordinates;
+      double*  vectorOldCoordinates;
+      double** matrixDisplacement;
+      size_t numAtoms;
+   private:
+      template<class vector>
+      vector Matrix2Vector(vector const* matrix){return matrix == NULL ? NULL : &matrix[0][0];}
+   public:
+      BFGSState(MolDS_base::Molecule& molecule);
+      virtual ~BFGSState();
+      double const* GetVectorForce         (){return this->Matrix2Vector(this->matrixForce);}
+      double*       GetVectorOldForce      (){return this->Matrix2Vector(this->matrixOldForce);}
+      double*       GetVectorStep          (){return this->Matrix2Vector(this->matrixStep);}
+      double const* GetVectorOldCoordinates(){return this->Matrix2Vector(this->matrixOldCoordinates);}
+      double**  GetMatrixHessian()          {return this->matrixHessian;}
+      double**  GetMatrixOldForce()         {return this->matrixOldForce;}
+      double**  GetMatrixStep()             {return this->matrixStep;}
+      double**  GetMatrixOldCoordinates()   {return this->matrixOldCoordinates;}
+      double**& GetMatrixOldCoordinatesRef(){return this->matrixOldCoordinates;}
+      double**  GetMatrixDisplacement()     {return this->matrixDisplacement;}
+   };
 public:
    BFGS();
    ~BFGS();
@@ -56,8 +83,8 @@ protected:
                     double const* vectorForce,
                     const double maxNormStep,
                     const int dimension) const;
-   void CalcDisplacement(double      *      *& matrixDisplacement,
-                         double const* const*  matrixOldCoordinates,
+   void CalcDisplacement(double      *      * matrixDisplacement,
+                         double const* const* matrixOldCoordinates,
                          const MolDS_base::Molecule& molecule)const;
    void UpdateHessian(double**      matrixHessian,
                       const int     dimension,

--- a/src/optimization/Optimizer.h
+++ b/src/optimization/Optimizer.h
@@ -24,7 +24,7 @@ namespace MolDS_optimization{
 class Optimizer : public MolDS_base::PrintController{
 protected:
    class OptimizerState{
-   private:
+   protected:
       double currentEnergy;
       double initialEnergy;
       double const* const* matrixForce;