From ef9dd5ae1d04419d452b6358c958881c9f124d6e Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 8 Jun 2026 21:09:12 +0800 Subject: [PATCH 01/16] format update --- source/source_relax/bfgs.cpp | 2 +- source/source_relax/bfgs_basic.cpp | 2 +- source/source_relax/ions_move_bfgs.cpp | 16 ++-- source/source_relax/lattice_change_basic.cpp | 36 ++++---- source/source_relax/line_search.cpp | 14 ++-- source/source_relax/relax_driver.cpp | 84 +++++++++---------- source/source_relax/relax_driver.h | 6 +- source/source_relax/relax_nsync.cpp | 2 +- .../source_relax/test/ions_move_bfgs_test.cpp | 14 ++-- 9 files changed, 88 insertions(+), 88 deletions(-) diff --git a/source/source_relax/bfgs.cpp b/source/source_relax/bfgs.cpp index 13b87f55097..b7fbbe821a1 100644 --- a/source/source_relax/bfgs.cpp +++ b/source/source_relax/bfgs.cpp @@ -320,7 +320,7 @@ void BFGS::UpdatePos(UnitCell& ucell) move_ion[iat * 3 + 2] = move_ion_dr.z ; } } - ucell.update_pos_taud(move_ion); + ucell.update_pos_taud(move_ion); pos = this->MAddM(pos, dpos);*/ } diff --git a/source/source_relax/bfgs_basic.cpp b/source/source_relax/bfgs_basic.cpp index d598519fc3e..bdd7ae67626 100644 --- a/source/source_relax/bfgs_basic.cpp +++ b/source/source_relax/bfgs_basic.cpp @@ -76,7 +76,7 @@ void BFGS_Basic::update_inverse_hessian(const double &lat0) 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; diff --git a/source/source_relax/ions_move_bfgs.cpp b/source/source_relax/ions_move_bfgs.cpp index ee8a9d97e98..d3ff58732ca 100644 --- a/source/source_relax/ions_move_bfgs.cpp +++ b/source/source_relax/ions_move_bfgs.cpp @@ -23,10 +23,10 @@ Ions_Move_BFGS::~Ions_Move_BFGS(){}; void Ions_Move_BFGS::allocate() { ModuleBase::TITLE("Ions_Move_BFGS", "init"); - if (init_done) - { - return; - } + if (init_done) + { + return; + } this->allocate_basic(); // initialize data members @@ -50,9 +50,9 @@ void Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, con Ions_Move_Basic::setup_gradient(ucell, force, this->pos, this->grad); first_step = false; } - else - { - std::vector pos_tmp(3 * ucell.nat); + else + { + std::vector pos_tmp(3 * ucell.nat); Ions_Move_Basic::setup_gradient(ucell, force, pos_tmp.data(), this->grad); } // use energy_in and istep to setup etot and etot_old. @@ -129,7 +129,7 @@ void Ions_Move_BFGS::restart_bfgs(const double& lat0) } else { - // bfgs initialization + // bfgs initialization ModuleBase::GlobalFunc::ZEROS(pos_p, dim); ModuleBase::GlobalFunc::ZEROS(grad_p, dim); ModuleBase::GlobalFunc::ZEROS(move_p, dim); diff --git a/source/source_relax/lattice_change_basic.cpp b/source/source_relax/lattice_change_basic.cpp index 47b0c283e19..1bc41ce6238 100644 --- a/source/source_relax/lattice_change_basic.cpp +++ b/source/source_relax/lattice_change_basic.cpp @@ -117,24 +117,24 @@ void Lattice_Change_Basic::change_lattice(UnitCell &ucell, double *move, double if (ModuleSymmetry::Symmetry::symm_flag && ucell.symm.nrotk > 0) { ModuleBase::matrix move_mat_t(3, 3); - for (int i = 0;i < 3;++i) - { - for (int j = 0;j < 3;++j) - { - move_mat_t(j, i) = move[i * 3 + j] / ucell.lat0; //transpose - } - } - ModuleBase::matrix symm_move_mat_t = (move_mat_t * ucell.G.to_matrix());//symmetrize (latvec^{-1} * move_mat)^T + for (int i = 0;i < 3;++i) + { + for (int j = 0;j < 3;++j) + { + move_mat_t(j, i) = move[i * 3 + j] / ucell.lat0; //transpose + } + } + ModuleBase::matrix symm_move_mat_t = (move_mat_t * ucell.G.to_matrix());//symmetrize (latvec^{-1} * move_mat)^T ucell.symm.symmetrize_mat3(symm_move_mat_t, ucell.lat); move_mat_t = symm_move_mat_t * ucell.latvec.Transpose().to_matrix();//G^{-1}=latvec^T - for (int i = 0;i < 3;++i) - { - for (int j = 0;j < 3;++j) - { - move[i * 3 + j] = move_mat_t(j, i) * ucell.lat0;//transpose back - } - } + for (int i = 0;i < 3;++i) + { + for (int j = 0;j < 3;++j) + { + move[i * 3 + j] = move_mat_t(j, i) * ucell.lat0;//transpose back + } + } } if (ucell.lc[0] != 0) @@ -252,9 +252,9 @@ void Lattice_Change_Basic::check_converged(const UnitCell &ucell, ModuleBase::ma for (int i = 0; i < 3; i++) { if (stress_ii_max < std::abs(stress(i, i))) - { - stress_ii_max = std::abs(stress(i, i)); - } + { + stress_ii_max = std::abs(stress(i, i)); + } for (int j = 0; j < 3; j++) { if (Lattice_Change_Basic::largest_grad < std::abs(stress(i, j))) diff --git a/source/source_relax/line_search.cpp b/source/source_relax/line_search.cpp index 1c07a59c3d4..813fa002cb9 100644 --- a/source/source_relax/line_search.cpp +++ b/source/source_relax/line_search.cpp @@ -15,7 +15,7 @@ bool Line_Search::line_search(const bool restart, { if (restart) { - ls_step = 0; + ls_step = 0; } if (ls_step == 0) // first point: make a trial step into trial direction @@ -103,7 +103,7 @@ bool Line_Search::third_order(const double x, const double y, const double f, do dmoveh = -fa / (fab - fa) / 2.0; if (dmoveh < 0) - { + { dmoveh = 4.0; } @@ -116,7 +116,7 @@ bool Line_Search::third_order(const double x, const double y, const double f, do if (dmove > 4.0) { - dmove = 4.0; + dmove = 4.0; } xnew = dmove + xa; @@ -143,7 +143,7 @@ void Line_Search::init_brent(const double x, const double y, const double f) yb = y; fb = f; if (fa * fb > 0) - { + { bracked = false; } fstart = fa; @@ -240,7 +240,7 @@ bool Line_Search::brent(const double x, const double y, const double f, double& qq = (qq - 1.0) * (r - 1.0) * (s - 1.0); } if (p > 0.0) - { + { qq = -qq; } p = std::abs(p); @@ -281,7 +281,7 @@ bool Line_Search::brent(const double x, const double y, const double f, double& xnew = xb; if (std::abs(dy) < conv_thr) - { + { return true; } if (ls_step == 4) // I'm not sure if this is a good choice, but the idea is there should not be so many line @@ -344,7 +344,7 @@ bool Line_Search::brent(const double x, const double y, const double f, double& qq = (qq - 1.0) * (r - 1.0) * (s - 1.0); } if (p > 0.0) - { + { qq = -qq; } p = std::abs(p); diff --git a/source/source_relax/relax_driver.cpp b/source/source_relax/relax_driver.cpp index ef97baf4dd7..82d96b7b9dd 100644 --- a/source/source_relax/relax_driver.cpp +++ b/source/source_relax/relax_driver.cpp @@ -10,9 +10,9 @@ #include "source_cell/print_cell.h" void Relax_Driver::relax_driver( - ModuleESolver::ESolver* p_esolver, - UnitCell& ucell, - const Input_para& inp) + ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp) { ModuleBase::TITLE("Relax_Driver", "relax_driver"); ModuleBase::timer::start("Relax_Driver", "relax_driver"); @@ -38,13 +38,13 @@ void Relax_Driver::relax_driver( { time_t estart = time(nullptr); - if (inp.out_level == "ie" - && (inp.calculation == "relax" - || inp.calculation == "cell-relax" - || inp.calculation == "scf" - || inp.calculation == "nscf") - && (inp.esolver_type != "lr")) - { + if (inp.out_level == "ie" + && (inp.calculation == "relax" + || inp.calculation == "cell-relax" + || inp.calculation == "scf" + || inp.calculation == "nscf") + && (inp.esolver_type != "lr")) + { ModuleIO::print_screen(stress_step, force_step, istep); } @@ -73,10 +73,10 @@ void Relax_Driver::relax_driver( { p_esolver->cal_force(ucell, force); } - else - { + else + { // do nothing - } + } // calculate and gather all parts of stress @@ -84,10 +84,10 @@ void Relax_Driver::relax_driver( { p_esolver->cal_stress(ucell, stress); } - else - { + else + { // do nothing - } + } if (inp.calculation == "relax" || inp.calculation == "cell-relax") { @@ -181,37 +181,37 @@ void Relax_Driver::relax_driver( "data_?"); } - if (inp.calculation == "relax" || inp.calculation == "cell-relax") - { - if (istep-1 == inp.relax_nmax) - { - std::cout << "\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; - std::cout << " Geometry relaxation stops here due to reaching the maximum " << std::endl; - std::cout << " relaxation steps. More steps are needed to converge the results " << std::endl; - std::cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; - } - else - { - std::cout << "\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; - std::cout << " Geometry relaxation thresholds are reached within " << istep-1 << " steps." << std::endl; - std::cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; - } - } - else - { - // do nothing - } - - if (inp.relax_nmax == 0) + if (inp.calculation == "relax" || inp.calculation == "cell-relax") + { + if (istep-1 == inp.relax_nmax) + { + std::cout << "\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; + std::cout << " Geometry relaxation stops here due to reaching the maximum " << std::endl; + std::cout << " relaxation steps. More steps are needed to converge the results " << std::endl; + std::cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; + } + else + { + std::cout << "\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; + std::cout << " Geometry relaxation thresholds are reached within " << istep-1 << " steps." << std::endl; + std::cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; + } + } + else + { + // do nothing + } + + if (inp.relax_nmax == 0) { std::cout << "-----------------------------------------------" << std::endl; std::cout << " relax_nmax = 0, DRY RUN TEST SUCCEEDS :)" << std::endl; std::cout << "-----------------------------------------------" << std::endl; } - else - { - // do nothing - } + else + { + // do nothing + } ModuleBase::timer::end("Relax_Driver", "relax_driver"); return; diff --git a/source/source_relax/relax_driver.h b/source/source_relax/relax_driver.h index 8063205caaf..78b343efd81 100644 --- a/source/source_relax/relax_driver.h +++ b/source/source_relax/relax_driver.h @@ -16,9 +16,9 @@ class Relax_Driver Relax_Driver(){}; ~Relax_Driver(){}; - void relax_driver(ModuleESolver::ESolver* p_esolver, - UnitCell& ucell, - const Input_para& inp); + void relax_driver(ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp); private: // mohan add 2021-01-28 diff --git a/source/source_relax/relax_nsync.cpp b/source/source_relax/relax_nsync.cpp index d50a07ded03..9541c9dc177 100644 --- a/source/source_relax/relax_nsync.cpp +++ b/source/source_relax/relax_nsync.cpp @@ -89,7 +89,7 @@ bool Relax_old::if_do_relax(const UnitCell& ucell) ModuleBase::WARNING("Ions", "No atom is allowed to move!"); return false; } - // if(!IMM.get_converged()) return 1; + // if(!IMM.get_converged()) return 1; else { assert(PARAM.inp.cal_force == 1); diff --git a/source/source_relax/test/ions_move_bfgs_test.cpp b/source/source_relax/test/ions_move_bfgs_test.cpp index a9094e20212..bbf9bfadf93 100644 --- a/source/source_relax/test/ions_move_bfgs_test.cpp +++ b/source/source_relax/test/ions_move_bfgs_test.cpp @@ -187,13 +187,13 @@ TEST_F(IonsMoveBFGSTest, RestartBfgsCase2) for (int j = 0; j < Ions_Move_Basic::dim; ++j) { if (i == j) - { - EXPECT_DOUBLE_EQ(bfgs.inv_hess(i, j), 1.0); - } - else - { - EXPECT_DOUBLE_EQ(bfgs.inv_hess(i, j), 0.0); - } + { + EXPECT_DOUBLE_EQ(bfgs.inv_hess(i, j), 1.0); + } + else + { + EXPECT_DOUBLE_EQ(bfgs.inv_hess(i, j), 0.0); + } } } } From cfdf83f5b50722e0b2329fa60f6dcb629d808461 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 8 Jun 2026 21:17:47 +0800 Subject: [PATCH 02/16] replace new with vector --- source/source_relax/ions_move_sd.cpp | 28 ++++++++-------------------- source/source_relax/ions_move_sd.h | 8 +++++--- 2 files changed, 13 insertions(+), 23 deletions(-) diff --git a/source/source_relax/ions_move_sd.cpp b/source/source_relax/ions_move_sd.cpp index e3831623050..133f809367c 100644 --- a/source/source_relax/ions_move_sd.cpp +++ b/source/source_relax/ions_move_sd.cpp @@ -4,32 +4,19 @@ #include "ions_move_basic.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" -#include using namespace Ions_Move_Basic; -Ions_Move_SD::Ions_Move_SD() +Ions_Move_SD::Ions_Move_SD() : energy_saved(1.0e10) { - this->energy_saved = 1.0e10; - this->grad_saved = nullptr; - this->pos_saved = nullptr; -} -Ions_Move_SD::~Ions_Move_SD() -{ - delete[] grad_saved; - delete[] pos_saved; } void Ions_Move_SD::allocate() { ModuleBase::TITLE("Ions_Move_SD", "allocate"); assert(dim > 0); - delete[] grad_saved; - delete[] pos_saved; - this->grad_saved = new double[dim]; - this->pos_saved = new double[dim]; - ModuleBase::GlobalFunc::ZEROS(grad_saved, dim); - ModuleBase::GlobalFunc::ZEROS(pos_saved, dim); + grad_saved.resize(dim, 0.0); + pos_saved.resize(dim, 0.0); } void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const double& etot_in) @@ -37,8 +24,8 @@ void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const ModuleBase::TITLE("Ions_Move_SD", "start"); assert(dim > 0); - assert(grad_saved != nullptr); - assert(pos_saved != nullptr); + assert(grad_saved.size() == static_cast(dim)); + assert(pos_saved.size() == static_cast(dim)); std::vector pos(dim); std::vector grad(dim); @@ -58,9 +45,10 @@ void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const printf("in cheak_converged"); printf("pos[0]: %f\n", pos[0]); energy_saved = etot_in; - for (int i = 0; i < dim; i++) { + for (int i = 0; i < dim; i++) + { pos_saved[i] = pos[i]; -} + } for (int i = 0; i < dim; i++) { grad_saved[i] = grad[i]; diff --git a/source/source_relax/ions_move_sd.h b/source/source_relax/ions_move_sd.h index 0fa56f8696a..e2eb53d587c 100644 --- a/source/source_relax/ions_move_sd.h +++ b/source/source_relax/ions_move_sd.h @@ -3,19 +3,21 @@ #include "source_base/matrix.h" #include "source_cell/unitcell.h" +#include + class Ions_Move_SD { public: Ions_Move_SD(); - ~Ions_Move_SD(); + ~Ions_Move_SD() = default; void allocate(void); void start(UnitCell& ucell, const ModuleBase::matrix& force, const double& etot); private: double energy_saved; - double* pos_saved = nullptr; - double* grad_saved = nullptr; + std::vector pos_saved; + std::vector grad_saved; void cal_tradius_sd(void) const; }; From 94eb11bb4a66cd2c85ff9cb009984f18ab2327f5 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 8 Jun 2026 21:24:16 +0800 Subject: [PATCH 03/16] format update --- source/source_relax/line_search.cpp | 32 ++++++++++++++--------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/source/source_relax/line_search.cpp b/source/source_relax/line_search.cpp index 813fa002cb9..db5a049382b 100644 --- a/source/source_relax/line_search.cpp +++ b/source/source_relax/line_search.cpp @@ -13,9 +13,9 @@ bool Line_Search::line_search(const bool restart, double& xnew, // the next point that we want to try const double conv_thr) { - if (restart) + if (restart) { - ls_step = 0; + ls_step = 0; } if (ls_step == 0) // first point: make a trial step into trial direction @@ -102,8 +102,8 @@ bool Line_Search::third_order(const double x, const double y, const double f, do } dmoveh = -fa / (fab - fa) / 2.0; - if (dmoveh < 0) - { + if (dmoveh < 0) + { dmoveh = 4.0; } @@ -114,14 +114,14 @@ bool Line_Search::third_order(const double x, const double y, const double f, do } } // end anharmonic case - if (dmove > 4.0) + if (dmove > 4.0) { - dmove = 4.0; + dmove = 4.0; } xnew = dmove + xa; double dy = (fb + (fab - fb) / (xa - xb) * (dmove - xb)) * (dmove - xb); - if (std::abs(dy) < conv_thr) + if (std::abs(dy) < conv_thr) { return true; } @@ -142,8 +142,8 @@ void Line_Search::init_brent(const double x, const double y, const double f) xb = x; yb = y; fb = f; - if (fa * fb > 0) - { + if (fa * fb > 0) + { bracked = false; } fstart = fa; @@ -239,8 +239,8 @@ bool Line_Search::brent(const double x, const double y, const double f, double& p = s * (2.0 * xm * qq * (qq - r) - (xb - xa) * (r - 1.0)); qq = (qq - 1.0) * (r - 1.0) * (s - 1.0); } - if (p > 0.0) - { + if (p > 0.0) + { qq = -qq; } p = std::abs(p); @@ -280,8 +280,8 @@ bool Line_Search::brent(const double x, const double y, const double f, double& } xnew = xb; - if (std::abs(dy) < conv_thr) - { + if (std::abs(dy) < conv_thr) + { return true; } if (ls_step == 4) // I'm not sure if this is a good choice, but the idea is there should not be so many line @@ -343,8 +343,8 @@ bool Line_Search::brent(const double x, const double y, const double f, double& p = s * (2.0 * xm * qq * (qq - r) - (xb - xa) * (r - 1.0)); qq = (qq - 1.0) * (r - 1.0) * (s - 1.0); } - if (p > 0.0) - { + if (p > 0.0) + { qq = -qq; } p = std::abs(p); @@ -390,7 +390,7 @@ bool Line_Search::brent(const double x, const double y, const double f, double& } xnew = xb; - if (std::abs(dy) < conv_thr) + if (std::abs(dy) < conv_thr) { return true; } From 6fcf5cf5900b6aeb5d9dbf057f8cdbe14f1c4849 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 8 Jun 2026 21:27:46 +0800 Subject: [PATCH 04/16] update --- source/source_relax/bfgs_basic.cpp | 56 +++++++----------------------- source/source_relax/bfgs_basic.h | 15 ++++---- 2 files changed, 20 insertions(+), 51 deletions(-) diff --git a/source/source_relax/bfgs_basic.cpp b/source/source_relax/bfgs_basic.cpp index bdd7ae67626..f542d209a30 100644 --- a/source/source_relax/bfgs_basic.cpp +++ b/source/source_relax/bfgs_basic.cpp @@ -1,10 +1,9 @@ #include "bfgs_basic.h" - #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 + using namespace Ions_Move_Basic; double BFGS_Basic::relax_bfgs_w1 = -1.0; // default is 0.01 @@ -12,50 +11,19 @@ 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); @@ -134,8 +102,8 @@ void BFGS_Basic::update_inverse_hessian(const double &lat0) void BFGS_Basic::check_wolfe_conditions(void) { - 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; @@ -263,7 +231,7 @@ 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; + GlobalV::ofs_running << " check the norm of new move " << dot_func(move.data(), move.data(), dim) << " (Bohr)" << std::endl; } else if (bfgs_ndim > 1) { @@ -313,13 +281,13 @@ void BFGS_Basic::compute_trust_radius(void) 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); diff --git a/source/source_relax/bfgs_basic.h b/source/source_relax/bfgs_basic.h index 85de10474ac..6e74cac9cdd 100644 --- a/source/source_relax/bfgs_basic.h +++ b/source/source_relax/bfgs_basic.h @@ -2,6 +2,7 @@ #define BFGS_BASIC #include "source_base/matrix.h" +#include // references // 1) Roger Fletcher, Practical Methods of Optimization, John Wiley and @@ -18,7 +19,7 @@ class BFGS_Basic public: BFGS_Basic(); - ~BFGS_Basic(); + ~BFGS_Basic() = default; protected: void allocate_basic(void); @@ -26,13 +27,13 @@ class BFGS_Basic 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 pos; // std::vector containing 3N coordinates of the system ( x ) + std::vector grad; // std::vector containing 3N components of ( grad( V(x) ) ) + std::vector move; // pos = pos_p + move. - double* pos_p = nullptr; // p: previous - double* grad_p = nullptr; // p: previous - double* move_p = nullptr; + std::vector pos_p; // p: previous + std::vector grad_p; // p: previous + std::vector move_p; public: // mohan update 2011-06-12 static double relax_bfgs_w1; // fixed: parameters for Wolfe conditions. From 01621e99051e4cf7bfc786e1b97c08f0ba095535 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 8 Jun 2026 21:29:41 +0800 Subject: [PATCH 05/16] update --- source/source_relax/ions_move_bfgs.cpp | 16 +++--- source/source_relax/lattice_change_cg.cpp | 63 ++++++++--------------- source/source_relax/lattice_change_cg.h | 11 ++-- 3 files changed, 35 insertions(+), 55 deletions(-) diff --git a/source/source_relax/ions_move_bfgs.cpp b/source/source_relax/ions_move_bfgs.cpp index d3ff58732ca..31c9f33f86d 100644 --- a/source/source_relax/ions_move_bfgs.cpp +++ b/source/source_relax/ions_move_bfgs.cpp @@ -47,19 +47,19 @@ void Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, con // In the following steps, the pos is updated by BFGS methods. if (first_step) { - Ions_Move_Basic::setup_gradient(ucell, force, this->pos, this->grad); + Ions_Move_Basic::setup_gradient(ucell, force, this->pos.data(), this->grad.data()); first_step = false; } else { std::vector pos_tmp(3 * ucell.nat); - Ions_Move_Basic::setup_gradient(ucell, force, pos_tmp.data(), this->grad); + Ions_Move_Basic::setup_gradient(ucell, force, pos_tmp.data(), this->grad.data()); } // use energy_in and istep to setup etot and etot_old. Ions_Move_Basic::setup_etot(energy_in, false); // use gradient and etot and etot_old to check // if the result is converged. - Ions_Move_Basic::check_converged(ucell, this->grad); + Ions_Move_Basic::check_converged(ucell, this->grad.data()); if (Ions_Move_Basic::converged) { @@ -86,7 +86,7 @@ void Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, con // even if the energy is higher, we save the information. this->save_bfgs(); - Ions_Move_Basic::move_atoms(ucell, move, pos); + Ions_Move_Basic::move_atoms(ucell, move.data(), pos.data()); } return; } @@ -130,9 +130,9 @@ void Ions_Move_BFGS::restart_bfgs(const double& lat0) else { // bfgs initialization - ModuleBase::GlobalFunc::ZEROS(pos_p, dim); - ModuleBase::GlobalFunc::ZEROS(grad_p, dim); - ModuleBase::GlobalFunc::ZEROS(move_p, dim); + ModuleBase::GlobalFunc::ZEROS(pos_p.data(), dim); + ModuleBase::GlobalFunc::ZEROS(grad_p.data(), dim); + ModuleBase::GlobalFunc::ZEROS(move_p.data(), dim); Ions_Move_Basic::update_iter = 0; @@ -304,7 +304,7 @@ void Ions_Move_BFGS::bfgs_routine(const double& lat0) ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "update iteration", Ions_Move_Basic::update_iter); // combine the direction and move length now - double norm = dot_func(this->move, this->move, dim); + double norm = dot_func(this->move.data(), this->move.data(), dim); norm = sqrt(norm); if (norm < 1.0e-16) diff --git a/source/source_relax/lattice_change_cg.cpp b/source/source_relax/lattice_change_cg.cpp index 1ed2449c26a..d0172ba7845 100644 --- a/source/source_relax/lattice_change_cg.cpp +++ b/source/source_relax/lattice_change_cg.cpp @@ -24,18 +24,7 @@ using namespace Lattice_Change_Basic; Lattice_Change_CG::Lattice_Change_CG() { - this->lat0 = nullptr; - this->grad0 = nullptr; - this->cg_grad0 = nullptr; - this->move0 = nullptr; -} - -Lattice_Change_CG::~Lattice_Change_CG() -{ - delete[] lat0; - delete[] grad0; - delete[] cg_grad0; - delete[] move0; + this->e0 = 0.0; } void Lattice_Change_CG::allocate(void) @@ -44,20 +33,10 @@ void Lattice_Change_CG::allocate(void) // mohan add 2021-02-07 assert(dim > 0); - delete[] lat0; - delete[] grad0; - delete[] cg_grad0; - delete[] move0; - - this->lat0 = new double[dim]; - this->grad0 = new double[dim]; - this->cg_grad0 = new double[dim]; - this->move0 = new double[dim]; - - ModuleBase::GlobalFunc::ZEROS(lat0, dim); - ModuleBase::GlobalFunc::ZEROS(grad0, dim); - ModuleBase::GlobalFunc::ZEROS(cg_grad0, dim); - ModuleBase::GlobalFunc::ZEROS(move0, dim); + this->lat0.resize(dim, 0.0); + this->grad0.resize(dim, 0.0); + this->cg_grad0.resize(dim, 0.0); + this->move0.resize(dim, 0.0); this->e0 = 0.0; } @@ -65,10 +44,10 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_ { ModuleBase::TITLE("Lattice_Change_CG", "start"); - assert(lat0 != 0); - assert(grad0 != 0); - assert(cg_grad0 != 0); - assert(move0 != 0); + assert(lat0.size() == static_cast(dim)); + assert(grad0.size() == static_cast(dim)); + assert(cg_grad0.size() == static_cast(dim)); + assert(move0.size() == static_cast(dim)); // sd , trial are two parameters, when sd=trial=true, @@ -158,16 +137,16 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_ { e0 = etot_in; setup_cg_grad(grad, - grad0, + grad0.data(), cg_grad, - cg_grad0, + cg_grad0.data(), ncggrad, flag); // we use the last direction ,the last grad and the grad now to get the direction now ncggrad++; normalize(cg_gradn, cg_grad, dim); - setup_move(move0, cg_gradn, steplength); // move the atom position - Lattice_Change_Basic::change_lattice(ucell, move0, lat); + setup_move(move0.data(), cg_gradn, steplength); // move the atom position + Lattice_Change_Basic::change_lattice(ucell, move0.data(), lat); for (int i = 0; i < dim; i++) // grad0 ,cg_grad0 are used to store the grad and cg_grad for the future using { @@ -175,8 +154,8 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_ cg_grad0[i] = cg_grad[i]; } - f_cal(move0, move0, dim, xb); // xb = trial steplength - f_cal(move0, grad, dim, fa); // fa is the projection force in this direction + f_cal(move0.data(), move0.data(), dim, xb); // xb = trial steplength + f_cal(move0.data(), grad, dim, fa); // fa is the projection force in this direction fmax = fa; sd = false; @@ -188,8 +167,8 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_ if (trial) { double e1 = etot_in; - f_cal(move0, grad, dim, fb); - f_cal(move0, move0, dim, xb); + f_cal(move0.data(), grad, dim, fb); + f_cal(move0.data(), move0.data(), dim, xb); if ((std::abs(fb) < std::abs((fa) / 10.0))) { @@ -200,7 +179,7 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_ goto CG_begin; } - normalize(cg_gradn, cg_grad0, dim); + normalize(cg_gradn, cg_grad0.data(), dim); third_order(e0, e1, fa, fb, xb, best_x); // cubic interpolation if (best_x > 6 * xb || best_x < (-xb)) @@ -213,7 +192,7 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_ trial = false; xa = 0; - f_cal(move0, move, dim, xc); + f_cal(move0.data(), move, dim, xc); xc = xb + xc; xpt = xc; @@ -222,7 +201,7 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_ else { double xtemp, ftemp; - f_cal(move0, grad, dim, fc); + f_cal(move0.data(), grad, dim, fc); fmin = std::abs(fc); nbrent++; @@ -251,7 +230,7 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_ goto CG_begin; } - normalize(cg_gradn, cg_grad0, dim); + normalize(cg_gradn, cg_grad0.data(), dim); setup_move(move, cg_gradn, best_x); Lattice_Change_Basic::change_lattice(ucell, move, lat); diff --git a/source/source_relax/lattice_change_cg.h b/source/source_relax/lattice_change_cg.h index 68fe1b69509..76d8dc70f8d 100644 --- a/source/source_relax/lattice_change_cg.h +++ b/source/source_relax/lattice_change_cg.h @@ -3,21 +3,22 @@ #include "source_base/matrix.h" #include "source_cell/unitcell.h" +#include class Lattice_Change_CG { public: Lattice_Change_CG(); - ~Lattice_Change_CG(); + ~Lattice_Change_CG() = default; void allocate(void); void start(UnitCell &ucell, const ModuleBase::matrix &stress_in, const double &etot); private: - double * lat0 = nullptr; - double * grad0 = nullptr; - double * cg_grad0 = nullptr; - double * move0 = nullptr; + std::vector lat0; + std::vector grad0; + std::vector cg_grad0; + std::vector move0; double e0=0.0; // setup gradients. From 1b8d93897500f7f420df4445b1e723b2476b705a Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 8 Jun 2026 21:36:13 +0800 Subject: [PATCH 06/16] update --- source/source_relax/CMakeLists.txt | 3 -- source/source_relax/ions_move_basic.cpp | 6 ++- source/source_relax/ions_move_basic.h | 50 +++++++++---------------- source/source_relax/ions_move_sd.cpp | 4 +- 4 files changed, 24 insertions(+), 39 deletions(-) diff --git a/source/source_relax/CMakeLists.txt b/source/source_relax/CMakeLists.txt index 30415edd3b5..136523a620f 100644 --- a/source/source_relax/CMakeLists.txt +++ b/source/source_relax/CMakeLists.txt @@ -1,12 +1,9 @@ add_library( relax OBJECT - relax_driver.cpp - relax_sync.cpp line_search.cpp - bfgs.cpp lbfgs.cpp relax_nsync.cpp diff --git a/source/source_relax/ions_move_basic.cpp b/source/source_relax/ions_move_basic.cpp index e4a2caef9db..5eea695ceaf 100644 --- a/source/source_relax/ions_move_basic.cpp +++ b/source/source_relax/ions_move_basic.cpp @@ -5,6 +5,7 @@ #include "source_base/global_variable.h" #include "source_cell/update_cell.h" #include "source_cell/print_cell.h" + int Ions_Move_Basic::dim = 0; bool Ions_Move_Basic::converged = false; double Ions_Move_Basic::largest_grad = 0.0; @@ -96,9 +97,10 @@ void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos) const double move_threshold = 1.0e-10; const int total_freedom = ucell.nat * 3; - if (ModuleSymmetry::Symmetry::symm_flag && ucell.symm.all_mbl && ucell.symm.nrotk > 0) { + if (ModuleSymmetry::Symmetry::symm_flag && ucell.symm.all_mbl && ucell.symm.nrotk > 0) + { ucell.symm.symmetrize_vec3_nat(move); -} + } for (int i = 0; i < total_freedom; i++) { diff --git a/source/source_relax/ions_move_basic.h b/source/source_relax/ions_move_basic.h index 8f030850286..e26acbba2fa 100644 --- a/source/source_relax/ions_move_basic.h +++ b/source/source_relax/ions_move_basic.h @@ -6,56 +6,42 @@ namespace Ions_Move_Basic { -extern int dim; // dimension of the free variables, -extern bool converged; // converged force or not, -extern double largest_grad; // largest gradient among the forces, -extern int update_iter; // number of sucesfully updated iterations, -extern int istep; // index of ionic steps, -extern double ediff; // energy difference compared to last step, -extern double etot; // total energy of this step, -extern double etot_p; // total energy of last step, - -extern double trust_radius; // trust radius now, -extern double trust_radius_old; // old trust radius, -extern double relax_bfgs_rmax; // max value of trust radius, -extern double relax_bfgs_rmin; // min value of trust radius, -extern double relax_bfgs_init; // initial value of trust radius, -extern double best_xxx; // the last step length of cg , we use it as bfgs`s initial step length -extern std::vector relax_method; // relaxation method, +extern int dim; // dimension of the free variables +extern bool converged; // converged force or not +extern double largest_grad; // largest gradient among the forces +extern int update_iter; // number of successfully updated iterations +extern int istep; // index of ionic steps +extern double ediff; // energy difference compared to last step +extern double etot; // total energy of this step +extern double etot_p; // total energy of last step + +extern double trust_radius; // trust radius now +extern double trust_radius_old; // old trust radius +extern double relax_bfgs_rmax; // max value of trust radius +extern double relax_bfgs_rmin; // min value of trust radius +extern double relax_bfgs_init; // initial value of trust radius +extern double best_xxx; // the last step length of cg, we use it as bfgs initial step length +extern std::vector relax_method; // relaxation method extern int out_stru; // output the structure or not -// funny way to pass this parameter, but nevertheless -//---------------------------------------------------------------------------- // setup the gradient, all the same for any geometry optimization methods. -//---------------------------------------------------------------------------- void setup_gradient(const UnitCell &ucell, const ModuleBase::matrix &force, double *pos, double *grad); -//---------------------------------------------------------------------------- // move the atom positions, considering the periodic boundary condition. -//---------------------------------------------------------------------------- void move_atoms(UnitCell &ucell, double *move, double *pos); -//---------------------------------------------------------------------------- -// check the converged conditions ( if largest gradient is smaller than -// the threshold) -//---------------------------------------------------------------------------- +// check the converged conditions ( if largest gradient is smaller than the threshold ) void check_converged(const UnitCell &ucell, const double *grad); -//---------------------------------------------------------------------------- // terminate the geometry optimization. -//---------------------------------------------------------------------------- void terminate(const UnitCell &ucell); -//---------------------------------------------------------------------------- // setup the total energy, keep the new energy or not. -//---------------------------------------------------------------------------- void setup_etot(const double &energy_in, const bool judgement); double dot_func(const double *a, const double *b, const int &dim); -//---------------------------------------------------------------------------- -// third order interpolation scheme, -//---------------------------------------------------------------------------- +// third order interpolation scheme void third_order(); } // namespace Ions_Move_Basic diff --git a/source/source_relax/ions_move_sd.cpp b/source/source_relax/ions_move_sd.cpp index 133f809367c..8302e509a46 100644 --- a/source/source_relax/ions_move_sd.cpp +++ b/source/source_relax/ions_move_sd.cpp @@ -55,7 +55,7 @@ void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const } // normalize the gradient, in convinience to // move atom. - double norm = dot_func(grad_saved, grad_saved, dim); + double norm = dot_func(grad_saved.data(), grad_saved.data(), dim); norm = sqrt(norm); for (int i = 0; i < dim; i++) { @@ -75,7 +75,7 @@ void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const { move[i] = -grad_saved[i] * trust_radius; } - move_atoms(ucell, move.data(), pos_saved); + move_atoms(ucell, move.data(), pos_saved.data()); Ions_Move_Basic::update_iter++; } From 2730b795462782c03bdc4c10f86dd5f616010a4d Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Mon, 8 Jun 2026 21:53:57 +0800 Subject: [PATCH 07/16] move bfgs.h to ions_move_bfgs2.h --- source/Makefile.Objects | 2 +- source/source_relax/CMakeLists.txt | 2 +- .../{bfgs.cpp => ions_move_bfgs2.cpp} | 22 +++++++++---------- .../{bfgs.h => ions_move_bfgs2.h} | 8 +++---- source/source_relax/ions_move_methods.cpp | 2 +- source/source_relax/ions_move_methods.h | 4 ++-- source/source_relax/relax_driver.h | 4 ++-- source/source_relax/test/bfgs_test.cpp | 6 ++--- 8 files changed, 25 insertions(+), 25 deletions(-) rename source/source_relax/{bfgs.cpp => ions_move_bfgs2.cpp} (94%) rename source/source_relax/{bfgs.h => ions_move_bfgs2.h} (95%) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 6e1adc3f448..e1d14749d8c 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -464,6 +464,7 @@ OBJS_RELAXATION=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,7 +473,6 @@ OBJS_RELAXATION=bfgs_basic.o\ lattice_change_methods.o\ relax_nsync.o\ relax_sync.o\ - bfgs.o\ lbfgs.o\ matrix_methods.o\ line_search.o\ diff --git a/source/source_relax/CMakeLists.txt b/source/source_relax/CMakeLists.txt index 136523a620f..87aa8d1c3b3 100644 --- a/source/source_relax/CMakeLists.txt +++ b/source/source_relax/CMakeLists.txt @@ -4,12 +4,12 @@ add_library( relax_driver.cpp relax_sync.cpp line_search.cpp - bfgs.cpp 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 diff --git a/source/source_relax/bfgs.cpp b/source/source_relax/ions_move_bfgs2.cpp similarity index 94% rename from source/source_relax/bfgs.cpp rename to source/source_relax/ions_move_bfgs2.cpp index b7fbbe821a1..79fedd84db9 100644 --- a/source/source_relax/bfgs.cpp +++ b/source/source_relax/ions_move_bfgs2.cpp @@ -1,4 +1,4 @@ -#include "bfgs.h" +#include "ions_move_bfgs2.h" #include "source_base/module_external/lapack_connector.h" #include "source_io/module_parameter/parameter.h" #include "ions_move_basic.h" @@ -7,7 +7,7 @@ #include "source_cell/print_cell.h" // lanshuyue add 2025-06-19 //! initialize H0、H、pos0、force0、force -void BFGS::allocate(const int _size) +void Ions_Move_BFGS2::allocate(const int _size) { assert(_size > 0); alpha=70;//default value in ase is 70 @@ -34,7 +34,7 @@ void BFGS::allocate(const int _size) } -void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell) +void Ions_Move_BFGS2::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell) { if(!is_initialized) { @@ -84,7 +84,7 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell) unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,GlobalV::ofs_running); } -void BFGS::GetPos(UnitCell& ucell,std::vector>& pos) +void Ions_Move_BFGS2::GetPos(UnitCell& ucell,std::vector>& pos) { assert(pos.size() == ucell.nat); int k=0; @@ -100,7 +100,7 @@ void BFGS::GetPos(UnitCell& ucell,std::vector>& pos) } } -void BFGS::GetPostaud(UnitCell& ucell, +void Ions_Move_BFGS2::GetPostaud(UnitCell& ucell, std::vector>& pos_taud) { assert(pos_taud.size() == ucell.nat); @@ -117,7 +117,7 @@ void BFGS::GetPostaud(UnitCell& ucell, } } -void BFGS::PrepareStep(std::vector>& force, +void Ions_Move_BFGS2::PrepareStep(std::vector>& force, std::vector>& pos, std::vector>& H, std::vector& pos0, @@ -177,7 +177,7 @@ void BFGS::PrepareStep(std::vector>& force, force0 = ReshapeMToV(force); } -void BFGS::Update(std::vector& pos, +void Ions_Move_BFGS2::Update(std::vector& pos, std::vector& force, std::vector>& H, UnitCell& ucell) @@ -253,7 +253,7 @@ void BFGS::Update(std::vector& pos, H = MSubM(H, term4); } -void BFGS::DetermineStep(std::vector& steplength, +void Ions_Move_BFGS2::DetermineStep(std::vector& steplength, std::vector>& dpos, double& maxstep) { @@ -273,7 +273,7 @@ void BFGS::DetermineStep(std::vector& steplength, } } -void BFGS::UpdatePos(UnitCell& ucell) +void Ions_Move_BFGS2::UpdatePos(UnitCell& ucell) { double a[3*size]; for(int i=0;iMAddM(pos, dpos);*/ } -void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell) +void Ions_Move_BFGS2::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell) { std::vector grad= std::vector(3*size, 0.0); int iat = 0; @@ -361,7 +361,7 @@ void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell } } -void BFGS::IsRestrain() +void Ions_Move_BFGS2::IsRestrain() { Ions_Move_Basic::converged = largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A #include @@ -10,7 +10,7 @@ #include "source_cell/unitcell.h" -class BFGS +class Ions_Move_BFGS2 { public: void allocate(const int _size);//initialize parameters @@ -46,4 +46,4 @@ class BFGS }; -#endif // BFGS_H +#endif // IONS_MOVE_BFGS2_H diff --git a/source/source_relax/ions_move_methods.cpp b/source/source_relax/ions_move_methods.cpp index 15f27d67d33..db485f4c09e 100644 --- a/source/source_relax/ions_move_methods.cpp +++ b/source/source_relax/ions_move_methods.cpp @@ -31,7 +31,7 @@ void Ions_Move_Methods::allocate(const int &natom) else if (Ions_Move_Basic::relax_method[0] == "cg_bfgs") { this->cg.allocate(); - this->bfgs.allocate(); // added by pengfei 13-8-8 + this->bfgs.allocate(); } else if(Ions_Move_Basic::relax_method[0] == "bfgs"&&Ions_Move_Basic::relax_method[1] == "1") { diff --git a/source/source_relax/ions_move_methods.h b/source/source_relax/ions_move_methods.h index cdddc4fdd61..b9b4e28611e 100644 --- a/source/source_relax/ions_move_methods.h +++ b/source/source_relax/ions_move_methods.h @@ -5,7 +5,7 @@ #include "ions_move_bfgs.h" #include "ions_move_cg.h" #include "ions_move_sd.h" -#include "bfgs.h" +#include "ions_move_bfgs2.h" #include "lbfgs.h" class Ions_Move_Methods @@ -47,7 +47,7 @@ class Ions_Move_Methods Ions_Move_BFGS bfgs; Ions_Move_CG cg; Ions_Move_SD sd; - BFGS bfgs_trad; + Ions_Move_BFGS2 bfgs_trad; LBFGS lbfgs; }; #endif diff --git a/source/source_relax/relax_driver.h b/source/source_relax/relax_driver.h index 78b343efd81..2c970e3a451 100644 --- a/source/source_relax/relax_driver.h +++ b/source/source_relax/relax_driver.h @@ -6,7 +6,7 @@ #include "source_esolver/esolver_ks.h" #include "relax_sync.h" #include "relax_nsync.h" -#include "bfgs.h" +#include "ions_move_bfgs2.h" #include "source_io/module_parameter/input_parameter.h" class Relax_Driver @@ -32,7 +32,7 @@ class Relax_Driver // old relaxation method Relax_old rl_old; - BFGS bfgs_trad; + Ions_Move_BFGS2 bfgs_trad; }; diff --git a/source/source_relax/test/bfgs_test.cpp b/source/source_relax/test/bfgs_test.cpp index a90c6b4ed75..81852def15f 100644 --- a/source/source_relax/test/bfgs_test.cpp +++ b/source/source_relax/test/bfgs_test.cpp @@ -4,7 +4,7 @@ #include "for_test.h" #define private public -#include "source_relax/bfgs.h" +#include "source_relax/ions_move_bfgs2.h" #undef private #define private public @@ -14,12 +14,12 @@ #include "source_relax/ions_move_basic.h" // for Ions_Move_Basic static members /************************************************ - * unit tests for BFGS (no MockUnitCell) + * unit tests for Ions_Move_BFGS2 (no MockUnitCell) ***********************************************/ class BFGSTest : public ::testing::Test { protected: - BFGS bfgs; + Ions_Move_BFGS2 bfgs; void SetUp() override { // Initialize variables before each test From 0ca3ecb899b802bb4d0e7442cad9f79e8be560c3 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 9 Jun 2026 13:48:48 +0800 Subject: [PATCH 08/16] fix --- source/source_relax/test/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_relax/test/CMakeLists.txt b/source/source_relax/test/CMakeLists.txt index 777e24d7cf9..2c19c41ff48 100644 --- a/source/source_relax/test/CMakeLists.txt +++ b/source/source_relax/test/CMakeLists.txt @@ -61,7 +61,7 @@ AddTest( AddTest( TARGET MODULE_RELAX_bfgs_test LIBS parameter ${math_libs} base device - SOURCES bfgs_test.cpp ../bfgs.cpp ../ions_move_basic.cpp ../matrix_methods.cpp ${cell_source_files} + SOURCES bfgs_test.cpp ../ions_move_bfgs2.cpp ../ions_move_basic.cpp ../matrix_methods.cpp ${cell_source_files} ) AddTest( From bdb057a07e7ac938eb2233300733057e1bfa88a1 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 9 Jun 2026 14:00:42 +0800 Subject: [PATCH 09/16] add relax_data --- source/Makefile.Objects | 3 +- source/source_relax/CMakeLists.txt | 1 + source/source_relax/relax_data.cpp | 18 +++++++++++ source/source_relax/relax_data.h | 51 ++++++++++++++++++++++++++++++ 4 files changed, 72 insertions(+), 1 deletion(-) create mode 100644 source/source_relax/relax_data.cpp create mode 100644 source/source_relax/relax_data.h diff --git a/source/Makefile.Objects b/source/Makefile.Objects index e1d14749d8c..e236a79fe2a 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -460,7 +460,8 @@ OBJS_PW=fft_bundle.o\ pw_transform.o\ pw_transform_k.o\ -OBJS_RELAXATION=bfgs_basic.o\ +OBJS_RELAXATION=relax_data.o\ + bfgs_basic.o\ relax_driver.o\ ions_move_basic.o\ ions_move_bfgs.o\ diff --git a/source/source_relax/CMakeLists.txt b/source/source_relax/CMakeLists.txt index 87aa8d1c3b3..b98170befdf 100644 --- a/source/source_relax/CMakeLists.txt +++ b/source/source_relax/CMakeLists.txt @@ -1,6 +1,7 @@ add_library( relax OBJECT + relax_data.cpp relax_driver.cpp relax_sync.cpp line_search.cpp diff --git a/source/source_relax/relax_data.cpp b/source/source_relax/relax_data.cpp new file mode 100644 index 00000000000..54c4e1a118f --- /dev/null +++ b/source/source_relax/relax_data.cpp @@ -0,0 +1,18 @@ +#include "relax_data.h" + +int Relax_Data::dim = 0; +bool Relax_Data::converged = false; +double Relax_Data::largest_grad = 0.0; +int Relax_Data::istep = 0; +double Relax_Data::ediff = 0.0; +double Relax_Data::etot = 0.0; +double Relax_Data::etot_p = 0.0; + +void Relax_Data::allocate(const int dim_in) { + pos.resize(dim_in, 0.0); + grad.resize(dim_in, 0.0); + move.resize(dim_in, 0.0); + pos_p.resize(dim_in, 0.0); + grad_p.resize(dim_in, 0.0); + move_p.resize(dim_in, 0.0); +} \ No newline at end of file diff --git a/source/source_relax/relax_data.h b/source/source_relax/relax_data.h new file mode 100644 index 00000000000..23caa55606e --- /dev/null +++ b/source/source_relax/relax_data.h @@ -0,0 +1,51 @@ +#ifndef RELAX_DATA_H +#define RELAX_DATA_H + +#include + +/** + * @brief Unified data structure for geometry optimization algorithms. + * + * This class provides a common data container for all relaxation methods + * (BFGS, CG, SD, etc.) to share and reuse key variables like positions, + * gradients, and move vectors. It serves as the single source of truth + * for optimization state across different algorithms. + */ +class Relax_Data { +public: + /** + * @brief Default constructor. + */ + Relax_Data() = default; + + /** + * @brief Default destructor. + */ + ~Relax_Data() = default; + + /** + * @brief Allocate memory for data vectors. + * @param dim_in Dimension of the optimization problem (3 * number of atoms). + */ + void allocate(const int dim_in); + + // Static members - shared global state across all relaxation instances + static int dim; ///< Dimension of free variables (3 * number of atoms) + static bool converged; ///< Convergence flag: true if optimization is converged + static double largest_grad; ///< Largest gradient component (force) in current step + static int istep; ///< Current iteration step index + static double ediff; ///< Energy difference from previous step (etot - etot_p) + static double etot; ///< Total energy of current step + static double etot_p; ///< Total energy of previous step + + // Instance members - per-iteration data + std::vector pos; ///< Current atomic positions in Bohr + std::vector grad; ///< Current gradient (negative force) in Ry/Bohr + std::vector move; ///< Displacement vector for next step + + std::vector pos_p; ///< Previous atomic positions + std::vector grad_p; ///< Previous gradient + std::vector move_p; ///< Previous displacement +}; + +#endif \ No newline at end of file From 73edf35c57a54c0e8e7252d01f541434b6cad507 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 9 Jun 2026 14:10:15 +0800 Subject: [PATCH 10/16] update --- source/source_relax/ions_move_basic.cpp | 11 +-- source/source_relax/ions_move_basic.h | 95 ++++++++++++++------ source/source_relax/lattice_change_basic.cpp | 13 +-- source/source_relax/lattice_change_basic.h | 82 ++++++++++------- 4 files changed, 124 insertions(+), 77 deletions(-) diff --git a/source/source_relax/ions_move_basic.cpp b/source/source_relax/ions_move_basic.cpp index 5eea695ceaf..c2b3d0cd9c8 100644 --- a/source/source_relax/ions_move_basic.cpp +++ b/source/source_relax/ions_move_basic.cpp @@ -6,23 +6,14 @@ #include "source_cell/update_cell.h" #include "source_cell/print_cell.h" -int Ions_Move_Basic::dim = 0; -bool Ions_Move_Basic::converged = false; -double Ions_Move_Basic::largest_grad = 0.0; +// Ions-specific parameters (shared variables are in Relax_Data) int Ions_Move_Basic::update_iter = 0; -int Ions_Move_Basic::istep = 0; - -double Ions_Move_Basic::ediff = 0.0; -double Ions_Move_Basic::etot = 0.0; -double Ions_Move_Basic::etot_p = 0.0; - double Ions_Move_Basic::trust_radius = 0.0; double Ions_Move_Basic::trust_radius_old = 0.0; double Ions_Move_Basic::relax_bfgs_rmax = -1.0; // default is 0.8 double Ions_Move_Basic::relax_bfgs_rmin = -1.0; // default is 1e-5 double Ions_Move_Basic::relax_bfgs_init = -1.0; // default is 0.5 double Ions_Move_Basic::best_xxx = 1.0; - int Ions_Move_Basic::out_stru = 0; std::vector Ions_Move_Basic::relax_method = {"bfgs","2"}; diff --git a/source/source_relax/ions_move_basic.h b/source/source_relax/ions_move_basic.h index e26acbba2fa..c852881a9af 100644 --- a/source/source_relax/ions_move_basic.h +++ b/source/source_relax/ions_move_basic.h @@ -1,48 +1,91 @@ #ifndef IONS_MOVE_BASIC_H #define IONS_MOVE_BASIC_H +#include "relax_data.h" #include "source_base/matrix.h" #include "source_cell/unitcell.h" +/** + * @namespace Ions_Move_Basic + * @brief Basic utilities and shared state for ionic relaxation algorithms. + * + * This namespace provides common functions and parameters used by all + * ion movement methods (BFGS, CG, SD, etc.). It shares core state variables + * through references to Relax_Data, ensuring consistent data across different + * optimization algorithms. + */ namespace Ions_Move_Basic { -extern int dim; // dimension of the free variables -extern bool converged; // converged force or not -extern double largest_grad; // largest gradient among the forces -extern int update_iter; // number of successfully updated iterations -extern int istep; // index of ionic steps -extern double ediff; // energy difference compared to last step -extern double etot; // total energy of this step -extern double etot_p; // total energy of last step - -extern double trust_radius; // trust radius now -extern double trust_radius_old; // old trust radius -extern double relax_bfgs_rmax; // max value of trust radius -extern double relax_bfgs_rmin; // min value of trust radius -extern double relax_bfgs_init; // initial value of trust radius -extern double best_xxx; // the last step length of cg, we use it as bfgs initial step length -extern std::vector relax_method; // relaxation method -extern int out_stru; // output the structure or not - -// setup the gradient, all the same for any geometry optimization methods. +// Shared state variables (referenced from Relax_Data for unified data sharing) +static int& dim = Relax_Data::dim; ///< Dimension of free variables (3 * number of atoms) +static bool& converged = Relax_Data::converged; ///< Convergence flag +static double& largest_grad = Relax_Data::largest_grad; ///< Largest gradient component +static int& istep = Relax_Data::istep; ///< Current ionic step index +static double& ediff = Relax_Data::ediff; ///< Energy difference from previous step +static double& etot = Relax_Data::etot; ///< Total energy of current step +static double& etot_p = Relax_Data::etot_p; ///< Total energy of previous step + +// Ions-specific parameters (not shared with lattice change) +extern int update_iter; ///< Number of successfully updated iterations +extern double trust_radius; ///< Current trust radius +extern double trust_radius_old; ///< Previous trust radius +extern double relax_bfgs_rmax; ///< Maximum trust radius (default: 0.8 Bohr) +extern double relax_bfgs_rmin; ///< Minimum trust radius (default: 1e-5 Bohr) +extern double relax_bfgs_init; ///< Initial trust radius (default: 0.5 Bohr) +extern double best_xxx; ///< Last step length from CG, used as BFGS initial guess +extern std::vector relax_method; ///< Relaxation method settings +extern int out_stru; ///< Structure output flag + +/** + * @brief Setup gradient from atomic forces. + * @param ucell Unit cell containing atomic information + * @param force Force matrix (nat x 3) + * @param pos Output position array (dimension: dim) + * @param grad Output gradient array (dimension: dim) + */ void setup_gradient(const UnitCell &ucell, const ModuleBase::matrix &force, double *pos, double *grad); -// move the atom positions, considering the periodic boundary condition. +/** + * @brief Move atoms according to displacement vector. + * @param ucell Unit cell to update + * @param move Displacement vector (dimension: dim) + * @param pos Current position array (dimension: dim) + */ void move_atoms(UnitCell &ucell, double *move, double *pos); -// check the converged conditions ( if largest gradient is smaller than the threshold ) +/** + * @brief Check convergence based on gradient threshold. + * @param ucell Unit cell containing lattice information + * @param grad Gradient array (dimension: dim) + */ void check_converged(const UnitCell &ucell, const double *grad); -// terminate the geometry optimization. +/** + * @brief Terminate geometry optimization and output results. + * @param ucell Unit cell to output + */ void terminate(const UnitCell &ucell); -// setup the total energy, keep the new energy or not. +/** + * @brief Update energy values and compute energy difference. + * @param energy_in Input energy value + * @param judgement Flag for SD method (true) or BFGS (false) + */ void setup_etot(const double &energy_in, const bool judgement); -double dot_func(const double *a, const double *b, const int &dim); +/** + * @brief Compute dot product of two vectors. + * @param a First vector + * @param b Second vector + * @param dim_in Dimension of vectors + * @return Dot product value + */ +double dot_func(const double *a, const double *b, const int &dim_in); -// third order interpolation scheme +/** + * @brief Third-order polynomial interpolation for line search. + */ void third_order(); } // namespace Ions_Move_Basic -#endif +#endif \ No newline at end of file diff --git a/source/source_relax/lattice_change_basic.cpp b/source/source_relax/lattice_change_basic.cpp index 1bc41ce6238..6c2ca4b68f2 100644 --- a/source/source_relax/lattice_change_basic.cpp +++ b/source/source_relax/lattice_change_basic.cpp @@ -6,19 +6,10 @@ #include "source_io/module_parameter/parameter.h" #include "source_cell/update_cell.h" -int Lattice_Change_Basic::dim = 0; -bool Lattice_Change_Basic::converged = true; -double Lattice_Change_Basic::largest_grad = 0.0; +// Lattice-specific parameters (shared variables are in Relax_Data) int Lattice_Change_Basic::update_iter = 0; -int Lattice_Change_Basic::istep = 0; int Lattice_Change_Basic::stress_step = 0; - -double Lattice_Change_Basic::ediff = 0.0; -double Lattice_Change_Basic::etot = 0.0; -double Lattice_Change_Basic::etot_p = 0.0; - -// double Lattice_Change_Basic::lattice_change_ini = 0.5; // default is 0.5 -double Lattice_Change_Basic::lattice_change_ini = 0.01; // default is 0.5 +double Lattice_Change_Basic::lattice_change_ini = 0.01; // default is 0.01 std::string Lattice_Change_Basic::fixed_axes = "None"; void Lattice_Change_Basic::setup_gradient(const UnitCell &ucell, double *lat, double *grad, ModuleBase::matrix &stress) diff --git a/source/source_relax/lattice_change_basic.h b/source/source_relax/lattice_change_basic.h index 263c7bd5ca0..db5d6de1a35 100644 --- a/source/source_relax/lattice_change_basic.h +++ b/source/source_relax/lattice_change_basic.h @@ -1,48 +1,70 @@ #ifndef LATTICE_CHANGE_BASIC_H #define LATTICE_CHANGE_BASIC_H +#include "relax_data.h" #include "source_base/matrix.h" #include "source_cell/unitcell.h" +/** + * @namespace Lattice_Change_Basic + * @brief Basic utilities and shared state for lattice relaxation algorithms. + * + * This namespace provides common functions and parameters used by lattice + * optimization methods. It shares core state variables through references to + * Relax_Data, ensuring consistent data across different optimization algorithms. + */ namespace Lattice_Change_Basic { -extern int dim; // dimension of the free variables, -extern bool converged; // converged force or not, -extern double largest_grad; // largest gradient among the forces, -extern int update_iter; // number of sucesfully updated iterations, -extern int istep; // index of ionic steps, -extern int stress_step; // index of stress step -extern double ediff; // energy difference compared to last step, -extern double etot; // total energy of this step, -extern double etot_p; // total energy of last step, - -extern double lattice_change_ini; // initial value of trust radius, -extern std::string fixed_axes; // convert from INPUT.fixed_axes - -//---------------------------------------------------------------------------- -// setup the gradient, all the same for any geometry optimization methods. -//---------------------------------------------------------------------------- +// Shared state variables (referenced from Relax_Data for unified data sharing) +static int& dim = Relax_Data::dim; ///< Dimension of free variables (9 for full lattice) +static bool& converged = Relax_Data::converged; ///< Convergence flag +static double& largest_grad = Relax_Data::largest_grad; ///< Largest gradient component +static int& istep = Relax_Data::istep; ///< Current ionic step index +static double& ediff = Relax_Data::ediff; ///< Energy difference from previous step +static double& etot = Relax_Data::etot; ///< Total energy of current step +static double& etot_p = Relax_Data::etot_p; ///< Total energy of previous step + +// Lattice-specific parameters (not shared with ion movement) +extern int update_iter; ///< Number of successfully updated iterations +extern int stress_step; ///< Index of stress optimization step +extern double lattice_change_ini; ///< Initial trust radius for lattice change (default: 0.01) +extern std::string fixed_axes; ///< Fixed axes constraint ("None", "shape", "volume", or specific axes) + +/** + * @brief Setup gradient from stress tensor. + * @param ucell Unit cell containing lattice information + * @param lat Output lattice vector array (9 elements: e11, e12, e13, e21, e22, e23, e31, e32, e33) + * @param grad Output gradient array (9 elements) + * @param stress Stress tensor (3x3 matrix) + */ void setup_gradient(const UnitCell &ucell, double *lat, double *grad, ModuleBase::matrix &stress); -//---------------------------------------------------------------------------- -// move the atom positions, considering the periodic boundary condition. -//---------------------------------------------------------------------------- +/** + * @brief Update lattice vectors according to displacement. + * @param ucell Unit cell to update + * @param move Displacement vector for lattice change (9 elements) + * @param lat Current lattice vectors (9 elements) + */ void change_lattice(UnitCell &ucell, double *move, double *lat); -//---------------------------------------------------------------------------- -// check the converged conditions ( if largest gradient is smaller than -// the threshold) -//---------------------------------------------------------------------------- +/** + * @brief Check convergence based on stress threshold. + * @param ucell Unit cell containing lattice constraints + * @param stress Stress tensor (3x3 matrix) + * @param grad Gradient array (9 elements) + */ void check_converged(const UnitCell &ucell, ModuleBase::matrix &stress, double *grad); -//---------------------------------------------------------------------------- -// terminate the geometry optimization. -//---------------------------------------------------------------------------- +/** + * @brief Terminate lattice optimization and output results. + */ void terminate(void); -//---------------------------------------------------------------------------- -// setup the total energy, keep the new energy or not. -//---------------------------------------------------------------------------- +/** + * @brief Update energy values and compute energy difference. + * @param energy_in Input energy value + * @param judgement Flag for SD method (true) or BFGS (false) + */ void setup_etot(const double &energy_in, const bool judgement); } // namespace Lattice_Change_Basic -#endif +#endif \ No newline at end of file From 44eb6410a9eec19767ca6e8599e98f6b85956341 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 9 Jun 2026 14:22:12 +0800 Subject: [PATCH 11/16] update --- source/Makefile.Objects | 1 + source/source_relax/CMakeLists.txt | 1 + source/source_relax/cg_base.cpp | 213 +++++++++++++++++++++++++++++ source/source_relax/cg_base.h | 83 +++++++++++ 4 files changed, 298 insertions(+) create mode 100644 source/source_relax/cg_base.cpp create mode 100644 source/source_relax/cg_base.h diff --git a/source/Makefile.Objects b/source/Makefile.Objects index e236a79fe2a..fc49ec2227a 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -461,6 +461,7 @@ OBJS_PW=fft_bundle.o\ pw_transform_k.o\ OBJS_RELAXATION=relax_data.o\ + cg_base.o\ bfgs_basic.o\ relax_driver.o\ ions_move_basic.o\ diff --git a/source/source_relax/CMakeLists.txt b/source/source_relax/CMakeLists.txt index b98170befdf..692e9315fb9 100644 --- a/source/source_relax/CMakeLists.txt +++ b/source/source_relax/CMakeLists.txt @@ -2,6 +2,7 @@ add_library( relax OBJECT relax_data.cpp + cg_base.cpp relax_driver.cpp relax_sync.cpp line_search.cpp diff --git a/source/source_relax/cg_base.cpp b/source/source_relax/cg_base.cpp new file mode 100644 index 00000000000..9fee5024d95 --- /dev/null +++ b/source/source_relax/cg_base.cpp @@ -0,0 +1,213 @@ +#include "cg_base.h" +#include +#include "source_base/global_function.h" + +void CG_Base::setup_cg_grad(const int dim, double *grad, const double *grad0, + double *cg_grad, const double *cg_grad0, + const int &ncggrad, int &flag) +{ + ModuleBase::TITLE("CG_Base", "setup_cg_grad"); + double gamma = 0.0; + + if (ncggrad % 10000 == 0 || flag == 2) + { + for (int i = 0; i < dim; i++) + { + cg_grad[i] = grad[i]; + } + } + else + { + double gp_gp = 0.0; // grad_p.grad_p + double gg = 0.0; // grad.grad + double g_gp = 0.0; // grad_p.grad + double cgp_gp = 0.0; // cg_grad_p.grad_p + double cgp_g = 0.0; // cg_grad_p.grad + for (int i = 0; i < dim; i++) + { + gp_gp += grad0[i] * grad0[i]; + gg += grad[i] * grad[i]; + g_gp += grad0[i] * grad[i]; + cgp_gp += cg_grad0[i] * grad0[i]; + cgp_g += cg_grad0[i] * grad[i]; + } + + const double gamma1 = gg / gp_gp; // FR + const double gamma2 = (gg - g_gp) / gp_gp; // PRP + + if (gamma1 < 0.5) + { + gamma = gamma1; + } + else + { + gamma = gamma2; + } + + for (int i = 0; i < dim; i++) + { + cg_grad[i] = grad[i] + gamma * cg_grad0[i]; + } + } +} + +void CG_Base::setup_move(const int dim, double *move, double *cg_gradn, const double &trust_radius) +{ + ModuleBase::TITLE("CG_Base", "setup_move"); + for (int i = 0; i < dim; ++i) + { + move[i] = -cg_gradn[i] * trust_radius; + } +} + +void CG_Base::Brent(double &fa, double &fb, double &fc, + double &xa, double &xb, double &xc, + double &best_x, double &xpt) +{ + ModuleBase::TITLE("CG_Base", "Brent"); + double dmove = 0.0; + double tmp = 0.0; + double k2 = 0.0; + double k1 = 0.0; + double k0 = 0.0; + double xnew1 = 0.0; + double xnew2 = 0.0; + double ecalnew1 = 0.0; + double ecalnew2 = 0.0; + + if ((fa * fb) > 0) + { + dmove = (xc * fa - xa * fc) / (fa - fc); + if (dmove > 4 * xc) + { + dmove = 4 * xc; + } + xb = xc; + fb = fc; + } + else + { + k2 = -((fb - fc) / (xb - xc) - (fa - fc) / (xa - xc)) / (xa - xb); + k1 = (fa - fc) / (xa - xc) - k2 * (xa + xc); + k0 = fa - k1 * xa - k2 * xa * xa; + xnew1 = (-k1 - sqrt(k1 * k1 - 4 * k2 * k0)) / (2 * k2); + xnew2 = (-k1 + sqrt(k1 * k1 - 4 * k2 * k0)) / (2 * k2); + + if (xnew1 > xnew2) + { + tmp = xnew2; + xnew2 = xnew1; + xnew1 = tmp; + } + + ecalnew1 = k2 * xnew1 * xnew1 * xnew1 / 3 + k1 * xnew1 * xnew1 / 2 + k0 * xnew1; + ecalnew2 = k2 * xnew2 * xnew2 * xnew2 / 3 + k1 * xnew2 * xnew2 / 2 + k0 * xnew2; + dmove = xnew1; + + if (ecalnew1 > ecalnew2) + { + dmove = xnew2; + } + if (dmove < 0) + { + dmove = 2 * xc; + } + if (fa * fc > 0) + { + xa = xc; + fa = fc; + } + if (fb * fc > 0) + { + xb = xc; + fb = fc; + } + } + + best_x = dmove - xpt; + xpt = dmove; + xc = dmove; +} + +void CG_Base::f_cal(const int dim, const double *g0, const double *g1, double &f_value) +{ + ModuleBase::TITLE("CG_Base", "f_cal"); + double hv0 = 0.0; + double hel = 0.0; + for (int i = 0; i < dim; i++) + { + hel += g0[i] * g1[i]; + hv0 += g0[i] * g0[i]; + } + f_value = hel / sqrt(hv0); +} + +void CG_Base::third_order(const double &e0, const double &e1, + const double &fa, const double &fb, + const double x, double &best_x) +{ + ModuleBase::TITLE("CG_Base", "third_order"); + double k3 = 0.0; + double k2 = 0.0; + double k1 = 0.0; + double dmoveh = 0.0; + double dmove1 = 0.0; + double dmove2 = 0.0; + double dmove = 0.0; + double ecal1 = 0.0; + double ecal2 = 0.0; + + k3 = 3 * ((fb + fa) * x - 2 * (e1 - e0)) / (x * x * x); + k2 = (fb - fa) / x - k3 * x; + k1 = fa; + + dmoveh = x * fb / (fa - fb); + dmove1 = -k2 * (1 - sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); + dmove2 = -k2 * (1 + sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); + + if ((std::abs(k3 / k1) < 0.01) || ((k1 * k3 / (k2 * k2)) >= 0.25)) + { + dmove = dmoveh; + } + else + { + dmove1 = -k2 * (1 - sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); + dmove2 = -k2 * (1 + sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); + ecal1 = k3 * dmove1 * dmove1 * dmove1 / 3 + k2 * dmove1 * dmove1 / 2 + k1 * dmove1; + ecal2 = k3 * dmove2 * dmove2 * dmove2 / 3 + k2 * dmove2 * dmove2 / 2 + k1 * dmove2; + if (ecal2 > ecal1) + { + dmove = dmove1 - x; + } + else + { + dmove = dmove2 - x; + } + + if (k3 < 0) + { + dmove = dmoveh; + } + } + + best_x = dmove; +} + +void CG_Base::normalize(const int dim, double *cg_gradn, const double *cg_grad) +{ + ModuleBase::TITLE("CG_Base", "normalize"); + double norm = 0.0; + for (int i = 0; i < dim; ++i) + { + norm += pow(cg_grad[i], 2); + } + norm = sqrt(norm); + + if (norm != 0.0) + { + for (int i = 0; i < dim; ++i) + { + cg_gradn[i] = cg_grad[i] / norm; + } + } +} \ No newline at end of file diff --git a/source/source_relax/cg_base.h b/source/source_relax/cg_base.h new file mode 100644 index 00000000000..ce245830ce4 --- /dev/null +++ b/source/source_relax/cg_base.h @@ -0,0 +1,83 @@ +#ifndef CG_BASE_H +#define CG_BASE_H + +/** + * @class CG_Base + * @brief Base class for conjugate gradient optimization algorithms. + * + * This class provides common CG utility functions that can be shared between + * different optimization implementations (ions movement and lattice change). + * The functions are designed to be independent of specific data sources. + */ +class CG_Base { +protected: + /** + * @brief Setup conjugate gradient direction. + * @param dim Dimension of the optimization problem + * @param grad Current gradient + * @param grad0 Previous gradient + * @param cg_grad Output CG direction + * @param cg_grad0 Previous CG direction + * @param ncggrad CG iteration counter + * @param flag Restart flag + */ + void setup_cg_grad(const int dim, double *grad, const double *grad0, + double *cg_grad, const double *cg_grad0, + const int &ncggrad, int &flag); + + /** + * @brief Calculate movement vector from CG direction. + * @param dim Dimension of the optimization problem + * @param move Output movement vector + * @param cg_gradn Normalized CG direction + * @param trust_radius Step size limit + */ + void setup_move(const int dim, double *move, double *cg_gradn, const double &trust_radius); + + /** + * @brief Brent's method for one-dimensional minimization. + * @param fa Function value at xa + * @param fb Function value at xb + * @param fc Function value at xc + * @param xa Left bound + * @param xb Middle point + * @param xc Right bound + * @param best_x Output best step length + * @param xpt Previous best point + */ + void Brent(double &fa, double &fb, double &fc, + double &xa, double &xb, double &xc, + double &best_x, double &xpt); + + /** + * @brief Calculate projection of gradient onto search direction. + * @param dim Dimension of vectors + * @param g0 First vector (usually search direction) + * @param g1 Second vector (usually gradient) + * @param f_value Output projection value + */ + void f_cal(const int dim, const double *g0, const double *g1, double &f_value); + + /** + * @brief Third-order polynomial interpolation for line search. + * @param e0 Energy at previous step + * @param e1 Energy at current step + * @param fa Gradient projection at previous step + * @param fb Gradient projection at current step + * @param x Current step length + * @param best_x Output optimal step length + */ + void third_order(const double &e0, const double &e1, + const double &fa, const double &fb, + const double x, double &best_x); + + /** + * @brief Normalize vector. + * @param dim Dimension of vector + * @param cg_gradn Output normalized vector + * @param cg_grad Input vector + */ + void normalize(const int dim, double *cg_gradn, const double *cg_grad); +}; + +#endif \ No newline at end of file From 499adf24bddc321c0be9f8df133ed75b10234bb9 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 9 Jun 2026 14:34:12 +0800 Subject: [PATCH 12/16] update --- source/source_relax/ions_move_cg.cpp | 488 +++++----------------- source/source_relax/ions_move_cg.h | 40 +- source/source_relax/lattice_change_cg.cpp | 473 +++++---------------- source/source_relax/lattice_change_cg.h | 31 +- 4 files changed, 217 insertions(+), 815 deletions(-) diff --git a/source/source_relax/ions_move_cg.cpp b/source/source_relax/ions_move_cg.cpp index 05e18c4b7f2..00a73f648d6 100644 --- a/source/source_relax/ions_move_cg.cpp +++ b/source/source_relax/ions_move_cg.cpp @@ -2,59 +2,26 @@ #include "ions_move_basic.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" +#include using namespace Ions_Move_Basic; double Ions_Move_CG::RELAX_CG_THR = -1.0; // default is 0.5 -//=================== NOTES ======================== -// in vasp, it's like -// contolled by POTIM. -// the approximate minimum of the total energy -// is calculated from a cubic (or quadratic) -// interpolation taking into account the change -// of the total energy and the changes of forces. -// ( 3 pieces of information ) -// 1. trial step -// 2. corrector step using cubic or quadratic interpolation -// 3. Check the force and see if we need brent method interpolation. -// 4. a new trial step.... -// Brent's method is a complicated but popular -// root-finding algorithm combining the bisection method, -// the secant method and inverse quadratic interpolation -//=================== NOTES ======================== Ions_Move_CG::Ions_Move_CG() { - this->pos0 = nullptr; - this->grad0 = nullptr; - this->cg_grad0 = nullptr; - this->move0 = nullptr; -} - -Ions_Move_CG::~Ions_Move_CG() -{ - delete[] pos0; - delete[] grad0; - delete[] cg_grad0; - delete[] move0; + this->e0 = 0.0; } void Ions_Move_CG::allocate(void) { ModuleBase::TITLE("Ions_Move_CG", "allocate"); assert(dim > 0); - delete[] pos0; - delete[] grad0; - delete[] cg_grad0; - delete[] move0; - this->pos0 = new double[dim]; - this->grad0 = new double[dim]; - this->cg_grad0 = new double[dim]; - this->move0 = new double[dim]; - ModuleBase::GlobalFunc::ZEROS(pos0, dim); - ModuleBase::GlobalFunc::ZEROS(grad0, dim); - ModuleBase::GlobalFunc::ZEROS(cg_grad0, dim); - ModuleBase::GlobalFunc::ZEROS(move0, dim); + + this->pos0.resize(dim, 0.0); + this->grad0.resize(dim, 0.0); + this->cg_grad0.resize(dim, 0.0); + this->move0.resize(dim, 0.0); this->e0 = 0.0; } @@ -62,20 +29,9 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const { ModuleBase::TITLE("Ions_Move_CG", "start"); assert(dim > 0); - assert(pos0 != 0); - assert(grad0 != 0); - assert(cg_grad0 != 0); - assert(move0 != 0); - // sd , trial are two parameters, when sd=trial=true ,a new direction begins, when sd = false trial =true static bool sd = false; - - // a cubic interpolation is used to make the third point , - // when sa = trial = false , we use Brent to get the minimum point in this direction. static bool trial = false; - - // ncggrad is a parameter to control the cg method , every ten cg directions, - // we change the direction back to the steepest descent method static int ncggrad = 0; static double fa = 0.0; static double fb = 0.0; @@ -87,87 +43,73 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const static double steplength = 0.0; static double fmax = 0.0; static int nbrent = 0; - - // some arrays - double *pos = new double[dim]; - double *grad = new double[dim]; - double *cg_gradn = new double[dim]; - double *move = new double[dim]; - double *cg_grad = new double[dim]; + + std::vector pos(dim, 0.0); + std::vector grad(dim, 0.0); + std::vector cg_gradn(dim, 0.0); + std::vector move(dim, 0.0); + std::vector cg_grad(dim, 0.0); double best_x = 0.0; double fmin = 0.0; int flag = 0; - ModuleBase::GlobalFunc::ZEROS(pos, dim); - ModuleBase::GlobalFunc::ZEROS(grad, dim); - ModuleBase::GlobalFunc::ZEROS(move, dim); - ModuleBase::GlobalFunc::ZEROS(cg_grad, dim); - -CG_begin: - - if (Ions_Move_Basic::istep == 1) + while (true) { - steplength = Ions_Move_Basic::relax_bfgs_init; // read in the init trust radius - sd = true; - trial = true; - ncggrad = 0; - fa = 0.0; - fb = 0.0; - fc = 0.0; - xa = 0.0; - xb = 0.0; - xc = 0.0; - xpt = 0.0; - fmax = 0.0; - nbrent = 0; - } + if (Ions_Move_Basic::istep == 1) + { + steplength = Ions_Move_Basic::relax_bfgs_init; + sd = true; + trial = true; + ncggrad = 0; + fa = 0.0; + fb = 0.0; + fc = 0.0; + xa = 0.0; + xb = 0.0; + xc = 0.0; + xpt = 0.0; + fmax = 0.0; + nbrent = 0; + } - Ions_Move_Basic::setup_gradient(ucell, force, pos, grad); - // use energy_in and istep to setup etot and etot_old. - Ions_Move_Basic::setup_etot(etot_in, 0); - // use gradient and etot and etot_old to check - // if the result is converged. + Ions_Move_Basic::setup_gradient(ucell, force, pos.data(), grad.data()); + Ions_Move_Basic::setup_etot(etot_in, 0); + + if (flag == 0) + { + Ions_Move_Basic::check_converged(ucell, grad.data()); + } + if (Ions_Move_Basic::converged) + { + Ions_Move_Basic::terminate(ucell); + break; + } - if (flag == 0) - { - Ions_Move_Basic::check_converged(ucell, grad); - } - if (Ions_Move_Basic::converged) - { - Ions_Move_Basic::terminate(ucell); - } - else - { if (sd) { e0 = etot_in; - setup_cg_grad(grad, - grad0, - cg_grad, - cg_grad0, - ncggrad, - flag); // we use the last direction ,the last grad and the grad now to get the direction now + CG_Base::setup_cg_grad(dim, grad.data(), grad0.data(), cg_grad.data(), cg_grad0.data(), ncggrad, flag); ncggrad++; - normalize(cg_gradn, cg_grad, dim); - setup_move(move0, cg_gradn, steplength); // move the atom position - Ions_Move_Basic::move_atoms(ucell, move0, pos); + CG_Base::normalize(dim, cg_gradn.data(), cg_grad.data()); + CG_Base::setup_move(dim, move0.data(), cg_gradn.data(), steplength); + Ions_Move_Basic::move_atoms(ucell, move0.data(), pos.data()); - for (int i = 0; i < dim; i++) // grad0 ,cg_grad0 are used to store the grad and cg_grad for the future using + for (int i = 0; i < dim; i++) { grad0[i] = grad[i]; cg_grad0[i] = cg_grad[i]; } - f_cal(move0, move0, dim, xb); // xb = trial steplength - f_cal(move0, grad, dim, fa); // fa is the projection force in this direction + CG_Base::f_cal(dim, move0.data(), move0.data(), xb); + CG_Base::f_cal(dim, move0.data(), grad.data(), fa); fmax = fa; sd = false; if (Ions_Move_Basic::relax_method[0] == "cg_bfgs") { if (Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A - < RELAX_CG_THR) // cg to bfgs by pengfei 13-8-8 + < RELAX_CG_THR) { Ions_Move_Basic::relax_method[0] = "bfgs"; Ions_Move_Basic::relax_method[1] = "1"; @@ -176,301 +118,75 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const } Ions_Move_Basic::relax_bfgs_init = xb; + break; } - else - { - if (trial) - { - double e1 = etot_in; - f_cal(move0, grad, dim, fb); - f_cal(move0, move0, dim, xb); - if ((std::abs(fb) < std::abs((fa) / 10.0))) - { - sd = true; - trial = true; - steplength = xb; - flag = 1; - goto CG_begin; - } - normalize(cg_gradn, cg_grad0, dim); - third_order(e0, e1, fa, fb, xb, best_x); // cubic interpolation - - if (best_x > 6 * xb || best_x < (-xb)) - { - best_x = 6 * xb; - } + if (trial) + { + double e1 = etot_in; + CG_Base::f_cal(dim, move0.data(), grad.data(), fb); + CG_Base::f_cal(dim, move0.data(), move0.data(), xb); - setup_move(move, cg_gradn, best_x); - Ions_Move_Basic::move_atoms(ucell, move, pos); - trial = false; - xa = 0; - f_cal(move0, move, dim, xc); - xc = xb + xc; - xpt = xc; - Ions_Move_Basic::relax_bfgs_init = xc; - } - else + if ((std::abs(fb) < std::abs((fa) / 10.0))) { - double xtemp, ftemp; - f_cal(move0, grad, dim, fc); - fmin = std::abs(fc); - nbrent++; - - if ((fmin < std::abs((fmax) / 10.0)) || (nbrent > 3)) - { - nbrent = 0; - sd = true; - trial = true; - steplength = xpt; - flag = 1; - goto CG_begin; - } - else - { - Brent(fa, fb, fc, xa, xb, xc, best_x, xpt); // Brent method - if (xc < 0) - { - sd = true; - trial = true; - steplength = xb; - flag = 2; - goto CG_begin; - } - - normalize(cg_gradn, cg_grad0, dim); - setup_move(move, cg_gradn, best_x); - Ions_Move_Basic::move_atoms(ucell, move, pos); - Ions_Move_Basic::relax_bfgs_init = xc; - } + sd = true; + trial = true; + steplength = xb; + flag = 1; + continue; } - } - } - delete[] cg_grad; - delete[] grad; - delete[] cg_gradn; - delete[] pos; - delete[] move; + CG_Base::normalize(dim, cg_gradn.data(), cg_grad0.data()); + CG_Base::third_order(e0, e1, fa, fb, xb, best_x); - return; -} - -void Ions_Move_CG::setup_cg_grad(double *grad, - const double *grad0, - double *cg_grad, - const double *cg_grad0, - const int &ncggrad, - int &flag) -{ - ModuleBase::TITLE("Ions_Move_CG", "setup_cg_grad"); - assert(Ions_Move_Basic::istep > 0); - double gamma; - double cg0_cg, cg0_cg0, cg0_g; + if (best_x > 6 * xb || best_x < (-xb)) + { + best_x = 6 * xb; + } - if (ncggrad % 10000 == 0 || flag == 2) - { - for (int i = 0; i < dim; i++) - { - cg_grad[i] = grad[i]; - } - } - else - { - double gp_gp = 0.0; // grad_p.grad_p - double gg = 0.0; // grad.grad - double g_gp = 0.0; // grad_p.grad - double cgp_gp = 0.0; // cg_grad_p.grad_p - double cgp_g = 0.0; // cg_grad_p.grad - for (int i = 0; i < dim; i++) - { - gp_gp += grad0[i] * grad0[i]; - gg += grad[i] * grad[i]; - g_gp += grad0[i] * grad[i]; - cgp_gp += cg_grad0[i] * grad0[i]; - cgp_g += cg_grad0[i] * grad[i]; + CG_Base::setup_move(dim, move.data(), cg_gradn.data(), best_x); + Ions_Move_Basic::move_atoms(ucell, move.data(), pos.data()); + trial = false; + xa = 0; + CG_Base::f_cal(dim, move0.data(), move.data(), xc); + xc = xb + xc; + xpt = xc; + Ions_Move_Basic::relax_bfgs_init = xc; + break; } - assert(g_gp != 0.0); - const double gamma1 = gg / gp_gp; // FR - // const double gamma2 = -(gg - g_gp)/(cgp_g - cgp_gp); //CW - const double gamma2 = (gg - g_gp) / gp_gp; // PRP - // const double gamma = gg/cgp_gp; //D - // const double gamma = -gg/(cgp_g - cgp_gp); //D-Y - if (gamma1 < 0.5) - { - gamma = gamma1; - } - else - { - gamma = gamma2; - } + double xtemp = 0.0; + double ftemp = 0.0; + CG_Base::f_cal(dim, move0.data(), grad.data(), fc); + fmin = std::abs(fc); + nbrent++; - for (int i = 0; i < dim; i++) + if ((fmin < std::abs((fmax) / 10.0)) || (nbrent > 3)) { - // we can consider step as modified gradient. - cg_grad[i] = grad[i] + gamma * cg_grad0[i]; + nbrent = 0; + sd = true; + trial = true; + steplength = xpt; + flag = 1; + continue; } - } - return; -} -void Ions_Move_CG::third_order(const double &e0, - const double &e1, - const double &fa, - const double &fb, - const double x, - double &best_x) -{ - double k3, k2, k1; - double dmoveh, dmove1, dmove2, dmove, ecal1, ecal2; - - k3 = 3 * ((fb + fa) * x - 2 * (e1 - e0)) / (x * x * x); - k2 = (fb - fa) / x - k3 * x; - k1 = fa; - - dmoveh = x * fb / (fa - fb); - dmove1 = -k2 * (1 - sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); - dmove2 = -k2 * (1 + sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); - - if ((std::abs(k3 / k1) < 0.01) || ((k1 * k3 / (k2 * k2)) >= 0.25)) // this condition may be wrong - { - dmove = dmoveh; - } - else - { - dmove1 = -k2 * (1 - sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); - dmove2 = -k2 * (1 + sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); - ecal1 = k3 * dmove1 * dmove1 * dmove1 / 3 + k2 * dmove1 * dmove1 / 2 + k1 * dmove1; - ecal2 = k3 * dmove2 * dmove2 * dmove2 / 3 + k2 * dmove2 * dmove2 / 2 + k1 * dmove2; - if (ecal2 > ecal1) - dmove = dmove1 - x; - else - dmove = dmove2 - x; - - if (k3 < 0) - dmove = dmoveh; - } - - best_x = dmove; - return; -} - -void Ions_Move_CG::Brent(double &fa, - double &fb, - double &fc, - double &xa, - double &xb, - double &xc, - double &best_x, - double &xpt) -{ - double dmove; - double tmp; - double k2, k1, k0; - double xnew1, xnew2; - double ecalnew1, ecalnew2; - - if ((fa * fb) > 0) - { - dmove = (xc * fa - xa * fc) / (fa - fc); - if (dmove > 4 * xc) - // if(dmove > 4 * xc || dmove < 0) - { - dmove = 4 * xc; - } - xb = xc; - fb = fc; - } - else - { - k2 = -((fb - fc) / (xb - xc) - (fa - fc) / (xa - xc)) / (xa - xb); - k1 = (fa - fc) / (xa - xc) - k2 * (xa + xc); - k0 = fa - k1 * xa - k2 * xa * xa; - xnew1 = (-k1 - sqrt(k1 * k1 - 4 * k2 * k0)) / (2 * k2); - xnew2 = (-k1 + sqrt(k1 * k1 - 4 * k2 * k0)) / (2 * k2); - - if (xnew1 > xnew2) - { - tmp = xnew2; - xnew2 = xnew1; - xnew1 = tmp; - } - - ecalnew1 = k2 * xnew1 * xnew1 * xnew1 / 3 + k1 * xnew1 * xnew1 / 2 + k0 * xnew1; - ecalnew2 = k2 * xnew2 * xnew2 * xnew2 / 3 + k1 * xnew2 * xnew2 / 2 + k0 * xnew2; - dmove = xnew1; - - if (ecalnew1 > ecalnew2) - { - dmove = xnew2; - } - if (dmove < 0) - { - dmove = 2 * xc; // pengfei 14-6-5 - } - if (fa * fc > 0) + CG_Base::Brent(fa, fb, fc, xa, xb, xc, best_x, xpt); + if (xc < 0) { - xa = xc; - fa = fc; + sd = true; + trial = true; + steplength = xb; + flag = 2; + continue; } - if (fb * fc > 0) - { - xb = xc; - fb = fc; - } - } - - best_x = dmove - xpt; - xpt = dmove; - xc = dmove; - return; -} - -void Ions_Move_CG::f_cal(const double *g0, const double *g1, const int &dim, double &f_value) -{ - double hv0, hel; - hel = 0; - hv0 = 0; - for (int i = 0; i < dim; i++) - { - hel += g0[i] * g1[i]; - } - for (int i = 0; i < dim; i++) - { - hv0 += g0[i] * g0[i]; + CG_Base::normalize(dim, cg_gradn.data(), cg_grad0.data()); + CG_Base::setup_move(dim, move.data(), cg_gradn.data(), best_x); + Ions_Move_Basic::move_atoms(ucell, move.data(), pos.data()); + Ions_Move_Basic::relax_bfgs_init = xc; + break; } - f_value = hel / sqrt(hv0); return; -} - -void Ions_Move_CG::setup_move(double *move, double *cg_gradn, const double &trust_radius) -{ - // movement using gradient and trust_radius. - for (int i = 0; i < dim; ++i) - { - move[i] = -cg_gradn[i] * trust_radius; - } - return; -} - -void Ions_Move_CG::normalize(double *cg_gradn, const double *cg_grad, int dim) -{ - double norm = 0.0; - for (int i = 0; i < dim; ++i) - { - norm += pow(cg_grad[i], 2); - } - norm = sqrt(norm); - - if (norm != 0.0) - { - for (int i = 0; i < dim; ++i) - { - cg_gradn[i] = cg_grad[i] / norm; - } - } - return; -} +} \ No newline at end of file diff --git a/source/source_relax/ions_move_cg.h b/source/source_relax/ions_move_cg.h index 6e77127ca52..d846219e8ec 100644 --- a/source/source_relax/ions_move_cg.h +++ b/source/source_relax/ions_move_cg.h @@ -3,42 +3,28 @@ #include "source_base/matrix.h" #include "source_cell/unitcell.h" -class Ions_Move_CG +#include "cg_base.h" +#include + +class Ions_Move_CG : public CG_Base { public: Ions_Move_CG(); - ~Ions_Move_CG(); + ~Ions_Move_CG() = default; void allocate(void); void start(UnitCell &ucell, const ModuleBase::matrix &force, const double &etot); static double RELAX_CG_THR; - int sd_step=0; - int cg_step=0; + int sd_step = 0; + int cg_step = 0; private: - double * pos0 = nullptr; - double * grad0 = nullptr; - double * cg_grad0 = nullptr; - double * move0 = nullptr; - double e0=0.0; - // setup gradients. - void setup_cg_grad(double *grad, - const double *grad0, - double *cg_grad, - const double *cg_grad0, - const int &ncggrad, - int &flag); // LiuXh fix bug of lpf, 20180515 - void setup_move(double *move, double *cg_gradn, const double &trust_radius); - void Brent(double &fa, double &fb, double &fc, double &xa, double &xb, double &xc, double &best_x, double &xpt); - void f_cal(const double *g0, const double *g1, const int &dim, double &f_value); - void third_order(const double &e0, - const double &e1, - const double &fa, - const double &fb, - const double x, - double &best_x); - void normalize(double *cg_gradn, const double *cg_grad, int dim); + std::vector pos0; + std::vector grad0; + std::vector cg_grad0; + std::vector move0; + double e0 = 0.0; }; -#endif +#endif \ No newline at end of file diff --git a/source/source_relax/lattice_change_cg.cpp b/source/source_relax/lattice_change_cg.cpp index d0172ba7845..601da813181 100644 --- a/source/source_relax/lattice_change_cg.cpp +++ b/source/source_relax/lattice_change_cg.cpp @@ -3,25 +3,10 @@ #include "lattice_change_basic.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" +#include -// the 'dim' variable is defined in Lattice_Change_Basic using namespace Lattice_Change_Basic; -//=================== NOTES ======================== -// the approximate minimum of the total energy -// is calculated from a cubic (or quadratic) -// interpolation taking into account the change -// of the total energy and the changes of forces. -// ( 3 pieces of information ) -// 1. trial step -// 2. corrector step using cubic or quadratic interpolation -// 3. check the force and see if we need brent method interpolation. -// 4. a new trial step.... -// Brent's method is a complicated but popular -// root-finding algorithm combining the bisection method, -// the secant method and inverse quadratic interpolation -//=================== NOTES ======================== - Lattice_Change_CG::Lattice_Change_CG() { this->e0 = 0.0; @@ -30,7 +15,6 @@ Lattice_Change_CG::Lattice_Change_CG() void Lattice_Change_CG::allocate(void) { ModuleBase::TITLE("Lattice_Change_CG", "allocate"); - // mohan add 2021-02-07 assert(dim > 0); this->lat0.resize(dim, 0.0); @@ -48,423 +32,160 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_ assert(grad0.size() == static_cast(dim)); assert(cg_grad0.size() == static_cast(dim)); assert(move0.size() == static_cast(dim)); - - // sd , trial are two parameters, when sd=trial=true, - // a new direction begins, when sd = false trial =true static bool sd = false; - - // a cubic interpolation is used to make the third point, - // when sa = trial = false, we use Brent to get the - // minimum point in this direction. - static bool trial = false; - - // ncggrad is a parameter to control the cg method, - // every ten cg direction, we change the direction back to - // the steepest descent method + static bool trial = false; static int ncggrad = 0; - static double fa = 0.0; static double fb = 0.0; static double fc = 0.0; - static double xa = 0.0; static double xb = 0.0; static double xc = 0.0; - static double xpt = 0.0; static double steplength = 0.0; static double fmax = 0.0; - static int nbrent = 0; - double *lat = new double[dim]; - double *grad = new double[dim]; - double *cg_gradn = new double[dim]; - double *move = new double[dim]; - double *cg_grad = new double[dim]; - + std::vector lat(dim, 0.0); + std::vector grad(dim, 0.0); + std::vector cg_gradn(dim, 0.0); + std::vector move(dim, 0.0); + std::vector cg_grad(dim, 0.0); double best_x = 0.0; double fmin = 0.0; - int flag = 0; - ModuleBase::GlobalFunc::ZEROS(lat, dim); - ModuleBase::GlobalFunc::ZEROS(grad, dim); - ModuleBase::GlobalFunc::ZEROS(cg_gradn, dim); - ModuleBase::GlobalFunc::ZEROS(move, dim); - ModuleBase::GlobalFunc::ZEROS(cg_grad, dim); - -CG_begin: - - if (Lattice_Change_Basic::stress_step == 1) + while (true) { - steplength = Lattice_Change_Basic::lattice_change_ini; // read in the init trust radius - // cout<<"Lattice_Change_Basic::lattice_change_ini = "< 6 * xb || best_x < (-xb)) - { - best_x = 6 * xb; - } - - setup_move(move, cg_gradn, best_x); - Lattice_Change_Basic::change_lattice(ucell, move, lat); - - trial = false; - xa = 0; - f_cal(move0.data(), move, dim, xc); - xc = xb + xc; - xpt = xc; - - Lattice_Change_Basic::lattice_change_ini = xc; + sd = true; + trial = true; + steplength = xb; + flag = 1; + continue; } - else + + CG_Base::normalize(dim, cg_gradn.data(), cg_grad0.data()); + CG_Base::third_order(e0, e1, fa, fb, xb, best_x); + + if (best_x > 6 * xb || best_x < (-xb)) { - double xtemp, ftemp; - f_cal(move0.data(), grad, dim, fc); - - fmin = std::abs(fc); - nbrent++; - // cout<<"nbrent = "< 0); - double gamma = 0.0; - double cg0_cg, cg0_cg0, cg0_g; - - if (ncggrad % 10000 == 0 || flag == 2) - { - for (int i = 0; i < dim; i++) - { - cg_grad[i] = grad[i]; - } - } - else - { - double gp_gp = 0.0; // grad_p.grad_p - double gg = 0.0; // grad.grad - double g_gp = 0.0; // grad_p.grad - double cgp_gp = 0.0; // cg_grad_p.grad_p - double cgp_g = 0.0; // cg_grad_p.grad - for (int i = 0; i < dim; i++) - { - gp_gp += grad0[i] * grad0[i]; - gg += grad[i] * grad[i]; - g_gp += grad0[i] * grad[i]; - cgp_gp += cg_grad0[i] * grad0[i]; - cgp_g += cg_grad0[i] * grad[i]; + Lattice_Change_Basic::lattice_change_ini = xc; + break; } - assert(g_gp != 0.0); - const double gamma1 = gg / gp_gp; // FR - // const double gamma2 = -(gg - g_gp)/(cgp_g - cgp_gp); //CW - const double gamma2 = (gg - g_gp) / gp_gp; // PRP - // const double gamma = gg/cgp_gp; //D - // const double gamma = -gg/(cgp_g - cgp_gp); //D-Y + double xtemp = 0.0; + double ftemp = 0.0; + CG_Base::f_cal(dim, move0.data(), grad.data(), fc); - if (gamma1 < 0.5) - { - gamma = gamma1; - } - else - { - gamma = gamma2; - } + fmin = std::abs(fc); + nbrent++; - for (int i = 0; i < dim; i++) + if ((fmin < std::abs((fmax) / 10.0)) || (nbrent > 3)) { - // we can consider step as modified gradient. - cg_grad[i] = grad[i] + gamma * cg_grad0[i]; + nbrent = 0; + sd = true; + trial = true; + steplength = xpt; + flag = 1; + continue; } - } - return; -} - -void Lattice_Change_CG::third_order(const double &e0, - const double &e1, - const double &fa, - const double &fb, - const double x, - double &best_x) -{ - double k3, k2, k1; - double dmoveh, dmove1, dmove2, dmove, ecal1, ecal2; - k3 = 3 * ((fb + fa) * x - 2 * (e1 - e0)) / (x * x * x); - k2 = (fb - fa) / x - k3 * x; - k1 = fa; - - dmoveh = x * fb / (fa - fb); - dmove1 = -k2 * (1 - sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); - dmove2 = -k2 * (1 + sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); - - if ((std::abs(k3 / k1) < 0.01) || ((k1 * k3 / (k2 * k2)) >= 0.25)) // this condition may be wrong - { - dmove = dmoveh; - } - else - { - dmove1 = -k2 * (1 - sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); - dmove2 = -k2 * (1 + sqrt(1 - 4 * k1 * k3 / (k2 * k2))) / (2 * k3); - ecal1 = k3 * dmove1 * dmove1 * dmove1 / 3 + k2 * dmove1 * dmove1 / 2 + k1 * dmove1; - ecal2 = k3 * dmove2 * dmove2 * dmove2 / 3 + k2 * dmove2 * dmove2 / 2 + k1 * dmove2; - if (ecal2 > ecal1) - dmove = dmove1 - x; - else - dmove = dmove2 - x; - - if (k3 < 0) - dmove = dmoveh; - } - - best_x = dmove; - return; -} - -void Lattice_Change_CG::Brent(double &fa, - double &fb, - double &fc, - double &xa, - double &xb, - double &xc, - double &best_x, - double &xpt) -{ - double dmove; - double tmp; - double k2, k1, k0; - double xnew1, xnew2; - double ecalnew1, ecalnew2; - - if ((fa * fb) > 0) - { - dmove = (xc * fa - xa * fc) / (fa - fc); - if (dmove > 4 * xc) - // if(dmove > 4 * xc || dmove < 0) - { - dmove = 4 * xc; - } - xb = xc; - fb = fc; - } - else - { - k2 = -((fb - fc) / (xb - xc) - (fa - fc) / (xa - xc)) / (xa - xb); - k1 = (fa - fc) / (xa - xc) - k2 * (xa + xc); - k0 = fa - k1 * xa - k2 * xa * xa; - xnew1 = (-k1 - sqrt(k1 * k1 - 4 * k2 * k0)) / (2 * k2); - xnew2 = (-k1 + sqrt(k1 * k1 - 4 * k2 * k0)) / (2 * k2); - - if (xnew1 > xnew2) + CG_Base::Brent(fa, fb, fc, xa, xb, xc, best_x, xpt); + if (xc < 0) { - tmp = xnew2; - xnew2 = xnew1; - xnew1 = tmp; + sd = true; + trial = true; + steplength = xb; + flag = 2; + continue; } - ecalnew1 = k2 * xnew1 * xnew1 * xnew1 / 3 + k1 * xnew1 * xnew1 / 2 + k0 * xnew1; - ecalnew2 = k2 * xnew2 * xnew2 * xnew2 / 3 + k1 * xnew2 * xnew2 / 2 + k0 * xnew2; - dmove = xnew1; + CG_Base::normalize(dim, cg_gradn.data(), cg_grad0.data()); + CG_Base::setup_move(dim, move.data(), cg_gradn.data(), best_x); + Lattice_Change_Basic::change_lattice(ucell, move.data(), lat.data()); - if (ecalnew1 > ecalnew2) - { - dmove = xnew2; - } - if (dmove < 0) - { - dmove = 2 * xc; // pengfei 14-6-5 - } - if (fa * fc > 0) - { - xa = xc; - fa = fc; - } - if (fb * fc > 0) - { - xb = xc; - fb = fc; - } + Lattice_Change_Basic::lattice_change_ini = xc; + break; } - best_x = dmove - xpt; - xpt = dmove; - xc = dmove; - return; -} - -void Lattice_Change_CG::f_cal(const double *g0, const double *g1, const int &dim, double &f_value) -{ - double hv0, hel; - hel = 0; - hv0 = 0; - for (int i = 0; i < dim; i++) - { - hel += g0[i] * g1[i]; - } - for (int i = 0; i < dim; i++) - { - hv0 += g0[i] * g0[i]; - } - - f_value = hel / sqrt(hv0); - return; -} - -void Lattice_Change_CG::setup_move(double *move, double *cg_gradn, const double &trust_radius) -{ - // movement using gradient and trust_radius. - for (int i = 0; i < dim; ++i) - { - move[i] = -cg_gradn[i] * trust_radius; - } - return; -} - -void Lattice_Change_CG::normalize(double *cg_gradn, const double *cg_grad, int dim) -{ - double norm = 0.0; - for (int i = 0; i < dim; ++i) - { - norm += pow(cg_grad[i], 2); - } - norm = sqrt(norm); - - if (norm != 0.0) - { - for (int i = 0; i < dim; ++i) - { - cg_gradn[i] = cg_grad[i] / norm; - } - } - return; -} +} \ No newline at end of file diff --git a/source/source_relax/lattice_change_cg.h b/source/source_relax/lattice_change_cg.h index 76d8dc70f8d..a2302cb26d9 100644 --- a/source/source_relax/lattice_change_cg.h +++ b/source/source_relax/lattice_change_cg.h @@ -3,8 +3,10 @@ #include "source_base/matrix.h" #include "source_cell/unitcell.h" +#include "cg_base.h" #include -class Lattice_Change_CG + +class Lattice_Change_CG : public CG_Base { public: @@ -19,30 +21,7 @@ class Lattice_Change_CG std::vector grad0; std::vector cg_grad0; std::vector move0; - double e0=0.0; - - // setup gradients. - void setup_cg_grad(double *grad, - const double *grad0, - double *cg_grad, - const double *cg_grad0, - const int &ncggrad, - int &flag); - - void setup_move(double *move, double *cg_gradn, const double &trust_radius); - - void Brent(double &fa, double &fb, double &fc, double &xa, double &xb, double &xc, double &best_x, double &xpt); - - void f_cal(const double *g0, const double *g1, const int &dim, double &f_value); - - void third_order(const double &e0, - const double &e1, - const double &fa, - const double &fb, - const double x, - double &best_x); - - void normalize(double *cg_gradn, const double *cg_grad, int dim); + double e0 = 0.0; }; -#endif +#endif \ No newline at end of file From e5fc0a9c5e2a2dc400f87b77d3e6f601c63beb86 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 9 Jun 2026 20:01:59 +0800 Subject: [PATCH 13/16] fix --- source/source_relax/test/lattice_change_methods_test.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/source/source_relax/test/lattice_change_methods_test.cpp b/source/source_relax/test/lattice_change_methods_test.cpp index 64173e110e1..dfc7fbddb5f 100644 --- a/source/source_relax/test/lattice_change_methods_test.cpp +++ b/source/source_relax/test/lattice_change_methods_test.cpp @@ -20,10 +20,6 @@ Lattice_Change_CG::Lattice_Change_CG() { } -Lattice_Change_CG::~Lattice_Change_CG() -{ -} - void Lattice_Change_CG::allocate(void) { } From 95c3ca7a6791bc2ad47f9b9e740d08cbec476f41 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 9 Jun 2026 20:55:22 +0800 Subject: [PATCH 14/16] fix test of relax module --- source/source_relax/cg_base.h | 40 +++++++------- source/source_relax/test/CMakeLists.txt | 21 +++++--- source/source_relax/test/bfgs_basic_test.cpp | 54 +++---------------- source/source_relax/test/bfgs_test.cpp | 1 + .../test/ions_move_basic_test.cpp | 12 +---- .../source_relax/test/ions_move_bfgs_test.cpp | 9 +--- .../source_relax/test/ions_move_cg_test.cpp | 36 +++++-------- .../source_relax/test/ions_move_sd_test.cpp | 14 ++--- .../test/lattice_change_basic_test.cpp | 10 +--- .../test/lattice_change_cg_test.cpp | 36 +++++-------- 10 files changed, 77 insertions(+), 156 deletions(-) diff --git a/source/source_relax/cg_base.h b/source/source_relax/cg_base.h index ce245830ce4..e2a3d929121 100644 --- a/source/source_relax/cg_base.h +++ b/source/source_relax/cg_base.h @@ -10,7 +10,7 @@ * The functions are designed to be independent of specific data sources. */ class CG_Base { -protected: +public: /** * @brief Setup conjugate gradient direction. * @param dim Dimension of the optimization problem @@ -25,15 +25,6 @@ class CG_Base { double *cg_grad, const double *cg_grad0, const int &ncggrad, int &flag); - /** - * @brief Calculate movement vector from CG direction. - * @param dim Dimension of the optimization problem - * @param move Output movement vector - * @param cg_gradn Normalized CG direction - * @param trust_radius Step size limit - */ - void setup_move(const int dim, double *move, double *cg_gradn, const double &trust_radius); - /** * @brief Brent's method for one-dimensional minimization. * @param fa Function value at xa @@ -49,15 +40,6 @@ class CG_Base { double &xa, double &xb, double &xc, double &best_x, double &xpt); - /** - * @brief Calculate projection of gradient onto search direction. - * @param dim Dimension of vectors - * @param g0 First vector (usually search direction) - * @param g1 Second vector (usually gradient) - * @param f_value Output projection value - */ - void f_cal(const int dim, const double *g0, const double *g1, double &f_value); - /** * @brief Third-order polynomial interpolation for line search. * @param e0 Energy at previous step @@ -71,6 +53,26 @@ class CG_Base { const double &fa, const double &fb, const double x, double &best_x); + /** + * @brief Calculate projection of gradient onto search direction. + * @param dim Dimension of vectors + * @param g0 First vector (usually search direction) + * @param g1 Second vector (usually gradient) + * @param f_value Output projection value + */ + void f_cal(const int dim, const double *g0, const double *g1, double &f_value); + + /** + * @brief Calculate movement vector from CG direction. + * @param dim Dimension of the optimization problem + * @param move Output movement vector + * @param cg_gradn Normalized CG direction + * @param trust_radius Step size limit + */ + void setup_move(const int dim, double *move, double *cg_gradn, const double &trust_radius); + +protected: + /** * @brief Normalize vector. * @param dim Dimension of vector diff --git a/source/source_relax/test/CMakeLists.txt b/source/source_relax/test/CMakeLists.txt index 2c19c41ff48..95b38fe3ba8 100644 --- a/source/source_relax/test/CMakeLists.txt +++ b/source/source_relax/test/CMakeLists.txt @@ -32,13 +32,13 @@ list(APPEND cell_source_files AddTest( TARGET MODULE_RELAX_lattice_change_methods_test LIBS parameter ${math_libs} base device - SOURCES lattice_change_methods_test.cpp ../lattice_change_methods.cpp ../lattice_change_basic.cpp mock_remake_cell.cpp + SOURCES lattice_change_methods_test.cpp ../lattice_change_methods.cpp ../lattice_change_basic.cpp ../relax_data.cpp mock_remake_cell.cpp ) AddTest( TARGET MODULE_RELAX_lattice_change_basic_test LIBS parameter ${math_libs} base device - SOURCES lattice_change_basic_test.cpp ../lattice_change_basic.cpp mock_remake_cell.cpp + SOURCES lattice_change_basic_test.cpp ../lattice_change_basic.cpp ../relax_data.cpp mock_remake_cell.cpp ) AddTest( @@ -47,27 +47,31 @@ AddTest( SOURCES lattice_change_cg_test.cpp ../lattice_change_cg.cpp ../lattice_change_basic.cpp + ../cg_base.cpp + ../relax_data.cpp mock_remake_cell.cpp ../../source_io/module_output/orb_io.cpp ) AddTest( TARGET MODULE_RELAX_bfgs_basic_test - LIBS parameter ${math_libs} base device - SOURCES bfgs_basic_test.cpp ../bfgs_basic.cpp + LIBS parameter ${math_libs} base device symmetry + SOURCES bfgs_basic_test.cpp ../bfgs_basic.cpp ../relax_data.cpp ../ions_move_basic.cpp + ../../source_io/module_output/orb_io.cpp + ${cell_source_files} ) AddTest( TARGET MODULE_RELAX_bfgs_test LIBS parameter ${math_libs} base device - SOURCES bfgs_test.cpp ../ions_move_bfgs2.cpp ../ions_move_basic.cpp ../matrix_methods.cpp ${cell_source_files} + SOURCES bfgs_test.cpp ../ions_move_bfgs2.cpp ../ions_move_basic.cpp ../matrix_methods.cpp ../relax_data.cpp ${cell_source_files} ) AddTest( TARGET MODULE_RELAX_ions_move_basic_test LIBS parameter ${math_libs} base device - SOURCES ions_move_basic_test.cpp ../ions_move_basic.cpp ${cell_source_files} + SOURCES ions_move_basic_test.cpp ../ions_move_basic.cpp ../relax_data.cpp ${cell_source_files} ) AddTest( @@ -77,6 +81,7 @@ AddTest( ../ions_move_bfgs.cpp ../ions_move_basic.cpp ../bfgs_basic.cpp + ../relax_data.cpp ../../source_io/module_output/orb_io.cpp ${cell_source_files} ) @@ -86,7 +91,9 @@ AddTest( LIBS parameter ${math_libs} base device SOURCES ions_move_cg_test.cpp ../ions_move_cg.cpp + ../cg_base.cpp ../ions_move_basic.cpp + ../relax_data.cpp ../../source_io/module_output/orb_io.cpp ${cell_source_files} ) @@ -94,5 +101,5 @@ AddTest( AddTest( TARGET MODULE_RELAX_ions_move_sd_test LIBS parameter ${math_libs} base device - SOURCES ions_move_sd_test.cpp ../ions_move_sd.cpp ../ions_move_basic.cpp ${cell_source_files} + SOURCES ions_move_sd_test.cpp ../ions_move_sd.cpp ../ions_move_basic.cpp ../relax_data.cpp ${cell_source_files} ) \ No newline at end of file diff --git a/source/source_relax/test/bfgs_basic_test.cpp b/source/source_relax/test/bfgs_basic_test.cpp index 91c028ec310..efd30e9829e 100644 --- a/source/source_relax/test/bfgs_basic_test.cpp +++ b/source/source_relax/test/bfgs_basic_test.cpp @@ -1,4 +1,5 @@ #include "source_relax/ions_move_basic.h" +#include "source_relax/relax_data.h" #include "gmock/gmock.h" #define private public #include "source_io/module_parameter/parameter.h" @@ -13,44 +14,6 @@ * unit tests of class BFGS_Basic ***********************************************/ -/** - * - Tested Functions: - * - BFGS_Basic::allocate_basic() - * - BFGS_Basic::new_step() - * - BFGS_Basic::reset_hessian() - * - BFGS_Basic::save_bfgs() - * - BFGS_Basic::check_move() - * - BFGS_Basic::update_inverse_hessian() - * - BFGS_Basic::check_wolfe_conditions() - * - BFGS_Basic::compute_trust_radius() - */ - -int Ions_Move_Basic::dim = 0; -bool Ions_Move_Basic::converged = false; -double Ions_Move_Basic::largest_grad = 0.0; -int Ions_Move_Basic::update_iter = 0; -int Ions_Move_Basic::istep = 0; -double Ions_Move_Basic::ediff = 0.0; -double Ions_Move_Basic::etot = 0.0; -double Ions_Move_Basic::etot_p = 0.0; -double Ions_Move_Basic::trust_radius = 0.0; -double Ions_Move_Basic::trust_radius_old = 0.0; -double Ions_Move_Basic::relax_bfgs_rmax = -1.0; -double Ions_Move_Basic::relax_bfgs_rmin = -1.0; -double Ions_Move_Basic::relax_bfgs_init = -1.0; -double Ions_Move_Basic::best_xxx = 1.0; -int Ions_Move_Basic::out_stru = 0; - -double Ions_Move_Basic::dot_func(const double *a, const double *b, const int &dim_in) -{ - double result = 0.0; - for (int i = 0; i < dim_in; i++) - { - result += a[i] * b[i]; - } - return result; -} - class BFGSBasicTest : public ::testing::Test { protected: @@ -73,14 +36,13 @@ TEST_F(BFGSBasicTest, TestAllocate) Ions_Move_Basic::dim = 4; bfgs.allocate_basic(); - // Check if allocated arrays are not empty - EXPECT_NE(nullptr, bfgs.pos); - EXPECT_NE(nullptr, bfgs.pos_p); - EXPECT_NE(nullptr, bfgs.grad); - EXPECT_NE(nullptr, bfgs.grad_p); - EXPECT_NE(nullptr, bfgs.move); - EXPECT_NE(nullptr, bfgs.move_p); - EXPECT_NE(nullptr, bfgs.inv_hess.c); + // Check if allocated vectors are not empty + EXPECT_EQ(bfgs.pos.size(), 4U); + EXPECT_EQ(bfgs.pos_p.size(), 4U); + EXPECT_EQ(bfgs.grad.size(), 4U); + EXPECT_EQ(bfgs.grad_p.size(), 4U); + EXPECT_EQ(bfgs.move.size(), 4U); + EXPECT_EQ(bfgs.move_p.size(), 4U); } // Test if a dimension less than or equal to 0 results in an assertion error diff --git a/source/source_relax/test/bfgs_test.cpp b/source/source_relax/test/bfgs_test.cpp index 81852def15f..081012dba39 100644 --- a/source/source_relax/test/bfgs_test.cpp +++ b/source/source_relax/test/bfgs_test.cpp @@ -12,6 +12,7 @@ #undef private #include "source_relax/ions_move_basic.h" // for Ions_Move_Basic static members +#include "source_relax/relax_data.h" /************************************************ * unit tests for Ions_Move_BFGS2 (no MockUnitCell) diff --git a/source/source_relax/test/ions_move_basic_test.cpp b/source/source_relax/test/ions_move_basic_test.cpp index d9e9caaf338..b6f02628aec 100644 --- a/source/source_relax/test/ions_move_basic_test.cpp +++ b/source/source_relax/test/ions_move_basic_test.cpp @@ -4,23 +4,13 @@ #include "source_io/module_parameter/parameter.h" #undef private #include "source_relax/ions_move_basic.h" +#include "source_relax/relax_data.h" #include "for_test.h" /************************************************ * unit tests of namespace Ions_Move_Basic ***********************************************/ -/** - * - Tested Functions: - * - Ions_Move_Basic::setup_gradient() - * - Ions_Move_Basic::move_atoms() - * - Ions_Move_Basic::check_converged() - * - Ions_Move_Basic::terminate() - * - Ions_Move_Basic::setup_etot() - * - Ions_Move_Basic::dot_func() - * - Ions_Move_Basic::third_order() - */ - // Define a fixture for the tests class IonsMoveBasicTest : public ::testing::Test { diff --git a/source/source_relax/test/ions_move_bfgs_test.cpp b/source/source_relax/test/ions_move_bfgs_test.cpp index bbf9bfadf93..f13ef8c3ea5 100644 --- a/source/source_relax/test/ions_move_bfgs_test.cpp +++ b/source/source_relax/test/ions_move_bfgs_test.cpp @@ -8,18 +8,11 @@ #include "source_relax/ions_move_bfgs.h" #undef private #undef protected + /************************************************ * unit tests of class Ions_Move_BFGS ***********************************************/ -/** - * - Tested Functions: - * - Ions_Move_BFGS::allocate() - * - Ions_Move_BFGS::start() - * - Ions_Move_BFGS::bfgs_routine() - * - Ions_Move_BFGS::restart_bfgs() - */ - // Define a fixture for the tests class IonsMoveBFGSTest : public ::testing::Test { diff --git a/source/source_relax/test/ions_move_cg_test.cpp b/source/source_relax/test/ions_move_cg_test.cpp index c0aa34dceca..83c4323be11 100644 --- a/source/source_relax/test/ions_move_cg_test.cpp +++ b/source/source_relax/test/ions_move_cg_test.cpp @@ -7,21 +7,11 @@ #include "source_relax/ions_move_basic.h" #include "source_relax/ions_move_cg.h" #undef private + /************************************************ * unit tests of class Ions_Move_CG ***********************************************/ -/** - * - Tested Functions: - * - Ions_Move_CG::allocate() - * - Ions_Move_CG::start() - * - Ions_Move_CG::setup_cg_grad() - * - Ions_Move_CG::setup_move() - * - Ions_Move_CG::Brent() - * - Ions_Move_CG::f_cal() - * - Ions_Move_CG::third_order() - */ - class IonsMoveCGTest : public ::testing::Test { protected: @@ -70,11 +60,11 @@ TEST_F(IonsMoveCGTest, TestAllocate) Ions_Move_Basic::dim = 4; im_cg.allocate(); - // Check if allocated arrays are not empty - EXPECT_NE(nullptr, im_cg.pos0); - EXPECT_NE(nullptr, im_cg.grad0); - EXPECT_NE(nullptr, im_cg.cg_grad0); - EXPECT_NE(nullptr, im_cg.move0); + // Check if allocated vectors are not empty + EXPECT_EQ(im_cg.pos0.size(), 4U); + EXPECT_EQ(im_cg.grad0.size(), 4U); + EXPECT_EQ(im_cg.cg_grad0.size(), 4U); + EXPECT_EQ(im_cg.move0.size(), 4U); } // Test if a dimension less than or equal to 0 results in an assertion error @@ -374,7 +364,7 @@ TEST_F(IonsMoveCGTest, SetupCgGradNcggradIsMultipleOf10000) int ncggrad = 50000; // multiple of 10000 int flag = 0; - im_cg.setup_cg_grad(grad, grad0, cggrad, cggrad0, ncggrad, flag); + im_cg.setup_cg_grad(Ions_Move_Basic::dim, grad, grad0, cggrad, cggrad0, ncggrad, flag); EXPECT_DOUBLE_EQ(cggrad[0], grad[0]); EXPECT_DOUBLE_EQ(cggrad[1], grad[1]); @@ -394,7 +384,7 @@ TEST_F(IonsMoveCGTest, SetupCgGradNcggradIsNotMultipleOf10000Case1) int ncggrad = 100; int flag = 0; - im_cg.setup_cg_grad(grad, grad0, cggrad, cggrad0, ncggrad, flag); + im_cg.setup_cg_grad(Ions_Move_Basic::dim, grad, grad0, cggrad, cggrad0, ncggrad, flag); EXPECT_DOUBLE_EQ(cggrad[0], 1.25); EXPECT_DOUBLE_EQ(cggrad[1], 0.0); @@ -414,7 +404,7 @@ TEST_F(IonsMoveCGTest, SetupCgGradNcggradIsNotMultipleOf10000Case2) int ncggrad = 100; int flag = 0; - im_cg.setup_cg_grad(grad, grad0, cggrad, cggrad0, ncggrad, flag); + im_cg.setup_cg_grad(Ions_Move_Basic::dim, grad, grad0, cggrad, cggrad0, ncggrad, flag); EXPECT_DOUBLE_EQ(cggrad[0], grad[0]); EXPECT_DOUBLE_EQ(cggrad[1], grad[1]); @@ -571,9 +561,9 @@ TEST_F(IonsMoveCGTest, Fcal) Ions_Move_Basic::dim = 9; double g0[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; double g1[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; - double f_value; + double f_value = 0.0; - im_cg.f_cal(g0, g1, Ions_Move_Basic::dim, f_value); + im_cg.f_cal(Ions_Move_Basic::dim, g0, g1, f_value); EXPECT_DOUBLE_EQ(f_value, 3.0); } @@ -584,9 +574,9 @@ TEST_F(IonsMoveCGTest, SetupMove) Ions_Move_Basic::dim = 9; double trust_radius = 1.0; double cg_gradn[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; - double move[9]; + double move[9] = {0.0}; - im_cg.setup_move(move, cg_gradn, trust_radius); + im_cg.setup_move(Ions_Move_Basic::dim, move, cg_gradn, trust_radius); EXPECT_DOUBLE_EQ(move[0], -1.0); EXPECT_DOUBLE_EQ(move[1], -1.0); diff --git a/source/source_relax/test/ions_move_sd_test.cpp b/source/source_relax/test/ions_move_sd_test.cpp index 26d5afe2332..647247ecd38 100644 --- a/source/source_relax/test/ions_move_sd_test.cpp +++ b/source/source_relax/test/ions_move_sd_test.cpp @@ -7,17 +7,11 @@ #include "source_relax/ions_move_basic.h" #include "source_relax/ions_move_sd.h" #undef private + /************************************************ * unit tests of class Ions_Move_SD ***********************************************/ -/** - * - Tested Functions: - * - Ions_Move_SD::allocate() - * - Ions_Move_SD::start() - * - Ions_Move_SD::cal_tradius_sd() - */ - class IonsMoveSDTest : public ::testing::Test { protected: @@ -44,9 +38,9 @@ TEST_F(IonsMoveSDTest, TestAllocate) Ions_Move_Basic::dim = 4; im_sd.allocate(); - // Check if allocated arrays are not empty - EXPECT_NE(nullptr, im_sd.grad_saved); - EXPECT_NE(nullptr, im_sd.pos_saved); + // Check if allocated vectors are not empty + EXPECT_EQ(im_sd.grad_saved.size(), 4U); + EXPECT_EQ(im_sd.pos_saved.size(), 4U); } // Test if a dimension less than or equal to 0 results in an assertion error diff --git a/source/source_relax/test/lattice_change_basic_test.cpp b/source/source_relax/test/lattice_change_basic_test.cpp index b40af6d88a9..2742d39f5b9 100644 --- a/source/source_relax/test/lattice_change_basic_test.cpp +++ b/source/source_relax/test/lattice_change_basic_test.cpp @@ -1,4 +1,5 @@ #include "source_relax/lattice_change_basic.h" +#include "source_relax/relax_data.h" #include "mock_remake_cell.h" #include "for_test.h" @@ -12,15 +13,6 @@ * unit tests of namespace Lattice_Change_Basic ***********************************************/ -/** - * - Tested Functions: - * - Lattice_Change_Basic::setup_gradient() - * - Lattice_Change_Basic::change_lattice() - * - Lattice_Change_Basic::check_converged() - * - Lattice_Change_Basic::terminate() - * - Lattice_Change_Basic::setup_etot() - */ - // Define a fixture for the tests class LatticeChangeBasicTest : public ::testing::Test { diff --git a/source/source_relax/test/lattice_change_cg_test.cpp b/source/source_relax/test/lattice_change_cg_test.cpp index 01c863b8eed..11a320d3712 100644 --- a/source/source_relax/test/lattice_change_cg_test.cpp +++ b/source/source_relax/test/lattice_change_cg_test.cpp @@ -5,21 +5,11 @@ #include "source_relax/lattice_change_basic.h" #include "source_relax/lattice_change_cg.h" #undef private + /************************************************ * unit tests of class Lattice_Change_CG ***********************************************/ -/** - * - Tested Functions: - * - Lattice_Change_CG::allocate() - * - Lattice_Change_CG::start() - * - Lattice_Change_CG::setup_cg_grad() - * - Lattice_Change_CG::setup_move() - * - Lattice_Change_CG::Brent() - * - Lattice_Change_CG::f_cal() - * - Lattice_Change_CG::third_order() - */ - class LatticeChangeCGTest : public ::testing::Test { protected: @@ -47,11 +37,11 @@ TEST_F(LatticeChangeCGTest, TestAllocate) Lattice_Change_Basic::dim = 4; lc_cg.allocate(); - // Check if allocated arrays are not empty - EXPECT_NE(nullptr, lc_cg.lat0); - EXPECT_NE(nullptr, lc_cg.grad0); - EXPECT_NE(nullptr, lc_cg.cg_grad0); - EXPECT_NE(nullptr, lc_cg.move0); + // Check if allocated vectors are not empty + EXPECT_EQ(lc_cg.lat0.size(), 4U); + EXPECT_EQ(lc_cg.grad0.size(), 4U); + EXPECT_EQ(lc_cg.cg_grad0.size(), 4U); + EXPECT_EQ(lc_cg.move0.size(), 4U); } // Test if a dimension less than or equal to 0 results in an assertion error @@ -300,7 +290,7 @@ TEST_F(LatticeChangeCGTest, SetupCgGradNcggradIsMultipleOf10000) int ncggrad = 50000; // multiple of 10000 int flag = 0; - lc_cg.setup_cg_grad(grad, grad0, cggrad, cggrad0, ncggrad, flag); + lc_cg.setup_cg_grad(Lattice_Change_Basic::dim, grad, grad0, cggrad, cggrad0, ncggrad, flag); EXPECT_DOUBLE_EQ(cggrad[0], grad[0]); EXPECT_DOUBLE_EQ(cggrad[1], grad[1]); @@ -323,7 +313,7 @@ TEST_F(LatticeChangeCGTest, SetupCgGradNcggradIsNotMultipleOf10000Case1) int ncggrad = 100; int flag = 0; - lc_cg.setup_cg_grad(grad, grad0, cggrad, cggrad0, ncggrad, flag); + lc_cg.setup_cg_grad(Lattice_Change_Basic::dim, grad, grad0, cggrad, cggrad0, ncggrad, flag); EXPECT_DOUBLE_EQ(cggrad[0], 1.25); EXPECT_DOUBLE_EQ(cggrad[1], 0.0); @@ -346,7 +336,7 @@ TEST_F(LatticeChangeCGTest, SetupCgGradNcggradIsNotMultipleOf10000Case2) int ncggrad = 100; int flag = 0; - lc_cg.setup_cg_grad(grad, grad0, cggrad, cggrad0, ncggrad, flag); + lc_cg.setup_cg_grad(Lattice_Change_Basic::dim, grad, grad0, cggrad, cggrad0, ncggrad, flag); EXPECT_DOUBLE_EQ(cggrad[0], grad[0]); EXPECT_DOUBLE_EQ(cggrad[1], grad[1]); @@ -506,9 +496,9 @@ TEST_F(LatticeChangeCGTest, Fcal) Lattice_Change_Basic::dim = 9; double g0[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; double g1[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; - double f_value; + double f_value = 0.0; - lc_cg.f_cal(g0, g1, Lattice_Change_Basic::dim, f_value); + lc_cg.f_cal(Lattice_Change_Basic::dim, g0, g1, f_value); EXPECT_DOUBLE_EQ(f_value, 3.0); } @@ -519,9 +509,9 @@ TEST_F(LatticeChangeCGTest, SetupMove) Lattice_Change_Basic::dim = 9; double trust_radius = 1.0; double cg_gradn[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; - double move[9]; + double move[9] = {0.0}; - lc_cg.setup_move(move, cg_gradn, trust_radius); + lc_cg.setup_move(Lattice_Change_Basic::dim, move, cg_gradn, trust_radius); EXPECT_DOUBLE_EQ(move[0], -1.0); EXPECT_DOUBLE_EQ(move[1], -1.0); From 18fbff9907953337efa4e8692dd27f50845baf92 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 9 Jun 2026 21:07:31 +0800 Subject: [PATCH 15/16] remove ZEROS function --- source/source_relax/bfgs_basic.cpp | 16 ++++++---------- source/source_relax/ions_move_basic.cpp | 5 +++-- source/source_relax/ions_move_bfgs.cpp | 7 ++++--- source/source_relax/ions_move_sd.cpp | 10 ++++------ 4 files changed, 17 insertions(+), 21 deletions(-) diff --git a/source/source_relax/bfgs_basic.cpp b/source/source_relax/bfgs_basic.cpp index f542d209a30..36fc12f0a15 100644 --- a/source/source_relax/bfgs_basic.cpp +++ b/source/source_relax/bfgs_basic.cpp @@ -1,4 +1,5 @@ #include "bfgs_basic.h" +#include #include "source_io/module_parameter/parameter.h" #include "ions_move_basic.h" #include "source_base/global_function.h" @@ -36,10 +37,8 @@ void BFGS_Basic::update_inverse_hessian(const double &lat0) // ModuleBase::TITLE("Ions_Move_BFGS","update_inverse_hessian"); assert(dim > 0); - std::vector s(dim); - std::vector y(dim); - ModuleBase::GlobalFunc::ZEROS(s.data(), dim); - ModuleBase::GlobalFunc::ZEROS(y.data(), dim); + std::vector s(dim, 0.0); + std::vector y(dim, 0.0); for (int i = 0; i < dim; i++) { @@ -64,12 +63,9 @@ void BFGS_Basic::update_inverse_hessian(const double &lat0) return; } - std::vector Hs(dim); - std::vector Hy(dim); - std::vector yH(dim); - ModuleBase::GlobalFunc::ZEROS(Hs.data(), dim); - ModuleBase::GlobalFunc::ZEROS(Hy.data(), dim); - ModuleBase::GlobalFunc::ZEROS(yH.data(), dim); + std::vector Hs(dim, 0.0); + std::vector Hy(dim, 0.0); + std::vector yH(dim, 0.0); for (int i = 0; i < dim; i++) { diff --git a/source/source_relax/ions_move_basic.cpp b/source/source_relax/ions_move_basic.cpp index c2b3d0cd9c8..08b45d29ce8 100644 --- a/source/source_relax/ions_move_basic.cpp +++ b/source/source_relax/ions_move_basic.cpp @@ -1,5 +1,6 @@ #include "ions_move_basic.h" +#include #include "source_io/module_parameter/parameter.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" @@ -26,8 +27,8 @@ void Ions_Move_Basic::setup_gradient(const UnitCell &ucell, const ModuleBase::ma assert(grad != nullptr); assert(dim == 3 * ucell.nat); - ModuleBase::GlobalFunc::ZEROS(pos, dim); - ModuleBase::GlobalFunc::ZEROS(grad, dim); + std::fill_n(pos, dim, 0.0); + std::fill_n(grad, dim, 0.0); // (1) init gradient // the unit of pos: Bohr. diff --git a/source/source_relax/ions_move_bfgs.cpp b/source/source_relax/ions_move_bfgs.cpp index 31c9f33f86d..11b712805ba 100644 --- a/source/source_relax/ions_move_bfgs.cpp +++ b/source/source_relax/ions_move_bfgs.cpp @@ -1,5 +1,6 @@ #include "ions_move_bfgs.h" +#include #include "source_io/module_parameter/parameter.h" #include "ions_move_basic.h" #include "source_base/global_function.h" @@ -130,9 +131,9 @@ void Ions_Move_BFGS::restart_bfgs(const double& lat0) else { // bfgs initialization - ModuleBase::GlobalFunc::ZEROS(pos_p.data(), dim); - ModuleBase::GlobalFunc::ZEROS(grad_p.data(), dim); - ModuleBase::GlobalFunc::ZEROS(move_p.data(), dim); + std::fill(pos_p.begin(), pos_p.end(), 0.0); + std::fill(grad_p.begin(), grad_p.end(), 0.0); + std::fill(move_p.begin(), move_p.end(), 0.0); Ions_Move_Basic::update_iter = 0; diff --git a/source/source_relax/ions_move_sd.cpp b/source/source_relax/ions_move_sd.cpp index 8302e509a46..2e0f4f7a322 100644 --- a/source/source_relax/ions_move_sd.cpp +++ b/source/source_relax/ions_move_sd.cpp @@ -1,5 +1,6 @@ #include "ions_move_sd.h" +#include #include "source_io/module_parameter/parameter.h" #include "ions_move_basic.h" #include "source_base/global_function.h" @@ -27,12 +28,9 @@ void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const assert(grad_saved.size() == static_cast(dim)); assert(pos_saved.size() == static_cast(dim)); - std::vector pos(dim); - std::vector grad(dim); - std::vector move(dim); - ModuleBase::GlobalFunc::ZEROS(pos.data(), dim); - ModuleBase::GlobalFunc::ZEROS(grad.data(), dim); - ModuleBase::GlobalFunc::ZEROS(move.data(), dim); + std::vector pos(dim, 0.0); + std::vector grad(dim, 0.0); + std::vector move(dim, 0.0); // 1: ediff = 0 // 0: ediff < 0 From fee96c227e0f0cd4a5be5b5ca5b1b79d28478f7f Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Tue, 9 Jun 2026 22:01:53 +0800 Subject: [PATCH 16/16] remove GlobalC, phase 1 --- source/source_relax/bfgs_basic.cpp | 63 +++++++------- source/source_relax/bfgs_basic.h | 10 ++- source/source_relax/ions_move_basic.cpp | 50 +++++------ source/source_relax/ions_move_basic.h | 14 +++- source/source_relax/ions_move_bfgs.cpp | 40 ++++----- source/source_relax/ions_move_bfgs.h | 8 +- source/source_relax/ions_move_cg.cpp | 14 ++-- source/source_relax/ions_move_cg.h | 4 +- source/source_relax/ions_move_methods.cpp | 11 +-- source/source_relax/ions_move_methods.h | 5 +- source/source_relax/ions_move_sd.cpp | 10 +-- source/source_relax/ions_move_sd.h | 4 +- source/source_relax/relax_nsync.cpp | 7 +- source/source_relax/relax_nsync.h | 4 +- source/source_relax/test/bfgs_basic_test.cpp | 49 +++++++---- .../test/ions_move_basic_test.cpp | 65 ++++++++------- .../source_relax/test/ions_move_bfgs_test.cpp | 78 ++++++++++-------- .../source_relax/test/ions_move_cg_test.cpp | 82 ++++++++++++------- .../source_relax/test/ions_move_sd_test.cpp | 20 ++--- 19 files changed, 306 insertions(+), 232 deletions(-) diff --git a/source/source_relax/bfgs_basic.cpp b/source/source_relax/bfgs_basic.cpp index 36fc12f0a15..54a8c889a1b 100644 --- a/source/source_relax/bfgs_basic.cpp +++ b/source/source_relax/bfgs_basic.cpp @@ -32,7 +32,7 @@ void BFGS_Basic::allocate_basic(void) 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); @@ -57,8 +57,8 @@ 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; } @@ -96,7 +96,7 @@ 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.data(), move_p.data(), dim); double dot = dot_func(grad.data(), move_p.data(), dim); @@ -111,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; @@ -135,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; } @@ -175,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"); @@ -203,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); } //-------------------------------------------------------------------- @@ -227,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.data(), move.data(), 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) { @@ -244,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]; @@ -264,7 +265,7 @@ 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 ="<<" "<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); @@ -320,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) @@ -336,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++) { diff --git a/source/source_relax/bfgs_basic.h b/source/source_relax/bfgs_basic.h index 6e74cac9cdd..8d101bce1e4 100644 --- a/source/source_relax/bfgs_basic.h +++ b/source/source_relax/bfgs_basic.h @@ -1,6 +1,8 @@ #ifndef BFGS_BASIC #define BFGS_BASIC +#include +#include #include "source_base/matrix.h" #include @@ -23,7 +25,7 @@ class BFGS_Basic 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); @@ -53,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 diff --git a/source/source_relax/ions_move_basic.cpp b/source/source_relax/ions_move_basic.cpp index 08b45d29ce8..03157644e56 100644 --- a/source/source_relax/ions_move_basic.cpp +++ b/source/source_relax/ions_move_basic.cpp @@ -18,7 +18,7 @@ double Ions_Move_Basic::best_xxx = 1.0; int Ions_Move_Basic::out_stru = 0; std::vector Ions_Move_Basic::relax_method = {"bfgs","2"}; -void Ions_Move_Basic::setup_gradient(const UnitCell &ucell, const ModuleBase::matrix &force, double *pos, double *grad) +void Ions_Move_Basic::setup_gradient(const UnitCell &ucell, const ModuleBase::matrix &force, double *pos, double *grad, std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_Basic", "setup_gradient"); @@ -55,7 +55,7 @@ void Ions_Move_Basic::setup_gradient(const UnitCell &ucell, const ModuleBase::ma return; } -void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos) +void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos, std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_Basic", "move_atoms"); @@ -68,8 +68,8 @@ void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos) if (PARAM.inp.test_relax_method) { int iat = 0; - GlobalV::ofs_running << "\n movement of ions (unit is Bohr) : " << std::endl; - GlobalV::ofs_running << " " << std::setw(12) << "Atom" << std::setw(15) << "x" << std::setw(15) << "y" + ofs << "\n movement of ions (unit is Bohr) : " << std::endl; + ofs << " " << std::setw(12) << "Atom" << std::setw(15) << "x" << std::setw(15) << "y" << std::setw(15) << "z" << std::endl; for (int it = 0; it < ucell.ntype; it++) { @@ -77,7 +77,7 @@ void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos) { std::stringstream ss; ss << "move_" << ucell.atoms[it].label << ia + 1; - GlobalV::ofs_running << " " << std::setw(12) << ss.str().c_str() << std::setw(15) << move[3 * iat + 0] + ofs << " " << std::setw(12) << ss.str().c_str() << std::setw(15) << move[3 * iat + 0] << std::setw(15) << move[3 * iat + 1] << std::setw(15) << move[3 * iat + 2] << std::endl; iat++; @@ -106,12 +106,12 @@ void Ions_Move_Basic::move_atoms(UnitCell &ucell, double *move, double *pos) //-------------------------------------------- // Print out the structure file. //-------------------------------------------- - unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,GlobalV::ofs_running); + unitcell::print_tau(ucell.atoms,ucell.Coordinate,ucell.ntype,ucell.lat0,ofs); return; } -void Ions_Move_Basic::check_converged(const UnitCell &ucell, const double *grad) +void Ions_Move_Basic::check_converged(const UnitCell &ucell, const double *grad, std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_Basic", "check_converged"); assert(dim > 0); @@ -132,10 +132,10 @@ void Ions_Move_Basic::check_converged(const UnitCell &ucell, const double *grad) if (PARAM.inp.test_relax_method) { - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "old total energy (ry)", etot_p); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "new total energy (ry)", etot); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "energy difference (ry)", Ions_Move_Basic::ediff); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "largest gradient (ry/bohr)", Ions_Move_Basic::largest_grad); + ModuleBase::GlobalFunc::OUT(ofs, "old total energy (ry)", etot_p); + ModuleBase::GlobalFunc::OUT(ofs, "new total energy (ry)", etot); + ModuleBase::GlobalFunc::OUT(ofs, "energy difference (ry)", Ions_Move_Basic::ediff); + ModuleBase::GlobalFunc::OUT(ofs, "largest gradient (ry/bohr)", Ions_Move_Basic::largest_grad); } if (PARAM.inp.out_level == "ie") @@ -144,7 +144,7 @@ void Ions_Move_Basic::check_converged(const UnitCell &ucell, const double *grad) std::cout << " LARGEST GRAD (eV/Angstrom) : " << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << std::endl; - GlobalV::ofs_running << "\n Largest force is " << largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A + ofs << "\n Largest force is " << largest_grad * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << " eV/Angstrom while threshold is " << PARAM.inp.force_thr_ev << " eV/Angstrom" << std::endl; } @@ -156,22 +156,22 @@ void Ions_Move_Basic::check_converged(const UnitCell &ucell, const double *grad) if (Ions_Move_Basic::largest_grad == 0.0) { - GlobalV::ofs_running << " largest force is 0, no movement is possible." << std::endl; - GlobalV::ofs_running << " it may converged, otherwise no movement of atom is allowed." << std::endl; + ofs << " largest force is 0, no movement is possible." << std::endl; + ofs << " it may converged, otherwise no movement of atom is allowed." << std::endl; Ions_Move_Basic::converged = true; } // mohan update 2011-04-21 else if (etot_diff < etot_thr && Ions_Move_Basic::largest_grad < PARAM.inp.force_thr ) { - GlobalV::ofs_running << "\n Ion relaxation is converged!" << std::endl; - GlobalV::ofs_running << "\n Energy difference (Ry) = " << etot_diff << std::endl; + ofs << "\n Ion relaxation is converged!" << std::endl; + ofs << "\n Energy difference (Ry) = " << etot_diff << std::endl; Ions_Move_Basic::converged = true; ++Ions_Move_Basic::update_iter; } else { - GlobalV::ofs_running << "\n Ion relaxation is not converged yet (threshold is " + ofs << "\n Ion relaxation is not converged yet (threshold is " << PARAM.inp.force_thr * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << ")" << std::endl; // std::cout << "\n etot_diff=" << etot_diff << " etot_thr=" << etot_thr //<< " largest_grad=" << largest_grad << " force_thr=" << PARAM.inp.force_thr << std::endl; @@ -181,16 +181,16 @@ void Ions_Move_Basic::check_converged(const UnitCell &ucell, const double *grad) return; } -void Ions_Move_Basic::terminate(const UnitCell &ucell) +void Ions_Move_Basic::terminate(const UnitCell &ucell, std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_Basic", "terminate"); if (Ions_Move_Basic::converged) { - GlobalV::ofs_running << " end of geometry optimization" << std::endl; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "istep", Ions_Move_Basic::istep); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "update iteration", Ions_Move_Basic::update_iter); + ofs << " end of geometry optimization" << std::endl; + ModuleBase::GlobalFunc::OUT(ofs, "istep", Ions_Move_Basic::istep); + ModuleBase::GlobalFunc::OUT(ofs, "update iteration", Ions_Move_Basic::update_iter); /* - GlobalV::ofs_running<<"Saving the approximate inverse hessian"< +#include #include "relax_data.h" #include "source_base/matrix.h" #include "source_cell/unitcell.h" @@ -42,29 +44,33 @@ extern int out_stru; ///< Structure output flag * @param force Force matrix (nat x 3) * @param pos Output position array (dimension: dim) * @param grad Output gradient array (dimension: dim) + * @param ofs Output stream for logging */ -void setup_gradient(const UnitCell &ucell, const ModuleBase::matrix &force, double *pos, double *grad); +void setup_gradient(const UnitCell &ucell, const ModuleBase::matrix &force, double *pos, double *grad, std::ofstream& ofs); /** * @brief Move atoms according to displacement vector. * @param ucell Unit cell to update * @param move Displacement vector (dimension: dim) * @param pos Current position array (dimension: dim) + * @param ofs Output stream for logging */ -void move_atoms(UnitCell &ucell, double *move, double *pos); +void move_atoms(UnitCell &ucell, double *move, double *pos, std::ofstream& ofs); /** * @brief Check convergence based on gradient threshold. * @param ucell Unit cell containing lattice information * @param grad Gradient array (dimension: dim) + * @param ofs Output stream for logging */ -void check_converged(const UnitCell &ucell, const double *grad); +void check_converged(const UnitCell &ucell, const double *grad, std::ofstream& ofs); /** * @brief Terminate geometry optimization and output results. * @param ucell Unit cell to output + * @param ofs Output stream for logging */ -void terminate(const UnitCell &ucell); +void terminate(const UnitCell &ucell, std::ofstream& ofs); /** * @brief Update energy values and compute energy difference. diff --git a/source/source_relax/ions_move_bfgs.cpp b/source/source_relax/ions_move_bfgs.cpp index 11b712805ba..3330e874aac 100644 --- a/source/source_relax/ions_move_bfgs.cpp +++ b/source/source_relax/ions_move_bfgs.cpp @@ -37,7 +37,7 @@ void Ions_Move_BFGS::allocate() return; } -void Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, const double& energy_in) +void Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, const double& energy_in, std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_BFGS", "start"); @@ -48,23 +48,23 @@ void Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, con // In the following steps, the pos is updated by BFGS methods. if (first_step) { - Ions_Move_Basic::setup_gradient(ucell, force, this->pos.data(), this->grad.data()); + Ions_Move_Basic::setup_gradient(ucell, force, this->pos.data(), this->grad.data(), ofs); first_step = false; } else { std::vector pos_tmp(3 * ucell.nat); - Ions_Move_Basic::setup_gradient(ucell, force, pos_tmp.data(), this->grad.data()); + Ions_Move_Basic::setup_gradient(ucell, force, pos_tmp.data(), this->grad.data(), ofs); } // use energy_in and istep to setup etot and etot_old. Ions_Move_Basic::setup_etot(energy_in, false); // use gradient and etot and etot_old to check // if the result is converged. - Ions_Move_Basic::check_converged(ucell, this->grad.data()); + Ions_Move_Basic::check_converged(ucell, this->grad.data(), ofs); if (Ions_Move_Basic::converged) { - Ions_Move_Basic::terminate(ucell); + Ions_Move_Basic::terminate(ucell, ofs); } else { @@ -73,7 +73,7 @@ void Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, con // [ if run from previous saved info ] // the BFGS file is read from previous run. // and the move_p is renormalized. - this->restart_bfgs(ucell.lat0); + this->restart_bfgs(ucell.lat0, ofs); //[ if etot>etot_p ] // interpolation @@ -81,18 +81,18 @@ void Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, con // calculate the new step -> the new move using hessian // matrix, and set the new trust radius. // [compute the move at last] - this->bfgs_routine(ucell.lat0); + this->bfgs_routine(ucell.lat0, ofs); // get prepared for the next try. // even if the energy is higher, we save the information. this->save_bfgs(); - Ions_Move_Basic::move_atoms(ucell, move.data(), pos.data()); + Ions_Move_Basic::move_atoms(ucell, move.data(), pos.data(), ofs); } return; } -void Ions_Move_BFGS::restart_bfgs(const double& lat0) +void Ions_Move_BFGS::restart_bfgs(const double& lat0, std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_BFGS", "restart_bfgs"); @@ -114,7 +114,7 @@ void Ions_Move_BFGS::restart_bfgs(const double& lat0) if (PARAM.inp.test_relax_method) { - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "trust_radius_old (bohr)", trust_radius_old); + ModuleBase::GlobalFunc::OUT(ofs, "trust_radius_old (bohr)", trust_radius_old); } // (2) @@ -168,7 +168,7 @@ void Ions_Move_BFGS::restart_bfgs(const double& lat0) return; } -void Ions_Move_BFGS::bfgs_routine(const double& lat0) +void Ions_Move_BFGS::bfgs_routine(const double& lat0, std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_BFGS", "bfgs_routine"); using namespace Ions_Move_Basic; @@ -229,9 +229,9 @@ void Ions_Move_BFGS::bfgs_routine(const double& lat0) if (PARAM.inp.test_relax_method) { - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "dE0s", dE0s); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "den", den); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "interpolated trust radius", trust_radius); + ModuleBase::GlobalFunc::OUT(ofs, "dE0s", dE0s); + ModuleBase::GlobalFunc::OUT(ofs, "den", den); + ModuleBase::GlobalFunc::OUT(ofs, "interpolated trust radius", trust_radius); } // std::cout << " Formula : " << etot << " * s^2 + " << dE0s << " * s + " << etot_p << std::endl; // std::cout << " Lowest point : " << trust_radius << std::endl; @@ -242,7 +242,7 @@ void Ions_Move_BFGS::bfgs_routine(const double& lat0) // then do is again, but smaller raidus. trust_radius = 0.5 * trust_radius_old; - GlobalV::ofs_running << " quadratic interpolation is impossible." << std::endl; + ofs << " quadratic interpolation is impossible." << std::endl; } // values from the last succeseful bfgs step are restored etot = etot_p; @@ -257,11 +257,11 @@ void Ions_Move_BFGS::bfgs_routine(const double& lat0) // we are trapped in this case..., so the algorithim must be restart // the history is reset // xiaohui add 2013-03-17 - GlobalV::ofs_running << "trust_radius = " << trust_radius << std::endl; - GlobalV::ofs_running << "relax_bfgs_rmin = " << relax_bfgs_rmin << std::endl; - GlobalV::ofs_running << "relax_bfgs_rmax = " << relax_bfgs_rmax << std::endl; + ofs << "trust_radius = " << trust_radius << std::endl; + ofs << "relax_bfgs_rmin = " << relax_bfgs_rmin << std::endl; + ofs << "relax_bfgs_rmax = " << relax_bfgs_rmax << std::endl; // xiaohui add 2013-03-17 - GlobalV::ofs_running << " trust_radius < relax_bfgs_rmin, reset bfgs history." << std::endl; + ofs << " trust_radius < relax_bfgs_rmin, reset bfgs history." << std::endl; if (tr_min_hit) { @@ -293,7 +293,7 @@ void Ions_Move_BFGS::bfgs_routine(const double& lat0) } else if (etot <= etot_p) { - this->new_step(lat0); + this->new_step(lat0, ofs); } if (PARAM.inp.out_level == "ie") diff --git a/source/source_relax/ions_move_bfgs.h b/source/source_relax/ions_move_bfgs.h index a9c030b0200..6d75868165a 100644 --- a/source/source_relax/ions_move_bfgs.h +++ b/source/source_relax/ions_move_bfgs.h @@ -1,6 +1,8 @@ #ifndef IONS_MOVE_BFGS_H #define IONS_MOVE_BFGS_H +#include +#include #include "bfgs_basic.h" #include "source_base/matrix.h" #include "source_cell/unitcell.h" @@ -11,12 +13,12 @@ class Ions_Move_BFGS : public BFGS_Basic ~Ions_Move_BFGS(); void allocate(void); - void start(UnitCell& ucell, const ModuleBase::matrix& force, const double& energy_in); + void start(UnitCell& ucell, const ModuleBase::matrix& force, const double& energy_in, std::ofstream& ofs); private: bool init_done; - void bfgs_routine(const double& lat0); - void restart_bfgs(const double& lat0); + void bfgs_routine(const double& lat0, std::ofstream& ofs); + void restart_bfgs(const double& lat0, std::ofstream& ofs); bool first_step=true; // If it is the first step of the relaxation. The pos is only generated from ucell in the first step, and in the following steps, the pos is generated from the previous step. }; diff --git a/source/source_relax/ions_move_cg.cpp b/source/source_relax/ions_move_cg.cpp index 00a73f648d6..432d340d86e 100644 --- a/source/source_relax/ions_move_cg.cpp +++ b/source/source_relax/ions_move_cg.cpp @@ -25,7 +25,7 @@ void Ions_Move_CG::allocate(void) this->e0 = 0.0; } -void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const double &etot_in) +void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const double &etot_in, std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_CG", "start"); assert(dim > 0); @@ -72,16 +72,16 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const nbrent = 0; } - Ions_Move_Basic::setup_gradient(ucell, force, pos.data(), grad.data()); + Ions_Move_Basic::setup_gradient(ucell, force, pos.data(), grad.data(), ofs); Ions_Move_Basic::setup_etot(etot_in, 0); if (flag == 0) { - Ions_Move_Basic::check_converged(ucell, grad.data()); + Ions_Move_Basic::check_converged(ucell, grad.data(), ofs); } if (Ions_Move_Basic::converged) { - Ions_Move_Basic::terminate(ucell); + Ions_Move_Basic::terminate(ucell, ofs); break; } @@ -93,7 +93,7 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const CG_Base::normalize(dim, cg_gradn.data(), cg_grad.data()); CG_Base::setup_move(dim, move0.data(), cg_gradn.data(), steplength); - Ions_Move_Basic::move_atoms(ucell, move0.data(), pos.data()); + Ions_Move_Basic::move_atoms(ucell, move0.data(), pos.data(), ofs); for (int i = 0; i < dim; i++) { @@ -145,7 +145,7 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const } CG_Base::setup_move(dim, move.data(), cg_gradn.data(), best_x); - Ions_Move_Basic::move_atoms(ucell, move.data(), pos.data()); + Ions_Move_Basic::move_atoms(ucell, move.data(), pos.data(), ofs); trial = false; xa = 0; CG_Base::f_cal(dim, move0.data(), move.data(), xc); @@ -183,7 +183,7 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const CG_Base::normalize(dim, cg_gradn.data(), cg_grad0.data()); CG_Base::setup_move(dim, move.data(), cg_gradn.data(), best_x); - Ions_Move_Basic::move_atoms(ucell, move.data(), pos.data()); + Ions_Move_Basic::move_atoms(ucell, move.data(), pos.data(), ofs); Ions_Move_Basic::relax_bfgs_init = xc; break; } diff --git a/source/source_relax/ions_move_cg.h b/source/source_relax/ions_move_cg.h index d846219e8ec..bd1366c8d01 100644 --- a/source/source_relax/ions_move_cg.h +++ b/source/source_relax/ions_move_cg.h @@ -1,6 +1,8 @@ #ifndef IONS_MOVE_CG_H #define IONS_MOVE_CG_H +#include +#include #include "source_base/matrix.h" #include "source_cell/unitcell.h" #include "cg_base.h" @@ -13,7 +15,7 @@ class Ions_Move_CG : public CG_Base ~Ions_Move_CG() = default; void allocate(void); - void start(UnitCell &ucell, const ModuleBase::matrix &force, const double &etot); + void start(UnitCell &ucell, const ModuleBase::matrix &force, const double &etot, std::ofstream& ofs); static double RELAX_CG_THR; int sd_step = 0; diff --git a/source/source_relax/ions_move_methods.cpp b/source/source_relax/ions_move_methods.cpp index db485f4c09e..d54f2b1b4c0 100644 --- a/source/source_relax/ions_move_methods.cpp +++ b/source/source_relax/ions_move_methods.cpp @@ -53,7 +53,8 @@ void Ions_Move_Methods::cal_movement(const int &istep, const int &force_step, const ModuleBase::matrix &f, const double &etot, - UnitCell &ucell) + UnitCell &ucell, + std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_Methods", "init"); // Ions_Move_Basic::istep = istep; @@ -63,19 +64,19 @@ void Ions_Move_Methods::cal_movement(const int &istep, // move_ions // output tau // check all symmery - bfgs.start(ucell, f, etot); + bfgs.start(ucell, f, etot, ofs); } else if (Ions_Move_Basic::relax_method[0] == "sd") { - sd.start(ucell, f, etot); + sd.start(ucell, f, etot, ofs); } else if (Ions_Move_Basic::relax_method[0] == "cg") { - cg.start(ucell, f, etot); + cg.start(ucell, f, etot, ofs); } else if (Ions_Move_Basic::relax_method[0] == "cg_bfgs") { - cg.start(ucell, f, etot); // added by pengfei 13-8-10 + cg.start(ucell, f, etot, ofs); // added by pengfei 13-8-10 } else if(Ions_Move_Basic::relax_method[0] == "bfgs"&&Ions_Move_Basic::relax_method[1] == "1") { diff --git a/source/source_relax/ions_move_methods.h b/source/source_relax/ions_move_methods.h index b9b4e28611e..a7b4ed84cc9 100644 --- a/source/source_relax/ions_move_methods.h +++ b/source/source_relax/ions_move_methods.h @@ -1,6 +1,8 @@ #ifndef IONS_MOVE_METHODS_H #define IONS_MOVE_METHODS_H +#include +#include #include "ions_move_basic.h" #include "ions_move_bfgs.h" #include "ions_move_cg.h" @@ -20,7 +22,8 @@ class Ions_Move_Methods const int &force_step, const ModuleBase::matrix &f, const double &etot, - UnitCell &ucell); + UnitCell &ucell, + std::ofstream& ofs); bool get_converged() const { diff --git a/source/source_relax/ions_move_sd.cpp b/source/source_relax/ions_move_sd.cpp index 2e0f4f7a322..8405d572d88 100644 --- a/source/source_relax/ions_move_sd.cpp +++ b/source/source_relax/ions_move_sd.cpp @@ -20,7 +20,7 @@ void Ions_Move_SD::allocate() pos_saved.resize(dim, 0.0); } -void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const double& etot_in) +void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const double& etot_in, std::ofstream& ofs) { ModuleBase::TITLE("Ions_Move_SD", "start"); @@ -36,7 +36,7 @@ void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const // 0: ediff < 0 bool judgement = false; setup_etot(etot_in, judgement); - setup_gradient(ucell, force, pos.data(), grad.data()); + setup_gradient(ucell, force, pos.data(), grad.data(), ofs); if (istep == 1 || etot_in <= energy_saved) { @@ -61,10 +61,10 @@ void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const } } - Ions_Move_Basic::check_converged(ucell, grad.data()); + Ions_Move_Basic::check_converged(ucell, grad.data(), ofs); if (Ions_Move_Basic::converged) { - Ions_Move_Basic::terminate(ucell); + Ions_Move_Basic::terminate(ucell, ofs); } else { @@ -73,7 +73,7 @@ void Ions_Move_SD::start(UnitCell& ucell, const ModuleBase::matrix& force, const { move[i] = -grad_saved[i] * trust_radius; } - move_atoms(ucell, move.data(), pos_saved.data()); + move_atoms(ucell, move.data(), pos_saved.data(), ofs); Ions_Move_Basic::update_iter++; } diff --git a/source/source_relax/ions_move_sd.h b/source/source_relax/ions_move_sd.h index e2eb53d587c..69b3c30e6a9 100644 --- a/source/source_relax/ions_move_sd.h +++ b/source/source_relax/ions_move_sd.h @@ -1,6 +1,8 @@ #ifndef IONS_MOVE_SD_H #define IONS_MOVE_SD_H +#include +#include #include "source_base/matrix.h" #include "source_cell/unitcell.h" #include @@ -12,7 +14,7 @@ class Ions_Move_SD ~Ions_Move_SD() = default; void allocate(void); - void start(UnitCell& ucell, const ModuleBase::matrix& force, const double& etot); + void start(UnitCell& ucell, const ModuleBase::matrix& force, const double& etot, std::ofstream& ofs); private: double energy_saved; diff --git a/source/source_relax/relax_nsync.cpp b/source/source_relax/relax_nsync.cpp index 9541c9dc177..81f64e47bf7 100644 --- a/source/source_relax/relax_nsync.cpp +++ b/source/source_relax/relax_nsync.cpp @@ -49,7 +49,7 @@ bool Relax_old::relax_step(const int& istep, { // do relax calculation and generate next structure bool converged = false; - converged = this->do_relax(istep, force, energy, ucell, force_step); + converged = this->do_relax(istep, force, energy, ucell, force_step, GlobalV::ofs_running); if (!converged) { ucell.ionic_position_updated = true; @@ -129,10 +129,11 @@ bool Relax_old::do_relax(const int& istep, const ModuleBase::matrix& ionic_force, const double& total_energy, UnitCell& ucell, - int& jstep) + int& jstep, + std::ofstream& ofs) { ModuleBase::TITLE("Relax_old", "do_relax"); - IMM.cal_movement(istep, jstep, ionic_force, total_energy, ucell); + IMM.cal_movement(istep, jstep, ionic_force, total_energy, ucell, ofs); ++jstep; return IMM.get_converged(); } diff --git a/source/source_relax/relax_nsync.h b/source/source_relax/relax_nsync.h index 33eca958020..d98525d536b 100644 --- a/source/source_relax/relax_nsync.h +++ b/source/source_relax/relax_nsync.h @@ -1,6 +1,7 @@ #ifndef RELAX_OLD_H #define RELAX_OLD_H +#include #include "ions_move_methods.h" #include "lattice_change_methods.h" #include "source_cell/unitcell.h" @@ -28,7 +29,8 @@ class Relax_old const ModuleBase::matrix& ionic_force, const double& total_energy, UnitCell& ucell, - int& jstep); + int& jstep, + std::ofstream& ofs); bool do_cellrelax(const int& istep, const int& stress_step, const ModuleBase::matrix& stress, diff --git a/source/source_relax/test/bfgs_basic_test.cpp b/source/source_relax/test/bfgs_basic_test.cpp index efd30e9829e..184632a8976 100644 --- a/source/source_relax/test/bfgs_basic_test.cpp +++ b/source/source_relax/test/bfgs_basic_test.cpp @@ -57,7 +57,10 @@ TEST_F(BFGSBasicTest, UpdateInverseHessianDeath) { Ions_Move_Basic::dim = 0; double lat0 = 1.0; - ASSERT_DEATH(bfgs.update_inverse_hessian(lat0), ""); + std::ofstream ofs("test_log_update_inverse_hessian_death.log"); + ASSERT_DEATH(bfgs.update_inverse_hessian(lat0, ofs), ""); + ofs.close(); + std::remove("test_log_update_inverse_hessian_death.log"); } // Test function update_inverse_hessian() when sdoty = 0 @@ -67,18 +70,18 @@ TEST_F(BFGSBasicTest, UpdateInverseHessianCase1) double lat0 = 1.0; bfgs.allocate_basic(); - GlobalV::ofs_running.open("log"); - bfgs.update_inverse_hessian(lat0); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_log_update_inverse_hessian_case1.log"); + bfgs.update_inverse_hessian(lat0, ofs); + ofs.close(); std::string expected_output = " WARINIG: unexpected behaviour in update_inverse_hessian\n Resetting bfgs history \n"; - std::ifstream ifs("log"); + std::ifstream ifs("test_log_update_inverse_hessian_case1.log"); std::string output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); + std::remove("test_log_update_inverse_hessian_case1.log"); EXPECT_EQ(expected_output, output); - std::remove("log"); } // Test function update_inverse_hessian() @@ -90,7 +93,9 @@ TEST_F(BFGSBasicTest, UpdateInverseHessianCase2) bfgs.pos[0] = 2.0; bfgs.grad[0] = 2.0; - bfgs.update_inverse_hessian(lat0); + std::ofstream ofs("test_log_update_inverse_hessian_case2.log"); + bfgs.update_inverse_hessian(lat0, ofs); + ofs.close(); EXPECT_DOUBLE_EQ(bfgs.inv_hess(0, 0), 0.5); EXPECT_DOUBLE_EQ(bfgs.inv_hess(0, 1), 0.0); @@ -101,6 +106,7 @@ TEST_F(BFGSBasicTest, UpdateInverseHessianCase2) EXPECT_DOUBLE_EQ(bfgs.inv_hess(2, 0), 0.0); EXPECT_DOUBLE_EQ(bfgs.inv_hess(2, 1), 0.0); EXPECT_DOUBLE_EQ(bfgs.inv_hess(2, 2), 0.0); + std::remove("test_log_update_inverse_hessian_case2.log"); } // Test function check_wolfe_conditions() @@ -114,9 +120,9 @@ TEST_F(BFGSBasicTest, CheckWolfeConditions) bfgs.grad[0] = 2.0; bfgs.move[0] = 1.0; - GlobalV::ofs_running.open("log"); - bfgs.check_wolfe_conditions(); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_log_check_wolfe_conditions.log"); + bfgs.check_wolfe_conditions(ofs); + ofs.close(); std::string expected_output = " etot - etot_p = 10\n relax_bfgs_w1 * dot_p = -0\n " @@ -125,13 +131,14 @@ TEST_F(BFGSBasicTest, CheckWolfeConditions) "wolfe1 = 0\n wolfe2 = 0\n etot - etot_p = 10\n " " relax_bfgs_w1 * dot_p = -0\n wolfe1 = 0\n " " wolfe2 = 0\n wolfe condition satisfied = 0\n"; - std::ifstream ifs("log"); + + std::ifstream ifs("test_log_check_wolfe_conditions.log"); std::string output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); + std::remove("test_log_check_wolfe_conditions.log"); EXPECT_EQ(bfgs.wolfe_flag, false); EXPECT_EQ(expected_output, output); - std::remove("log"); } // Test function reset_hessian() @@ -195,7 +202,8 @@ TEST_F(BFGSBasicTest, NewStepCase1) bfgs.inv_hess(1, 1) = -6.0; double lat0 = 1.0; - bfgs.new_step(lat0); + std::ofstream ofs("test_log.log"); + bfgs.new_step(lat0, ofs); EXPECT_EQ(Ions_Move_Basic::update_iter, 1); EXPECT_EQ(bfgs.tr_min_hit, false); @@ -228,7 +236,8 @@ TEST_F(BFGSBasicTest, NewStepCase2) bfgs.inv_hess(1, 1) = -6.0; double lat0 = 1.0; - bfgs.new_step(lat0); + std::ofstream ofs("test_log.log"); + bfgs.new_step(lat0, ofs); EXPECT_EQ(Ions_Move_Basic::update_iter, 3); EXPECT_DOUBLE_EQ(Ions_Move_Basic::trust_radius, -1.0); @@ -247,9 +256,10 @@ TEST_F(BFGSBasicTest, NewStepWarningQuit) bfgs.bfgs_ndim = 2; bfgs.allocate_basic(); double lat0 = 1.0; + std::ofstream ofs("test_log.log"); testing::internal::CaptureStdout(); - EXPECT_EXIT(bfgs.new_step(lat0), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(bfgs.new_step(lat0, ofs), ::testing::ExitedWithCode(1), ""); std::string output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("bfgs_ndim > 1 not implemented yet")); } @@ -270,7 +280,8 @@ TEST_F(BFGSBasicTest, ComputeTrustRadiusCase1) bfgs.wolfe_flag = true; bfgs.relax_bfgs_w1 = 1.0; - bfgs.compute_trust_radius(); + std::ofstream ofs("test_log.log"); + bfgs.compute_trust_radius(ofs); EXPECT_EQ(bfgs.tr_min_hit, false); EXPECT_DOUBLE_EQ(Ions_Move_Basic::trust_radius, -1.0); @@ -303,7 +314,8 @@ TEST_F(BFGSBasicTest, ComputeTrustRadiusCase2) bfgs.relax_bfgs_w1 = 1.0; bfgs.tr_min_hit = false; - bfgs.compute_trust_radius(); + std::ofstream ofs("test_log.log"); + bfgs.compute_trust_radius(ofs); EXPECT_EQ(bfgs.tr_min_hit, true); EXPECT_DOUBLE_EQ(Ions_Move_Basic::trust_radius, 100.0); @@ -336,8 +348,9 @@ TEST_F(BFGSBasicTest, ComputeTrustRadiusWarningQuit) bfgs.relax_bfgs_w1 = 1.0; bfgs.tr_min_hit = true; + std::ofstream ofs("test_log.log"); testing::internal::CaptureStdout(); - EXPECT_EXIT(bfgs.compute_trust_radius(), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(bfgs.compute_trust_radius(ofs), ::testing::ExitedWithCode(1), ""); std::string output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("bfgs history already reset at previous step, we got trapped!")); } \ No newline at end of file diff --git a/source/source_relax/test/ions_move_basic_test.cpp b/source/source_relax/test/ions_move_basic_test.cpp index b6f02628aec..dd980d57f5a 100644 --- a/source/source_relax/test/ions_move_basic_test.cpp +++ b/source/source_relax/test/ions_move_basic_test.cpp @@ -42,7 +42,10 @@ TEST_F(IonsMoveBasicTest, SetupGradient) { // Call the function being tested Ions_Move_Basic::dim = 6; - Ions_Move_Basic::setup_gradient(ucell, force, pos, grad); + std::ofstream ofs("test_setup_gradient.log"); + Ions_Move_Basic::setup_gradient(ucell, force, pos, grad, ofs); + ofs.close(); + std::remove("test_setup_gradient.log"); // Check that the expected positions and gradients were generated EXPECT_DOUBLE_EQ(pos[0], 0.0); @@ -72,18 +75,18 @@ TEST_F(IonsMoveBasicTest, MoveAtoms) } // Call the function being tested - GlobalV::ofs_running.open("log"); - Ions_Move_Basic::move_atoms(ucell, move, pos); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_move_atoms.log"); + Ions_Move_Basic::move_atoms(ucell, move, pos, ofs); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_move_atoms.log"); std::string expected_output = "\n movement of ions (unit is Bohr) : \n Atom x y " " z\n move_1 0 1 2\n " "move_2 3 4 5\n"; std::string output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_move_atoms.log"); EXPECT_THAT(output , ::testing::HasSubstr(expected_output)); EXPECT_DOUBLE_EQ(pos[0], 0.0); @@ -108,17 +111,17 @@ TEST_F(IonsMoveBasicTest, CheckConvergedCase1) } // Call the function being tested - GlobalV::ofs_running.open("log"); + std::ofstream ofs("test_check_converged_case1.log"); testing::internal::CaptureStdout(); - Ions_Move_Basic::check_converged(ucell, grad); + Ions_Move_Basic::check_converged(ucell, grad, ofs); std::string std_outout = testing::internal::GetCapturedStdout(); - GlobalV::ofs_running.close(); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_check_converged_case1.log"); std::string ofs_output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_check_converged_case1.log"); std::string expected_ofs = " old total energy (ry) = 0\n new total energy (ry) = 0\n " @@ -147,17 +150,17 @@ TEST_F(IonsMoveBasicTest, CheckConvergedCase2) grad[0] = 1.0; // Call the function being tested - GlobalV::ofs_running.open("log"); + std::ofstream ofs("test_check_converged_case2.log"); testing::internal::CaptureStdout(); - Ions_Move_Basic::check_converged(ucell, grad); + Ions_Move_Basic::check_converged(ucell, grad, ofs); std::string std_outout = testing::internal::GetCapturedStdout(); - GlobalV::ofs_running.close(); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_check_converged_case2.log"); std::string ofs_output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_check_converged_case2.log"); std::string expected_ofs = " old total energy (ry) = 0\n new total energy (ry) = 0\n " @@ -186,17 +189,17 @@ TEST_F(IonsMoveBasicTest, CheckConvergedCase3) grad[0] = 1.0; // Call the function being tested - GlobalV::ofs_running.open("log"); + std::ofstream ofs("test_check_converged_case3.log"); testing::internal::CaptureStdout(); - Ions_Move_Basic::check_converged(ucell, grad); + Ions_Move_Basic::check_converged(ucell, grad, ofs); std::string std_outout = testing::internal::GetCapturedStdout(); - GlobalV::ofs_running.close(); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_check_converged_case3.log"); std::string ofs_output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_check_converged_case3.log"); std::string expected_ofs = " old total energy (ry) = 0\n new total energy (ry) = 0\n " @@ -221,15 +224,15 @@ TEST_F(IonsMoveBasicTest, TerminateConverged) Ions_Move_Basic::update_iter = 5; // Call the function being tested - GlobalV::ofs_running.open("log"); - Ions_Move_Basic::terminate(ucell); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_terminate_converged.log"); + Ions_Move_Basic::terminate(ucell, ofs); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_terminate_converged.log"); std::string ofs_output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_terminate_converged.log"); std::string expected_ofs = " end of geometry optimization\n istep = 2\n " " update iteration = 5\n"; @@ -244,15 +247,15 @@ TEST_F(IonsMoveBasicTest, TerminateNotConverged) Ions_Move_Basic::converged = false; // Call the function being tested - GlobalV::ofs_running.open("log"); - Ions_Move_Basic::terminate(ucell); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_terminate_not_converged.log"); + Ions_Move_Basic::terminate(ucell, ofs); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_terminate_not_converged.log"); std::string ofs_output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_terminate_not_converged.log"); std::string expected_ofs = " the maximum number of steps has been reached.\n end of geometry optimization.\n"; diff --git a/source/source_relax/test/ions_move_bfgs_test.cpp b/source/source_relax/test/ions_move_bfgs_test.cpp index f13ef8c3ea5..8ad773bbf55 100644 --- a/source/source_relax/test/ions_move_bfgs_test.cpp +++ b/source/source_relax/test/ions_move_bfgs_test.cpp @@ -73,15 +73,15 @@ TEST_F(IonsMoveBFGSTest, StartCase1) // Call the function being tested bfgs.allocate(); - GlobalV::ofs_running.open("log"); - bfgs.start(ucell, force, energy_in); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_start_case1.log"); + bfgs.start(ucell, force, energy_in, ofs); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_start_case1.log"); std::string output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_start_case1.log"); EXPECT_THAT(output, testing::HasSubstr("update iteration")); } @@ -99,15 +99,18 @@ TEST_F(IonsMoveBFGSTest, StartCase2) // Call the function being tested bfgs.allocate(); - GlobalV::ofs_running.open("log"); - EXPECT_EXIT(bfgs.start(ucell, force, energy_in) , ::testing::ExitedWithCode(1), ""); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_start_case2.log"); + testing::internal::CaptureStderr(); + EXPECT_EXIT(bfgs.start(ucell, force, energy_in, ofs) , ::testing::ExitedWithCode(1), ""); + std::string stderr_output = testing::internal::GetCapturedStderr(); + ofs.close(); + std::remove("test_start_case2.log"); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_start_case2.log"); std::string output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_start_case2.log"); EXPECT_THAT(output, testing::HasSubstr("Ion relaxation is not converged yet")); } @@ -129,15 +132,15 @@ TEST_F(IonsMoveBFGSTest, RestartBfgsCase1) } // Call the function being tested - GlobalV::ofs_running.open("log"); - bfgs.restart_bfgs(lat0); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_restart_bfgs_case1.log"); + bfgs.restart_bfgs(lat0, ofs); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_restart_bfgs_case1.log"); std::string output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_restart_bfgs_case1.log"); std::string expected_output = " trust_radius_old (bohr) = 2.44949\n"; @@ -167,7 +170,10 @@ TEST_F(IonsMoveBFGSTest, RestartBfgsCase2) } // Call the function being tested - bfgs.restart_bfgs(lat0); + std::ofstream ofs("test_restart_bfgs_case2.log"); + bfgs.restart_bfgs(lat0, ofs); + ofs.close(); + std::remove("test_restart_bfgs_case2.log"); // Check the results EXPECT_DOUBLE_EQ(Ions_Move_Basic::update_iter, 0.0); @@ -212,17 +218,17 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineCase1) } // Call the function being tested - GlobalV::ofs_running.open("log"); + std::ofstream ofs("test_bfgs_routine_case1.log"); testing::internal::CaptureStdout(); - bfgs.bfgs_routine(lat0); + bfgs.bfgs_routine(lat0, ofs); std::string std_outout = testing::internal::GetCapturedStdout(); - GlobalV::ofs_running.close(); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_bfgs_routine_case1.log"); std::string ofs_output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_bfgs_routine_case1.log"); std::string expected_ofs = " dE0s = 0\n den = 0.1\n " @@ -278,17 +284,17 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineCase2) } // Call the function being tested - GlobalV::ofs_running.open("log"); + std::ofstream ofs("test_bfgs_routine_case2.log"); testing::internal::CaptureStdout(); - bfgs.bfgs_routine(lat0); + bfgs.bfgs_routine(lat0, ofs); std::string std_outout = testing::internal::GetCapturedStdout(); - GlobalV::ofs_running.close(); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_bfgs_routine_case2.log"); std::string ofs_output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_bfgs_routine_case2.log"); std::string expected_ofs = " quadratic interpolation is impossible.\n istep = " "0\n update iteration = 0\n"; @@ -339,15 +345,15 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineCase3) bfgs.inv_hess(1, 1) = -6.0; // Call the function being tested - GlobalV::ofs_running.open("log"); - bfgs.bfgs_routine(lat0); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_bfgs_routine_case3.log"); + bfgs.bfgs_routine(lat0, ofs); + ofs.close(); // Check the results - std::ifstream ifs("log"); + std::ifstream ifs("test_bfgs_routine_case3.log"); std::string ofs_output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_bfgs_routine_case3.log"); std::string expected_ofs = " check the norm of new move 410 (Bohr)\n Uphill move : resetting bfgs history\n " " istep = 0\n update iteration = 1\n"; @@ -397,9 +403,12 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineWarningQuit1) } // Check the results + std::ofstream ofs("test_bfgs_routine_warning_quit1.log"); testing::internal::CaptureStdout(); - EXPECT_EXIT(bfgs.bfgs_routine(lat0), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(bfgs.bfgs_routine(lat0, ofs), ::testing::ExitedWithCode(1), ""); std::string output = testing::internal::GetCapturedStdout(); + ofs.close(); + std::remove("test_bfgs_routine_warning_quit1.log"); EXPECT_THAT(output, testing::HasSubstr("trust radius is too small! Break down.")); } @@ -418,8 +427,11 @@ TEST_F(IonsMoveBFGSTest, BfgsRoutineWarningQuit2) Ions_Move_Basic::relax_bfgs_rmin = 1.0; // Check the results + std::ofstream ofs("test_bfgs_routine_warning_quit2.log"); testing::internal::CaptureStdout(); - EXPECT_EXIT(bfgs.bfgs_routine(lat0), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(bfgs.bfgs_routine(lat0, ofs), ::testing::ExitedWithCode(1), ""); std::string output = testing::internal::GetCapturedStdout(); + ofs.close(); + std::remove("test_bfgs_routine_warning_quit2.log"); EXPECT_THAT(output, testing::HasSubstr("BFGS: move-length unreasonably short")); } diff --git a/source/source_relax/test/ions_move_cg_test.cpp b/source/source_relax/test/ions_move_cg_test.cpp index 83c4323be11..a811e66c9c4 100644 --- a/source/source_relax/test/ions_move_cg_test.cpp +++ b/source/source_relax/test/ions_move_cg_test.cpp @@ -99,9 +99,9 @@ TEST_F(IonsMoveCGTest, TestStartConverged) double etot = 0.0; // call function - GlobalV::ofs_running.open("TestStartConverged.log"); - im_cg.start(ucell, force, etot); - GlobalV::ofs_running.close(); + std::ofstream ofs("TestStartConverged.log"); + im_cg.start(ucell, force, etot, ofs); + ofs.close(); // Check output std::string expected_output = "\n Largest force is 0 eV/Angstrom while threshold is -1 eV/Angstrom\n" @@ -137,9 +137,9 @@ TEST_F(IonsMoveCGTest, TestStartSd) double etot = 0.0; // call function - GlobalV::ofs_running.open("TestStartSd.log"); - im_cg.start(ucell, force, etot); - GlobalV::ofs_running.close(); + std::ofstream ofs("TestStartSd.log"); + im_cg.start(ucell, force, etot, ofs); + ofs.close(); // Check output std::string expected_output = "\n Largest force is 0.257111 eV/Angstrom while threshold is -1 eV/Angstrom\n\n" @@ -172,13 +172,16 @@ TEST_F(IonsMoveCGTest, TestStartTrialGoto) // call function im_cg.move0[0] = 1.0; - im_cg.start(ucell, force, etot); + std::ofstream ofs1("TestStartTrialGoto_temp1.log"); + im_cg.start(ucell, force, etot, ofs1); + ofs1.close(); + std::remove("TestStartTrialGoto_temp1.log"); Ions_Move_Basic::istep = 2; im_cg.move0[0] = 10.0; force(0, 0) = 0.001; - GlobalV::ofs_running.open("TestStartTrialGoto.log"); - im_cg.start(ucell, force, etot); - GlobalV::ofs_running.close(); + std::ofstream ofs("TestStartTrialGoto.log"); + im_cg.start(ucell, force, etot, ofs); + ofs.close(); // Check output std::string expected_output = "\n Largest force is 0.0257111 eV/Angstrom while threshold is -1 eV/Angstrom\n\n" @@ -211,12 +214,15 @@ TEST_F(IonsMoveCGTest, TestStartTrial) // call function im_cg.move0[0] = 1.0; - im_cg.start(ucell, force, etot); + std::ofstream ofs1("TestStartTrial_temp1.log"); + im_cg.start(ucell, force, etot, ofs1); + ofs1.close(); + std::remove("TestStartTrial_temp1.log"); Ions_Move_Basic::istep = 2; im_cg.move0[0] = 10.0; - GlobalV::ofs_running.open("TestStartTrial.log"); - im_cg.start(ucell, force, etot); - GlobalV::ofs_running.close(); + std::ofstream ofs("TestStartTrial.log"); + im_cg.start(ucell, force, etot, ofs); + ofs.close(); // Check output std::string expected_output = "\n Largest force is 0.257111 eV/Angstrom while threshold is -1 eV/Angstrom\n\n" @@ -249,14 +255,20 @@ TEST_F(IonsMoveCGTest, TestStartNoTrialGotoCase1) // call function im_cg.move0[0] = 1.0; - im_cg.start(ucell, force, etot); + std::ofstream ofs1("TestStartNoTrialGotoCase1_temp1.log"); + im_cg.start(ucell, force, etot, ofs1); + ofs1.close(); + std::remove("TestStartNoTrialGotoCase1_temp1.log"); Ions_Move_Basic::istep = 2; - im_cg.start(ucell, force, etot); + std::ofstream ofs2("TestStartNoTrialGotoCase1_temp2.log"); + im_cg.start(ucell, force, etot, ofs2); + ofs2.close(); + std::remove("TestStartNoTrialGotoCase1_temp2.log"); im_cg.move0[0] = 1.0; force(0, 0) = 0.001; - GlobalV::ofs_running.open("TestStartNoTrialGotoCase1.log"); - im_cg.start(ucell, force, etot); - GlobalV::ofs_running.close(); + std::ofstream ofs("TestStartNoTrialGotoCase1.log"); + im_cg.start(ucell, force, etot, ofs); + ofs.close(); // Check output std::string expected_output = "\n Largest force is 0.0257111 eV/Angstrom while threshold is -1 eV/Angstrom\n\n" @@ -289,13 +301,19 @@ TEST_F(IonsMoveCGTest, TestStartNoTrialGotoCase2) // call function im_cg.move0[0] = 1.0; - im_cg.start(ucell, force, etot); + std::ofstream ofs1("TestStartNoTrialGotoCase2_temp1.log"); + im_cg.start(ucell, force, etot, ofs1); + ofs1.close(); + std::remove("TestStartNoTrialGotoCase2_temp1.log"); Ions_Move_Basic::istep = 2; im_cg.move0[0] = 10.0; - im_cg.start(ucell, force, etot); - GlobalV::ofs_running.open("TestStartNoTrialGotoCase2.log"); - im_cg.start(ucell, force, etot); - GlobalV::ofs_running.close(); + std::ofstream ofs2("TestStartNoTrialGotoCase2_temp2.log"); + im_cg.start(ucell, force, etot, ofs2); + ofs2.close(); + std::remove("TestStartNoTrialGotoCase2_temp2.log"); + std::ofstream ofs("TestStartNoTrialGotoCase2.log"); + im_cg.start(ucell, force, etot, ofs); + ofs.close(); // Check output std::string expected_output = "\n Largest force is 0.257111 eV/Angstrom while threshold is -1 eV/Angstrom\n\n" @@ -328,14 +346,20 @@ TEST_F(IonsMoveCGTest, TestStartNoTrial) // call function im_cg.move0[0] = 1.0; - im_cg.start(ucell, force, etot); + std::ofstream ofs1("TestStartNoTrial_temp1.log"); + im_cg.start(ucell, force, etot, ofs1); + ofs1.close(); + std::remove("TestStartNoTrial_temp1.log"); Ions_Move_Basic::istep = 2; im_cg.move0[0] = 1.0; force(0, 0) = 0.001; - im_cg.start(ucell, force, etot); - GlobalV::ofs_running.open("TestStartNoTrial.log"); - im_cg.start(ucell, force, etot); - GlobalV::ofs_running.close(); + std::ofstream ofs2("TestStartNoTrial_temp2.log"); + im_cg.start(ucell, force, etot, ofs2); + ofs2.close(); + std::remove("TestStartNoTrial_temp2.log"); + std::ofstream ofs("TestStartNoTrial.log"); + im_cg.start(ucell, force, etot, ofs); + ofs.close(); // Check output std::string expected_output = "\n Largest force is 0.0257111 eV/Angstrom while threshold is -1 eV/Angstrom\n\n" diff --git a/source/source_relax/test/ions_move_sd_test.cpp b/source/source_relax/test/ions_move_sd_test.cpp index 647247ecd38..328cb831c40 100644 --- a/source/source_relax/test/ions_move_sd_test.cpp +++ b/source/source_relax/test/ions_move_sd_test.cpp @@ -74,19 +74,19 @@ TEST_F(IonsMoveSDTest, TestStartConverged) double etot = 0.0; // call function - GlobalV::ofs_running.open("log"); - im_sd.start(ucell, force, etot); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_sd_start_converged.log"); + im_sd.start(ucell, force, etot, ofs); + ofs.close(); // Check output std::string expected_output = "\n Largest force is 0 eV/Angstrom while threshold is -1 eV/Angstrom\n" " largest force is 0, no movement is possible.\n it may converged, otherwise no " "movement of atom is allowed.\n end of geometry optimization\n " " istep = 1\n update iteration = 5\n"; - std::ifstream ifs("log"); + std::ifstream ifs("test_sd_start_converged.log"); std::string output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_sd_start_converged.log"); std::regex pattern(R"(==> .*::.*\t[\d\.]+ GB\t\d+ s\n )"); output = std::regex_replace(output, pattern, ""); @@ -127,17 +127,17 @@ TEST_F(IonsMoveSDTest, TestStartNotConverged) } // call function - GlobalV::ofs_running.open("log"); - im_sd.start(ucell, force, etot); - GlobalV::ofs_running.close(); + std::ofstream ofs("test_sd_start_not_converged.log"); + im_sd.start(ucell, force, etot, ofs); + ofs.close(); // Check output std::string expected_output = "\n Largest force is 25.7111 eV/Angstrom while threshold is -1 eV/Angstrom\n\n" " Ion relaxation is not converged yet (threshold is 0.0257111)\n"; - std::ifstream ifs("log"); + std::ifstream ifs("test_sd_start_not_converged.log"); std::string output((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); ifs.close(); - std::remove("log"); + std::remove("test_sd_start_not_converged.log"); EXPECT_THAT(output, testing::HasSubstr(expected_output)); EXPECT_EQ(Ions_Move_Basic::converged, false);