Commit: 5675d3b6b291c6af7e507734baa7babba1fa8c2a - MolDS (git)

Commit MetaInfo

• Révision: 5675d3b6b291c6af7e507734baa7babba1fa8c2a (tree)
• l'heure: 2014-01-17 11:24:16
• Auteur: Katsuhiko Nishimra <ktns.87@gmai...>
• Commiter: Katsuhiko Nishimra

Message de Log

Factor out Optimizer::UpdateState. #32881

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

Change Summary

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

Modification

--- a/src/optimization/BFGS.cpp
+++ b/src/optimization/BFGS.cpp
@@ -133,14 +133,12 @@ void BFGS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStruct
                          Molecule& molecule,
                          double* lineSearchedEnergy,
                          bool* obtainesOptimizedStructure) const {
-   const int dimension = molecule.GetAtomVect().size()*CartesianType_end;
    BFGSState state(molecule, electronicStructure);
 
    // initial calculation
    bool requireGuess = true;
    this->UpdateElectronicStructure(electronicStructure, molecule, requireGuess, this->CanOutputLogs());
    state.SetCurrentEnergy(electronicStructure->GetElectronicEnergy(state.GetElecState()));
-
    state.SetMatrixForce(electronicStructure->GetForce(state.GetElecState()));
 
    this->InitializeState(state, molecule);
@@ -152,10 +150,11 @@ void BFGS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStruct
 
       this->CalcNextStepGeometry(molecule, state, electronicStructure, state.GetElecState(), state.GetDeltaT());
 
-      this->UpdateTrustRadius(state.GetTrustRadiusRef(), state.GetApproximateChange(), state.GetInitialEnergy(), state.GetCurrentEnergy());
-
+      state.SetCurrentEnergy(electronicStructure->GetElectronicEnergy(state.GetElecState()));
       state.SetMatrixForce(electronicStructure->GetForce(state.GetElecState()));
 
+      this->UpdateState(state);
+
       // check convergence
       if(this->SatisfiesConvergenceCriterion(state.GetMatrixForce(),
                molecule,
@@ -166,25 +165,6 @@ void BFGS::SearchMinimum(boost::shared_ptr<ElectronicStructure> electronicStruct
          *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());
-      }
-
-      // 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 = state.GetCurrentEnergy();
 }
@@ -255,6 +235,33 @@ void BFGS::CalcNextStepGeometry(Molecule &molecule,
    this->OutputMoleculeElectronicStructure(electronicStructure, molecule, this->CanOutputLogs());
 }
 
+void BFGS::UpdateState(OptimizerState& stateOrig) const{
+   BFGSState& state = stateOrig.CastRef<BFGSState>();
+   const MolDS_wrappers::molds_blas_int dimension = state.GetMolecule().GetAtomVect().size() * CartesianType_end;
+
+   this->UpdateTrustRadius(state.GetTrustRadiusRef(), state.GetApproximateChange(), state.GetInitialEnergy(), state.GetCurrentEnergy());
+
+   //Calculate displacement (K_k at Eq. (15) in [SJTO_1983])
+   this->CalcDisplacement(state.GetMatrixDisplacement(), state.GetMatrixOldCoordinates(), state.GetMolecule());
+
+   //Rollback geometry and energy if energy increases
+   bool isHillClimbing = state.GetCurrentEnergy() > state.GetInitialEnergy();
+   if(isHillClimbing){
+      this->OutputLog(this->messageHillClimbing);
+      this->RollbackMolecularGeometry(state.GetMolecule(), state.GetMatrixOldCoordinates());
+      state.SetCurrentEnergy(state.GetInitialEnergy());
+   }
+
+   state.SetMatrixForce(state.GetElectronicStructure()->GetForce(state.GetElecState()));
+
+   // 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());
+   }
+}
 void BFGS::CalcRFOStep(double* vectorStep,
                        double const* const* matrixHessian,
                        double const* vectorForce,

--- a/src/optimization/BFGS.h
+++ b/src/optimization/BFGS.h
@@ -102,6 +102,7 @@ protected:
                              boost::shared_ptr<MolDS_base::ElectronicStructure> electronicStructure,
                              const int elecState,
                              const double dt) const;
