Refactor the relax module

python/pyabacus/src/ModuleDriver/py_driver.cpp

@@ -442,7 +442,7 @@ CalculationResult PyDriver::run(
         if (cal == "scf" || cal == "relax" || cal == "cell-relax" || cal == "nscf")
         {
             Relax_Driver rl_driver;
-            rl_driver.relax_driver(impl_->p_esolver_, *impl_->ucell_, PARAM.inp);
+            rl_driver.relax_driver(impl_->p_esolver_, *impl_->ucell_, PARAM.inp, GlobalV::ofs_running);
         }
         else if (cal == "get_s")
         {

source/Makefile.Objects

@@ -460,10 +460,13 @@ OBJS_PW=fft_bundle.o\
     pw_transform.o\
     pw_transform_k.o\
 
-OBJS_RELAXATION=bfgs_basic.o\
+OBJS_RELAXATION=relax_data.o\
+    cg_base.o\
+    bfgs_basic.o\
     relax_driver.o\
     ions_move_basic.o\
     ions_move_bfgs.o\
+    ions_move_bfgs2.o\
     ions_move_cg.o\
     ions_move_methods.o\
     ions_move_sd.o\
@@ -472,8 +475,7 @@ OBJS_RELAXATION=bfgs_basic.o\
     lattice_change_methods.o\
     relax_nsync.o\
     relax_sync.o\
-    bfgs.o\
-    lbfgs.o\
+    ions_move_lbfgs.o\
     matrix_methods.o\
     line_search.o\
 

source/source_main/driver_run.cpp

@@ -81,7 +81,7 @@ void Driver::driver_run()
     else if (cal == "scf" || cal == "relax" || cal == "cell-relax" || cal == "nscf")
     {
         Relax_Driver rl_driver;
-        rl_driver.relax_driver(p_esolver, ucell, PARAM.inp);
+        rl_driver.relax_driver(p_esolver, ucell, PARAM.inp, GlobalV::ofs_running);
     }
     else if (cal == "get_s")
     {

source/source_relax/CMakeLists.txt

@@ -1,18 +1,17 @@
 add_library(
     relax
     OBJECT
-
+    relax_data.cpp
+    cg_base.cpp
     relax_driver.cpp
-
     relax_sync.cpp
     line_search.cpp
-
-    bfgs.cpp
-    lbfgs.cpp
+    ions_move_lbfgs.cpp
     relax_nsync.cpp
     bfgs_basic.cpp
     ions_move_basic.cpp
     ions_move_bfgs.cpp
+    ions_move_bfgs2.cpp
     ions_move_cg.cpp
     ions_move_sd.cpp
     ions_move_methods.cpp

source/source_relax/bfgs_basic.cpp

@@ -1,82 +1,49 @@
 #include "bfgs_basic.h"
-
+#include <algorithm>
 #include "source_io/module_parameter/parameter.h"
 #include "ions_move_basic.h"
 #include "source_base/global_function.h"
 #include "source_base/global_variable.h"
-#include<vector>
+
 using namespace Ions_Move_Basic;
 
 double BFGS_Basic::relax_bfgs_w1 = -1.0; // default is 0.01
 double BFGS_Basic::relax_bfgs_w2 = -1.0; // defalut is 0.05
 
 BFGS_Basic::BFGS_Basic()
 {
-    pos = nullptr;
-    pos_p = nullptr;
-    grad = nullptr;
-    grad_p = nullptr;
-    move = nullptr;
-    move_p = nullptr;
-
     bfgs_ndim = 1;
 }
 
-BFGS_Basic::~BFGS_Basic()
-{
-    delete[] pos;
-    delete[] pos_p;
-    delete[] grad;
-    delete[] grad_p;
-    delete[] move;
-    delete[] move_p;
-}
-
 void BFGS_Basic::allocate_basic(void)
 {
     assert(dim > 0);
 
-    delete[] pos;
-    delete[] pos_p;
-    delete[] grad;
-    delete[] grad_p;
-    delete[] move;
-    delete[] move_p;
-
-    pos = new double[dim];
-    pos_p = new double[dim];
-    grad = new double[dim];
-    grad_p = new double[dim];
-    move = new double[dim];
-    move_p = new double[dim];
-
-    ModuleBase::GlobalFunc::ZEROS(pos, dim);
-    ModuleBase::GlobalFunc::ZEROS(grad, dim);
-    ModuleBase::GlobalFunc::ZEROS(pos_p, dim);
-    ModuleBase::GlobalFunc::ZEROS(grad_p, dim);
-    ModuleBase::GlobalFunc::ZEROS(move, dim);
-    ModuleBase::GlobalFunc::ZEROS(move_p, dim);
+    pos.resize(dim, 0.0);
+    pos_p.resize(dim, 0.0);
+    grad.resize(dim, 0.0);
+    grad_p.resize(dim, 0.0);
+    move.resize(dim, 0.0);
+    move_p.resize(dim, 0.0);
 
     // init inverse Hessien matrix.
     inv_hess.create(dim, dim);
 
     return;
 }
 
-void BFGS_Basic::update_inverse_hessian(const double &lat0)
+void BFGS_Basic::update_inverse_hessian(const double &lat0, std::ofstream& ofs)
 {
     //  ModuleBase::TITLE("Ions_Move_BFGS","update_inverse_hessian");
     assert(dim > 0);
 
-    std::vector<double> s(dim);
-    std::vector<double> y(dim);
-    ModuleBase::GlobalFunc::ZEROS(s.data(), dim);
-    ModuleBase::GlobalFunc::ZEROS(y.data(), dim);
+    std::vector<double> s(dim, 0.0);
+    std::vector<double> y(dim, 0.0);
 
     for (int i = 0; i < dim; i++)
     {
         //      s[i] = this->pos[i] - this->pos_p[i];
-        //		mohan update 2010-07-27
+        //        mohan update 2010-07-27
         s[i] = this->check_move(lat0, pos[i], pos_p[i]);
         s[i] *= lat0;
 
@@ -90,18 +57,15 @@ void BFGS_Basic::update_inverse_hessian(const double &lat0)
     }
     if (std::abs(sdoty) < 1.0e-16)
     {
-        GlobalV::ofs_running << " WARINIG: unexpected behaviour in update_inverse_hessian" << std::endl;
-        GlobalV::ofs_running << " Resetting bfgs history " << std::endl;
+        ofs << " WARINIG: unexpected behaviour in update_inverse_hessian" << std::endl;
+        ofs << " Resetting bfgs history " << std::endl;
         this->reset_hessian();
         return;
     }
 
-    std::vector<double> Hs(dim);
-    std::vector<double> Hy(dim);
-    std::vector<double> yH(dim);
-    ModuleBase::GlobalFunc::ZEROS(Hs.data(), dim);
-    ModuleBase::GlobalFunc::ZEROS(Hy.data(), dim);
-    ModuleBase::GlobalFunc::ZEROS(yH.data(), dim);
+    std::vector<double> Hs(dim, 0.0);
+    std::vector<double> Hy(dim, 0.0);
+    std::vector<double> yH(dim, 0.0);
 
     for (int i = 0; i < dim; i++)
     {
@@ -132,10 +96,10 @@ void BFGS_Basic::update_inverse_hessian(const double &lat0)
     return;
 }
 
-void BFGS_Basic::check_wolfe_conditions(void)
+void BFGS_Basic::check_wolfe_conditions(std::ofstream& ofs)
 {
-    double dot_p = dot_func(grad_p, move_p, dim);
-    double dot = dot_func(grad, move_p, dim);
+    double dot_p = dot_func(grad_p.data(), move_p.data(), dim);
+    double dot = dot_func(grad.data(), move_p.data(), dim);
 
     // if the total energy falls rapidly, enlarge the trust radius.
     bool wolfe1 = (etot - etot_p) < this->relax_bfgs_w1 * dot_p;
@@ -147,14 +111,14 @@ void BFGS_Basic::check_wolfe_conditions(void)
 
     if (PARAM.inp.test_relax_method)
     {
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "etot - etot_p", etot - etot_p);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "relax_bfgs_w1 * dot_p", relax_bfgs_w1 * dot_p);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "dot", dot);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "relax_bfgs_w2 * dot_p", relax_bfgs_w2 * dot_p);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "relax_bfgs_w1", relax_bfgs_w1);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "relax_bfgs_w2", relax_bfgs_w2);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wolfe1", wolfe1);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wolfe2", wolfe2);
+        ModuleBase::GlobalFunc::OUT(ofs, "etot - etot_p", etot - etot_p);
+        ModuleBase::GlobalFunc::OUT(ofs, "relax_bfgs_w1 * dot_p", relax_bfgs_w1 * dot_p);
+        ModuleBase::GlobalFunc::OUT(ofs, "dot", dot);
+        ModuleBase::GlobalFunc::OUT(ofs, "relax_bfgs_w2 * dot_p", relax_bfgs_w2 * dot_p);
+        ModuleBase::GlobalFunc::OUT(ofs, "relax_bfgs_w1", relax_bfgs_w1);
+        ModuleBase::GlobalFunc::OUT(ofs, "relax_bfgs_w2", relax_bfgs_w2);
+        ModuleBase::GlobalFunc::OUT(ofs, "wolfe1", wolfe1);
+        ModuleBase::GlobalFunc::OUT(ofs, "wolfe2", wolfe2);
     }
 
     this->wolfe_flag = wolfe1 && wolfe2;
@@ -171,13 +135,13 @@ void BFGS_Basic::check_wolfe_conditions(void)
        }
      */
 
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "etot - etot_p", etot - etot_p);
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "relax_bfgs_w1 * dot_p", relax_bfgs_w1 * dot_p);
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wolfe1", wolfe1);
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wolfe2", wolfe2);
+    ModuleBase::GlobalFunc::OUT(ofs, "etot - etot_p", etot - etot_p);
+    ModuleBase::GlobalFunc::OUT(ofs, "relax_bfgs_w1 * dot_p", relax_bfgs_w1 * dot_p);
+    ModuleBase::GlobalFunc::OUT(ofs, "wolfe1", wolfe1);
+    ModuleBase::GlobalFunc::OUT(ofs, "wolfe2", wolfe2);
     //  ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"dot = ",dot);
     //  ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"dot_p = ",dot_p);
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wolfe condition satisfied", wolfe_flag);
+    ModuleBase::GlobalFunc::OUT(ofs, "wolfe condition satisfied", wolfe_flag);
     return;
 }
 
@@ -211,7 +175,7 @@ void BFGS_Basic::save_bfgs(void)
 // a new bfgs step is done
 // we have already done well in the previous direction
 // we should get a new direction in this case
-void BFGS_Basic::new_step(const double &lat0)
+void BFGS_Basic::new_step(const double &lat0, std::ofstream& ofs)
 {
     ModuleBase::TITLE("BFGS_Basic", "new_step");
 
@@ -239,8 +203,8 @@ void BFGS_Basic::new_step(const double &lat0)
     }
     else if (Ions_Move_Basic::update_iter > 1)
     {
-        this->check_wolfe_conditions();
-        this->update_inverse_hessian(lat0);
+        this->check_wolfe_conditions(ofs);
+        this->update_inverse_hessian(lat0, ofs);
     }
 
     //--------------------------------------------------------------------