+   void UpdateState(OptimizerState& state) const;
    void CalcRFOStep(double* vectorStep,
                     double const* const* matrixHessian,
                     double const* vectorForce,

--- a/src/optimization/ConjugateGradient.cpp
+++ b/src/optimization/ConjugateGradient.cpp
@@ -91,13 +91,10 @@ void ConjugateGradient::SearchMinimum(boost::shared_ptr<ElectronicStructure> ele
    bool requireGuess = true;
    this->UpdateElectronicStructure(electronicStructure, molecule, requireGuess, this->CanOutputLogs());
    state.SetCurrentEnergy(electronicStructure->GetElectronicEnergy(state.GetElecState()));
-
-   requireGuess = false;
    state.SetMatrixForce(electronicStructure->GetForce(state.GetElecState()));
 
    this->InitializeState(state, molecule);
 
-   // conugate gradient loop
    for(int s=0; s<state.GetTotalSteps(); s++){
       this->OutputOptimizationStepMessage(s);
 
@@ -105,7 +102,10 @@ void ConjugateGradient::SearchMinimum(boost::shared_ptr<ElectronicStructure> ele
 
       this->CalcNextStepGeometry(molecule, state, electronicStructure, state.GetElecState(), state.GetDeltaT());
 
-      this->UpdateSearchDirection(state, electronicStructure, molecule, state.GetElecState());
+      state.SetCurrentEnergy(electronicStructure->GetElectronicEnergy(state.GetElecState()));
+      state.SetMatrixForce(electronicStructure->GetForce(state.GetElecState()));
+
+      this->UpdateState(state);
 
       // check convergence
       if(this->SatisfiesConvergenceCriterion(state.GetMatrixForce(),
@@ -118,7 +118,6 @@ void ConjugateGradient::SearchMinimum(boost::shared_ptr<ElectronicStructure> ele
          break;
       }
    }
-
    *lineSearchedEnergy = state.GetCurrentEnergy();
 }
 
@@ -131,6 +130,19 @@ void ConjugateGradient::InitializeState(OptimizerState &stateOrig, const Molecul
    }
 }
 
+void ConjugateGradient::PrepareState(OptimizerState& stateOrig,
+                        const MolDS_base::Molecule& molecule,
+                        const boost::shared_ptr<MolDS_base::ElectronicStructure> electronicStructure,
+                        const int elecState) const{
+   ConjugateGradientState& state = stateOrig.CastRef<ConjugateGradientState>();
+
+   for(int a=0;a<molecule.GetAtomVect().size();a++){
+      for(int i=0; i<CartesianType_end; i++){
+         state.GetOldMatrixForce()[a][i] = state.GetMatrixForce()[a][i];
+      }
+   }
+}
+
 void ConjugateGradient::CalcNextStepGeometry(Molecule &molecule,
                                              OptimizerState& stateOrig,
                                              boost::shared_ptr<ElectronicStructure> electronicStructure,
@@ -143,17 +155,16 @@ void ConjugateGradient::CalcNextStepGeometry(Molecule &molecule,
    this->LineSearch(electronicStructure, molecule, state.GetCurrentEnergyRef(), state.GetMatrixSearchDirection(), elecState, dt);
 }
 
+void ConjugateGradient::UpdateState(OptimizerState& state) const{
+   this->UpdateSearchDirection(state, state.GetElectronicStructure(), state.GetMolecule(), state.GetElecState());
+}
+
+
 void ConjugateGradient::UpdateSearchDirection(OptimizerState& stateOrig,
                                               boost::shared_ptr<ElectronicStructure> electronicStructure, 
                                               const MolDS_base::Molecule& molecule,
                                               int elecState) const{
    ConjugateGradientState& state = stateOrig.CastRef<ConjugateGradientState>();
-   for(int a=0;a<molecule.GetAtomVect().size();a++){
-      for(int i=0; i<CartesianType_end; i++){
-         state.GetOldMatrixForce()[a][i] = state.GetMatrixForce()[a][i];
-      }
-   }
-   state.SetMatrixForce(electronicStructure->GetForce(elecState));
    double beta=0.0;
    double temp=0.0;
    for(int a=0;a<molecule.GetAtomVect().size();a++){

--- a/src/optimization/ConjugateGradient.h
+++ b/src/optimization/ConjugateGradient.h
@@ -53,12 +53,13 @@ private:
    virtual void PrepareState(OptimizerState& state,
                              const MolDS_base::Molecule& molecule,
                              const boost::shared_ptr<MolDS_base::ElectronicStructure> electronicStructure,
-                             const int elecState) const{};
+                             const int elecState) const;
    void CalcNextStepGeometry(MolDS_base::Molecule &molecule,
                              OptimizerState& state,
                              boost::shared_ptr<MolDS_base::ElectronicStructure> electronicStructure,
                              const int elecState,
                              const double dt) const;
+   void UpdateState(OptimizerState& state) const;
    void UpdateSearchDirection(OptimizerState& state,
                               boost::shared_ptr<MolDS_base::ElectronicStructure> electronicStructure,
                               const MolDS_base::Molecule& molecule,

--- a/src/optimization/Optimizer.h
+++ b/src/optimization/Optimizer.h
@@ -129,6 +129,7 @@ private:
                                      boost::shared_ptr<MolDS_base::ElectronicStructure> electronicStructure,
                                      const int elecState,
                                      const double dt) const= 0;
+   virtual void UpdateState(OptimizerState& state) const = 0;
    virtual const std::string& OptimizationStepMessage() const = 0;
 };
 

--- a/src/optimization/SteepestDescent.cpp
+++ b/src/optimization/SteepestDescent.cpp
@@ -74,17 +74,21 @@ void SteepestDescent::SearchMinimum(boost::shared_ptr<ElectronicStructure> elect
    // initial calculation
    bool requireGuess = true;
    this->UpdateElectronicStructure(electronicStructure, molecule, requireGuess, this->CanOutputLogs());
-
-   requireGuess = false;
-   state.SetMatrixForce(electronicStructure->GetForce(state.GetElecState()));
    state.SetCurrentEnergy(electronicStructure->GetElectronicEnergy(state.GetElecState()));
+   state.SetMatrixForce(electronicStructure->GetForce(state.GetElecState()));
+
+   this->InitializeState(state, molecule);
+
    for(int s=0; s<state.GetTotalSteps(); s++){
       this->OutputLog(boost::format("%s%d\n\n") % this->messageStartSteepestDescentStep.c_str() % (s+1));
       state.SetInitialEnergy(state.GetCurrentEnergy());
 
       this->CalcNextStepGeometry(molecule, state, electronicStructure, state.GetElecState(), state.GetDeltaT());
 
-      this->UpdateSearchDirection(state, electronicStructure, molecule, state.GetElecState());
+      state.SetCurrentEnergy(electronicStructure->GetElectronicEnergy(state.GetElecState()));
+      state.SetMatrixForce(electronicStructure->GetForce(state.GetElecState()));
+
+      this->UpdateState(state);
 
       // check convergence
       if(this->SatisfiesConvergenceCriterion(state.GetMatrixForce(),
@@ -97,7 +101,6 @@ void SteepestDescent::SearchMinimum(boost::shared_ptr<ElectronicStructure> elect
          break;
       }
    }
-
    *lineSearchedEnergy = state.GetCurrentEnergy();
 }
 
@@ -111,6 +114,10 @@ void SteepestDescent::CalcNextStepGeometry(Molecule &molecule,
    this->LineSearch(electronicStructure, molecule, state.GetCurrentEnergyRef(), state.GetMatrixForce(), elecState, dt);
 }
 
+void SteepestDescent::UpdateState(OptimizerState& state) const{
+   this->UpdateSearchDirection(state, state.GetElectronicStructure(), state.GetMolecule(), state.GetElecState());
+}
+
 void SteepestDescent::UpdateSearchDirection(OptimizerState& state,
                                             boost::shared_ptr<ElectronicStructure> electronicStructure,
                                             const MolDS_base::Molecule& molecule,

--- a/src/optimization/SteepestDescent.h
+++ b/src/optimization/SteepestDescent.h
@@ -46,6 +46,7 @@ private:
                              boost::shared_ptr<MolDS_base::ElectronicStructure> electronicStructure,
                              const int elecState,
                              const double dt) const;
+   void UpdateState(OptimizerState& state) const;
    void UpdateSearchDirection(OptimizerState& state,
                               boost::shared_ptr<MolDS_base::ElectronicStructure> electronicStructure,
                               const MolDS_base::Molecule& molecule,