@@ -263,7 +227,8 @@ void BFGS_Basic::new_step(const double &lat0)
 
             // std::cout << " move after hess " << move[i] << std::endl;
         }
-        GlobalV::ofs_running << " check the norm of new move " << dot_func(move, move, dim) << " (Bohr)" << std::endl;
+
+        ofs << " check the norm of new move " << dot_func(move.data(), move.data(), dim) << " (Bohr)" << std::endl;
     }
     else if (bfgs_ndim > 1)
     {
@@ -280,7 +245,7 @@ void BFGS_Basic::new_step(const double &lat0)
 
     if (dot > 0.0)
     {
-        GlobalV::ofs_running << " Uphill move : resetting bfgs history" << std::endl;
+        ofs << " Uphill move : resetting bfgs history" << std::endl;
         for (int i = 0; i < dim; i++)
         {
             move[i] = -grad[i];
@@ -300,28 +265,28 @@ void BFGS_Basic::new_step(const double &lat0)
     else if (Ions_Move_Basic::update_iter > 1)
     {
         trust_radius = trust_radius_old;
-        this->compute_trust_radius();
+        this->compute_trust_radius(ofs);
     }
     // std::cout<<"trust_radius ="<<" "<<trust_radius;
     return;
 }
 
 // trust radius is computed in this function
 // trust radius determine the step length
-void BFGS_Basic::compute_trust_radius(void)
+void BFGS_Basic::compute_trust_radius(std::ofstream& ofs)
 {
     ModuleBase::TITLE("BFGS_Basic", "compute_trust_radius");
 
     // (1) judge 1
-    double dot = dot_func(grad_p, move_p, dim);
+    double dot = dot_func(grad_p.data(), move_p.data(), dim);
     bool ltest = (etot - etot_p) < this->relax_bfgs_w1 * dot;
 
     // (2) judge 2
     // calculate the norm of move, which
     // is used to compare to trust_radius_old.
-    double norm_move = dot_func(this->move, this->move, dim);
+    double norm_move = dot_func(this->move.data(), this->move.data(), dim);
     norm_move = std::sqrt(norm_move);
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "move(norm)", norm_move);
+    ModuleBase::GlobalFunc::OUT(ofs, "move(norm)", norm_move);
 
     ltest = ltest && (norm_move > trust_radius_old);
 
@@ -356,11 +321,11 @@ void BFGS_Basic::compute_trust_radius(void)
 
     if (PARAM.inp.test_relax_method)
     {
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "wolfe_flag", wolfe_flag);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "trust_radius_old", trust_radius_old);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "2*a*trust_radius_old", 2.0 * a * trust_radius_old);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "norm_move", norm_move);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Trust_radius (Bohr)", trust_radius);
+        ModuleBase::GlobalFunc::OUT(ofs, "wolfe_flag", wolfe_flag);
+        ModuleBase::GlobalFunc::OUT(ofs, "trust_radius_old", trust_radius_old);
+        ModuleBase::GlobalFunc::OUT(ofs, "2*a*trust_radius_old", 2.0 * a * trust_radius_old);
+        ModuleBase::GlobalFunc::OUT(ofs, "norm_move", norm_move);
+        ModuleBase::GlobalFunc::OUT(ofs, "Trust_radius (Bohr)", trust_radius);
     }
 
     if (trust_radius < relax_bfgs_rmin)
@@ -372,7 +337,7 @@ void BFGS_Basic::compute_trust_radius(void)
             // something is going wrongsomething is going wrong
             ModuleBase::WARNING_QUIT("bfgs", "bfgs history already reset at previous step, we got trapped!");
         }
-        GlobalV::ofs_running << " Resetting BFGS history." << std::endl;
+        ofs << " Resetting BFGS history." << std::endl;
         this->reset_hessian();
         for (int i = 0; i < dim; i++)
         {

source/source_relax/bfgs_basic.h

@@ -1,7 +1,10 @@
 #ifndef BFGS_BASIC
 #define BFGS_BASIC
 
+#include <fstream>
+#include <iostream>
 #include "source_base/matrix.h"
+#include <vector>
 
 // references
 // 1) Roger Fletcher, Practical Methods of Optimization, John Wiley and
@@ -18,21 +21,21 @@ class BFGS_Basic
 
   public:
     BFGS_Basic();
-    ~BFGS_Basic();
+    ~BFGS_Basic() = default;
 
   protected:
     void allocate_basic(void);
-    void new_step(const double& lat0);
+    void new_step(const double& lat0, std::ofstream& ofs);
     void reset_hessian(void);
     void save_bfgs(void);
 
-    double* pos = nullptr;  // std::vector containing 3N coordinates of the system ( x )
-    double* grad = nullptr; // std::vector containing 3N components of ( grad( V(x) ) )
-    double* move = nullptr; // pos = pos_p + move.
+    std::vector<double> pos;  // std::vector containing 3N coordinates of the system ( x )
+    std::vector<double> grad; // std::vector containing 3N components of ( grad( V(x) ) )
+    std::vector<double> move; // pos = pos_p + move.
 
-    double* pos_p = nullptr;  // p: previous
-    double* grad_p = nullptr; // p: previous
-    double* move_p = nullptr;
+    std::vector<double> pos_p;  // p: previous
+    std::vector<double> grad_p; // p: previous
+    std::vector<double> move_p;
 
   public:                        // mohan update 2011-06-12
     static double relax_bfgs_w1; // fixed: parameters for Wolfe conditions.
@@ -52,9 +55,9 @@ class BFGS_Basic
 
     int bfgs_ndim;
 
-    void update_inverse_hessian(const double& lat0);
-    void check_wolfe_conditions(void);
-    void compute_trust_radius(void);
+    void update_inverse_hessian(const double& lat0, std::ofstream& ofs);
+    void check_wolfe_conditions(std::ofstream& ofs);
+    void compute_trust_radius(std::ofstream& ofs);
 };
 
 #endif