From 2a79ecb7fa17117b47a67213f1cb3aa5a593a74f Mon Sep 17 00:00:00 2001 From: feiyang <2100011095@stu.pku.edu.cn> Date: Mon, 2 Dec 2024 07:39:10 +0000 Subject: [PATCH 1/7] update bfgs method --- docs/advanced/input_files/input-main.md | 1 + docs/advanced/opt.md | 4 ++++ source/module_relax/relax_old/bfgs.cpp | 6 +++--- source/module_relax/relax_old/bfgs.h | 4 ++-- 4 files changed, 10 insertions(+), 5 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index ced6c56c54..81d077f2db 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -1369,6 +1369,7 @@ These variables are used to control the geometry relaxation. - **Description**: The methods to do geometry optimization. - cg: using the conjugate gradient (CG) algorithm. Note that there are two implementations of the conjugate gradient (CG) method, see [relax_new](#relax_new). - bfgs: using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. + - bfgs_trad: using BFGS algorithm consistent with the BFGS algorithm in ASE. - cg_bfgs: using the CG method for the initial steps, and switching to BFGS method when the force convergence is smaller than [relax_cg_thr](#relax_cg_thr). - sd: using the steepest descent (SD) algorithm. - fire: the Fast Inertial Relaxation Engine method (FIRE), a kind of molecular-dynamics-based relaxation algorithm, is implemented in the molecular dynamics (MD) module. The algorithm can be used by setting [calculation](#calculation) to `md` and [md_type](#md_type) to `fire`. Also ionic velocities should be set in this case. See [fire](../md.md#fire) for more details. diff --git a/docs/advanced/opt.md b/docs/advanced/opt.md index e27800aa4d..d4f687f275 100644 --- a/docs/advanced/opt.md +++ b/docs/advanced/opt.md @@ -24,6 +24,10 @@ The [BFGS method](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%9 In ABACUS, we implemented the BFGS method for doing fixed-cell structural relaxation. +### BFGS_TRAD method + +The BFGS_TRAD method is an algorithm implemented in ABACUS, referencing the BFGS method from ASE. The previous BFGS method in ABACUS did not perform well for some tests and the BFGS_TRAD method now produces results that are consistant with the BFGS method in ASE. In cases where the previous BFGS method could not converge within a limited number of steps, the BFGS_TRAD method can successfully converge. + ### SD method The [SD (steepest descent) method](https://en.wikipedia.org/wiki/Gradient_descent) is one of the simplest first-order optimization methods, where in each step the motion is along the direction of the gradient, where the function descents the fastest. diff --git a/source/module_relax/relax_old/bfgs.cpp b/source/module_relax/relax_old/bfgs.cpp index 5b8dbc03d5..2915dfa84b 100644 --- a/source/module_relax/relax_old/bfgs.cpp +++ b/source/module_relax/relax_old/bfgs.cpp @@ -6,7 +6,7 @@ void BFGS::allocate(const int _size) // initialize H0、H、pos0、force0、force { - alpha=70; + alpha=70;//default value in ase is 70 maxstep=PARAM.inp.relax_bfgs_rmax; size=_size; sign =true; @@ -24,7 +24,7 @@ void BFGS::allocate(const int _size) // initialize H0、H、pos0、force0、forc steplength = std::vector(size, 0.0); } -void BFGS::relax_step(ModuleBase::matrix _force,UnitCell& ucell) +void BFGS::relax_step(const ModuleBase::matrix _force,UnitCell& ucell) { GetPos(ucell,pos); GetPostaud(ucell,pos_taud); @@ -359,7 +359,7 @@ void BFGS::IsRestrain(std::vector>& dpos) Ions_Move_Basic::converged = Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177 grad= std::vector(3*size, 0.0); int iat = 0; diff --git a/source/module_relax/relax_old/bfgs.h b/source/module_relax/relax_old/bfgs.h index 5b7ae415c0..44d806afa5 100644 --- a/source/module_relax/relax_old/bfgs.h +++ b/source/module_relax/relax_old/bfgs.h @@ -49,14 +49,14 @@ class BFGS * @param _size */ void allocate(const int _size);//initialize parameters - void relax_step(ModuleBase::matrix _force,UnitCell& ucell);// + void relax_step(const ModuleBase::matrix _force,UnitCell& ucell);// void PrepareStep(std::vector>& force,std::vector>& pos,std::vector>& H,std::vector& pos0,std::vector& force0,std::vector& steplength,std::vector>& dpos,UnitCell& ucell); void IsRestrain(std::vector>& dpos); private: bool sign; - void CalculateLargestGrad(ModuleBase::matrix& _force,UnitCell& ucell); + void CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell); void GetPos(UnitCell& ucell,std::vector>& pos); void GetPostaud(UnitCell& ucell,std::vector>& pos_taud); void Update(std::vector& pos, std::vector& force,std::vector>& H,UnitCell& ucell); From 0b52719c0ce3c983566c8e67879b162f738d74d3 Mon Sep 17 00:00:00 2001 From: feiyang <2100011095@stu.pku.edu.cn> Date: Mon, 2 Dec 2024 09:17:45 +0000 Subject: [PATCH 2/7] update bfgs method and modify the parameter force in relax_step to be passed by reference --- docs/advanced/input_files/input-main.md | 2 +- docs/advanced/opt.md | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 81d077f2db..4957eea134 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -1369,7 +1369,7 @@ These variables are used to control the geometry relaxation. - **Description**: The methods to do geometry optimization. - cg: using the conjugate gradient (CG) algorithm. Note that there are two implementations of the conjugate gradient (CG) method, see [relax_new](#relax_new). - bfgs: using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. - - bfgs_trad: using BFGS algorithm consistent with the BFGS algorithm in ASE. + - bfgs_trad: using the standard Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. - cg_bfgs: using the CG method for the initial steps, and switching to BFGS method when the force convergence is smaller than [relax_cg_thr](#relax_cg_thr). - sd: using the steepest descent (SD) algorithm. - fire: the Fast Inertial Relaxation Engine method (FIRE), a kind of molecular-dynamics-based relaxation algorithm, is implemented in the molecular dynamics (MD) module. The algorithm can be used by setting [calculation](#calculation) to `md` and [md_type](#md_type) to `fire`. Also ionic velocities should be set in this case. See [fire](../md.md#fire) for more details. diff --git a/docs/advanced/opt.md b/docs/advanced/opt.md index d4f687f275..e6d605deb2 100644 --- a/docs/advanced/opt.md +++ b/docs/advanced/opt.md @@ -24,9 +24,9 @@ The [BFGS method](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%9 In ABACUS, we implemented the BFGS method for doing fixed-cell structural relaxation. -### BFGS_TRAD method +We have alse implemented an alternative BFGS method, which can be called by using the keyword 'bfgs_trad'. -The BFGS_TRAD method is an algorithm implemented in ABACUS, referencing the BFGS method from ASE. The previous BFGS method in ABACUS did not perform well for some tests and the BFGS_TRAD method now produces results that are consistant with the BFGS method in ASE. In cases where the previous BFGS method could not converge within a limited number of steps, the BFGS_TRAD method can successfully converge. +The bfgs_trad method is a quasi-Newton method that substitute an approximate matrix B for the Hessian matrix. The optimization direction is determined by the inverse of B, therefore, only the inverse of B is iteratively updated and no time-consuming operations such as matrix inversion are involved. ### SD method From f79ab761b146e9e08bebd55f2f8ef98f11bd017e Mon Sep 17 00:00:00 2001 From: feiyang <2100011095@stu.pku.edu.cn> Date: Mon, 2 Dec 2024 09:33:13 +0000 Subject: [PATCH 3/7] update bfgs method and modify the parameter force in relax_step to be passed by reference --- docs/advanced/opt.md | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/docs/advanced/opt.md b/docs/advanced/opt.md index e6d605deb2..31cad063bd 100644 --- a/docs/advanced/opt.md +++ b/docs/advanced/opt.md @@ -22,11 +22,9 @@ In the nested procedure mentioned above, we used CG method to perform cell relax The [BFGS method](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm) is a quasi-Newton method for solving nonlinear optimization problem. It belongs to the class of quasi-Newton method where the Hessian matrix is approximated during the optimization process. If the initial point is not far from the extrema, BFGS tends to work better than gradient-based methods. -In ABACUS, we implemented the BFGS method for doing fixed-cell structural relaxation. - -We have alse implemented an alternative BFGS method, which can be called by using the keyword 'bfgs_trad'. +There is an alternative standard BFGS method, which can be called by using the keyword 'bfgs_trad'. The bfgs_trad method is a quasi-Newton method that substitute an approximate matrix B for the Hessian matrix. The optimization direction is determined by the inverse of B, therefore, only the inverse of B is iteratively updated and no time-consuming operations such as matrix inversion are involved. -The bfgs_trad method is a quasi-Newton method that substitute an approximate matrix B for the Hessian matrix. The optimization direction is determined by the inverse of B, therefore, only the inverse of B is iteratively updated and no time-consuming operations such as matrix inversion are involved. +In ABACUS, we implemented the BFGS method for doing fixed-cell structural relaxation. ### SD method From 808feb5e793517c0a850bfbbb488b8ce953e1c46 Mon Sep 17 00:00:00 2001 From: feiyang <2100011095@stu.pku.edu.cn> Date: Tue, 10 Dec 2024 05:33:54 +0000 Subject: [PATCH 4/7] Introduce the differences between the two BFGS methods --- docs/advanced/input_files/input-main.md | 2 +- docs/advanced/opt.md | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 4957eea134..68b75bb610 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -1369,7 +1369,7 @@ These variables are used to control the geometry relaxation. - **Description**: The methods to do geometry optimization. - cg: using the conjugate gradient (CG) algorithm. Note that there are two implementations of the conjugate gradient (CG) method, see [relax_new](#relax_new). - bfgs: using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. - - bfgs_trad: using the standard Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. + - bfgs_trad: using the traditional Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. - cg_bfgs: using the CG method for the initial steps, and switching to BFGS method when the force convergence is smaller than [relax_cg_thr](#relax_cg_thr). - sd: using the steepest descent (SD) algorithm. - fire: the Fast Inertial Relaxation Engine method (FIRE), a kind of molecular-dynamics-based relaxation algorithm, is implemented in the molecular dynamics (MD) module. The algorithm can be used by setting [calculation](#calculation) to `md` and [md_type](#md_type) to `fire`. Also ionic velocities should be set in this case. See [fire](../md.md#fire) for more details. diff --git a/docs/advanced/opt.md b/docs/advanced/opt.md index 31cad063bd..9c2023b655 100644 --- a/docs/advanced/opt.md +++ b/docs/advanced/opt.md @@ -22,9 +22,9 @@ In the nested procedure mentioned above, we used CG method to perform cell relax The [BFGS method](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm) is a quasi-Newton method for solving nonlinear optimization problem. It belongs to the class of quasi-Newton method where the Hessian matrix is approximated during the optimization process. If the initial point is not far from the extrema, BFGS tends to work better than gradient-based methods. -There is an alternative standard BFGS method, which can be called by using the keyword 'bfgs_trad'. The bfgs_trad method is a quasi-Newton method that substitute an approximate matrix B for the Hessian matrix. The optimization direction is determined by the inverse of B, therefore, only the inverse of B is iteratively updated and no time-consuming operations such as matrix inversion are involved. +There is an alternative traditional BFGS method, which can be called by using the keyword 'bfgs_trad'. The bfgs_trad method is a quasi-Newton method that substitute an approximate matrix B for the Hessian matrix. The main difference between 'bfgs' and 'bfgs_trad' is that 'bfgs' updates the inverse of matrix B while 'bfgs_trad' updates matrix B and obtains the inverse of B by solving the matrix eigenvalues and taking the reciprocal of the eigenvalues. Both methods are mathematically equivalent, but in some cases, 'bfgs_trad' performs better. -In ABACUS, we implemented the BFGS method for doing fixed-cell structural relaxation. +In ABACUS, we implemented the BFGS method for doing fixed-cell structural relaxation. Users can choose which implementation of BFGS to call by adding the 'bfgs_trad' or 'bfgs' parameter. ### SD method From ebdac5899085b6d42391a23d9ccc2b517d5a2cf8 Mon Sep 17 00:00:00 2001 From: feiyang <2100011095@stu.pku.edu.cn> Date: Wed, 11 Dec 2024 13:06:39 +0000 Subject: [PATCH 5/7] Add & in relax_step input parameters --- source/module_relax/relax_old/bfgs.cpp | 2 +- source/module_relax/relax_old/bfgs.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/source/module_relax/relax_old/bfgs.cpp b/source/module_relax/relax_old/bfgs.cpp index 2915dfa84b..44f4bf8bf4 100644 --- a/source/module_relax/relax_old/bfgs.cpp +++ b/source/module_relax/relax_old/bfgs.cpp @@ -24,7 +24,7 @@ void BFGS::allocate(const int _size) // initialize H0、H、pos0、force0、forc steplength = std::vector(size, 0.0); } -void BFGS::relax_step(const ModuleBase::matrix _force,UnitCell& ucell) +void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell) { GetPos(ucell,pos); GetPostaud(ucell,pos_taud); diff --git a/source/module_relax/relax_old/bfgs.h b/source/module_relax/relax_old/bfgs.h index 44d806afa5..fba422c886 100644 --- a/source/module_relax/relax_old/bfgs.h +++ b/source/module_relax/relax_old/bfgs.h @@ -49,7 +49,7 @@ class BFGS * @param _size */ void allocate(const int _size);//initialize parameters - void relax_step(const ModuleBase::matrix _force,UnitCell& ucell);// + void relax_step(const ModuleBase::matrix& _force,UnitCell& ucell);// void PrepareStep(std::vector>& force,std::vector>& pos,std::vector>& H,std::vector& pos0,std::vector& force0,std::vector& steplength,std::vector>& dpos,UnitCell& ucell); void IsRestrain(std::vector>& dpos); From 093b2eb69d0fa0577ffb6272f71b04b41eacdf92 Mon Sep 17 00:00:00 2001 From: feiyang <2100011095@stu.pku.edu.cn> Date: Sat, 15 Feb 2025 06:08:51 +0000 Subject: [PATCH 6/7] add lbfgs method --- source/module_io/read_input_item_relax.cpp | 2 +- source/module_relax/CMakeLists.txt | 1 + source/module_relax/relax_driver.cpp | 7 +- source/module_relax/relax_old/bfgs.h | 28 +- .../relax_old/ions_move_methods.cpp | 12 +- .../relax_old/ions_move_methods.h | 7 +- source/module_relax/relax_old/lbfgs.cpp | 1007 +++++++++++++++++ source/module_relax/relax_old/lbfgs.h | 121 ++ source/module_relax/relax_old/relax_old.cpp | 10 +- source/module_relax/relax_old/relax_old.h | 8 +- .../module_relax/relax_old/test/bfgs_test.cpp | 95 +- 11 files changed, 1229 insertions(+), 69 deletions(-) create mode 100644 source/module_relax/relax_old/lbfgs.cpp create mode 100644 source/module_relax/relax_old/lbfgs.h diff --git a/source/module_io/read_input_item_relax.cpp b/source/module_io/read_input_item_relax.cpp index 2995843e7f..1d3e0c62a1 100644 --- a/source/module_io/read_input_item_relax.cpp +++ b/source/module_io/read_input_item_relax.cpp @@ -12,7 +12,7 @@ void ReadInput::item_relax() item.annotation = "cg; bfgs; sd; cg; cg_bfgs;"; read_sync_string(input.relax_method); item.check_value = [](const Input_Item& item, const Parameter& para) { - const std::vector relax_methods = {"cg", "bfgs", "sd", "cg_bfgs","bfgs_trad"}; + const std::vector relax_methods = {"cg", "bfgs", "sd", "cg_bfgs","bfgs_trad","lbfgs"}; if (std::find(relax_methods.begin(),relax_methods.end(), para.input.relax_method)==relax_methods.end()) { const std::string warningstr = nofound_str(relax_methods, "relax_method"); diff --git a/source/module_relax/CMakeLists.txt b/source/module_relax/CMakeLists.txt index ab2925e700..b66ac719d1 100644 --- a/source/module_relax/CMakeLists.txt +++ b/source/module_relax/CMakeLists.txt @@ -8,6 +8,7 @@ add_library( relax_new/line_search.cpp relax_old/bfgs.cpp + relax_old/lbfgs.cpp relax_old/relax_old.cpp relax_old/bfgs_basic.cpp relax_old/ions_move_basic.cpp diff --git a/source/module_relax/relax_driver.cpp b/source/module_relax/relax_driver.cpp index 752ce8687e..73604f61ab 100644 --- a/source/module_relax/relax_driver.cpp +++ b/source/module_relax/relax_driver.cpp @@ -10,10 +10,9 @@ #include "module_io/write_wfc_r.h" #include "module_parameter/parameter.h" void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& ucell) -{ +{ ModuleBase::TITLE("Ions", "opt_ions"); ModuleBase::timer::tick("Ions", "opt_ions"); - if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax" ) { if (!PARAM.inp.relax_new) @@ -25,7 +24,6 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& uce rl.init_relax(ucell.nat); } } - this->istep = 1; int force_step = 1; // pengfei Li 2018-05-14 int stress_step = 1; @@ -89,7 +87,8 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& uce force, stress, force_step, - stress_step); // pengfei Li 2018-05-14 + stress_step, + p_esolver); // pengfei Li 2018-05-14 } // print structure // changelog 20240509 diff --git a/source/module_relax/relax_old/bfgs.h b/source/module_relax/relax_old/bfgs.h index fba422c886..95b29d5e0a 100644 --- a/source/module_relax/relax_old/bfgs.h +++ b/source/module_relax/relax_old/bfgs.h @@ -3,7 +3,7 @@ /** * @file bfgs.h - * @author your name (you@domain.com) + * @author 19hello (you@domain.com) * @brief * @version 0.1 * @date 2024-11-28 @@ -30,16 +30,16 @@ class BFGS double alpha;//initialize H,diagonal element is alpha double maxstep;//every movement smaller than maxstep - int size;//number of etoms + int size;//number of atoms - std::vector steplength; - std::vector> H; - std::vector force0; + std::vector steplength;//the length of atoms displacement + std::vector> H;//Hessian matrix + std::vector force0;//force in previous step std::vector> force; - std::vector pos0; + std::vector pos0;//atom pos in previous step(cartesian coordinates) std::vector> pos; - std::vector pos_taud0; + std::vector pos_taud0;//atom pos in previous step(relative coordinates) std::vector> pos_taud; std::vector> dpos; @@ -49,19 +49,19 @@ class BFGS * @param _size */ void allocate(const int _size);//initialize parameters - void relax_step(const ModuleBase::matrix& _force,UnitCell& ucell);// - void PrepareStep(std::vector>& force,std::vector>& pos,std::vector>& H,std::vector& pos0,std::vector& force0,std::vector& steplength,std::vector>& dpos,UnitCell& ucell); - void IsRestrain(std::vector>& dpos); + void relax_step(const ModuleBase::matrix& _force,UnitCell& ucell);//a full iteration step + void PrepareStep(std::vector>& force,std::vector>& pos,std::vector>& H,std::vector& pos0,std::vector& force0,std::vector& steplength,std::vector>& dpos,UnitCell& ucell);//calculate the atomic displacement in one iteration step + void IsRestrain(std::vector>& dpos);//check if converged private: - bool sign; + bool sign;//check if this is the first iteration void CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell); void GetPos(UnitCell& ucell,std::vector>& pos); void GetPostaud(UnitCell& ucell,std::vector>& pos_taud); - void Update(std::vector& pos, std::vector& force,std::vector>& H,UnitCell& ucell); - void DetermineStep(std::vector& steplength,std::vector>& dpos,double& maxstep); - void UpdatePos(UnitCell& ucell); + void Update(std::vector& pos, std::vector& force,std::vector>& H,UnitCell& ucell);//update hessian matrix + void DetermineStep(std::vector& steplength,std::vector>& dpos,double& maxstep);//normalize large atomic displacements based on maxstep + void UpdatePos(UnitCell& ucell);//update ucell with the new coordinates // matrix method diff --git a/source/module_relax/relax_old/ions_move_methods.cpp b/source/module_relax/relax_old/ions_move_methods.cpp index 60c71fd554..2565a305d4 100644 --- a/source/module_relax/relax_old/ions_move_methods.cpp +++ b/source/module_relax/relax_old/ions_move_methods.cpp @@ -4,6 +4,7 @@ #include "module_base/global_function.h" #include "module_base/global_variable.h" + Ions_Move_Methods::Ions_Move_Methods() { } @@ -36,6 +37,10 @@ void Ions_Move_Methods::allocate(const int &natom) { this->bfgs_trad.allocate(natom); } + else if(Ions_Move_Basic::relax_method == "lbfgs") + { + this->lbfgs.allocate(natom); + } else { ModuleBase::WARNING("Ions_Move_Methods::init", "the parameter Ions_Move_Basic::relax_method is not correct."); @@ -48,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, + ModuleESolver::ESolver* p_esolver) { ModuleBase::TITLE("Ions_Move_Methods", "init"); @@ -78,6 +84,10 @@ void Ions_Move_Methods::cal_movement(const int &istep, { bfgs_trad.relax_step(f,ucell); } + else if(Ions_Move_Basic::relax_method == "lbfgs") + { + lbfgs.relax_step(f,ucell,etot,p_esolver); + } else { ModuleBase::WARNING("Ions_Move_Methods::init", "the parameter Ions_Move_Basic::relax_method is not correct."); diff --git a/source/module_relax/relax_old/ions_move_methods.h b/source/module_relax/relax_old/ions_move_methods.h index 475d4a4083..3364eae3c8 100644 --- a/source/module_relax/relax_old/ions_move_methods.h +++ b/source/module_relax/relax_old/ions_move_methods.h @@ -6,6 +6,9 @@ #include "ions_move_cg.h" #include "ions_move_sd.h" #include "bfgs.h" +#include "lbfgs.h" +#include "module_esolver/esolver.h" +#include "module_esolver/esolver_ks.h" class Ions_Move_Methods { @@ -19,7 +22,8 @@ class Ions_Move_Methods const int &force_step, const ModuleBase::matrix &f, const double &etot, - UnitCell &ucell); + UnitCell &ucell, + ModuleESolver::ESolver* p_esolver); bool get_converged() const { @@ -47,5 +51,6 @@ class Ions_Move_Methods Ions_Move_CG cg; Ions_Move_SD sd; BFGS bfgs_trad; + LBFGS lbfgs; }; #endif diff --git a/source/module_relax/relax_old/lbfgs.cpp b/source/module_relax/relax_old/lbfgs.cpp new file mode 100644 index 0000000000..10fd930c49 --- /dev/null +++ b/source/module_relax/relax_old/lbfgs.cpp @@ -0,0 +1,1007 @@ +#include "lbfgs.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_base/matrix3.h" +#include "module_parameter/parameter.h" +#include "ions_move_basic.h" + +void LBFGS::allocate(const int _size) // initialize H0、H、pos0、force0、force +{ + alpha=70;//default value in ase is 70 + maxstep=PARAM.inp.relax_bfgs_rmax; + size=_size; + sign =true; + memory=100; + damping=70.0; + iteration=0; + xtol=1e-14; + H = std::vector>(3*size, std::vector(3*size, 0.0)); + H0=1/alpha; + pos = std::vector> (size, std::vector(3, 0.0)); + pos0 = std::vector(3*size, 0.0); + pos_taud = std::vector> (size, std::vector(3, 0.0)); + pos_taud0 = std::vector(3*size, 0.0); + dpos = std::vector>(size, std::vector(3, 0.0)); + force0 = std::vector(3*size, 0.0); + force = std::vector>(size, std::vector(3, 0.0)); + steplength = std::vector(size, 0.0); + isave = std::vector(2, 0); + dsave = std::vector(13, 0.0); + old_stp=0; + task=0; +} + +void LBFGS::relax_step(const ModuleBase::matrix _force,UnitCell& ucell,const double &etot,ModuleESolver::ESolver* p_esolver) +{ + GetPos(ucell,pos); + GetPostaud(ucell,pos_taud); + solver=p_esolver; + ucell.ionic_position_updated = true; + for(int i = 0; i < _force.nr; i++) + { + for(int j=0;j<_force.nc;j++) + { + force[i][j]=_force(i,j)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A; + } + } + int k=0; + for(int i=0;iPrepareStep(force,pos,H,pos0,force0,dpos,ucell,etot); + //std::cout<<"enterstep2"<UpdatePos(ucell); + this->CalculateLargestGrad(_force,ucell); + this->IsRestrain(dpos); +} + +void LBFGS::GetPos(UnitCell& ucell,std::vector>& pos) +{ + int k=0; + for(int i=0;i>& pos_taud) +{ + int k=0; + for(int i=0;i>& force, +std::vector>& pos, +std::vector>& H, +std::vector& pos0, +std::vector& force0, +std::vector>& dpos, +UnitCell& ucell, +const double &etot) +{ + std::cout<<"force"< changedforce = this->ReshapeMToV(force); + std::vector changedpos = this->ReshapeMToV(pos); + this->Update(pos_taud,pos_taud0,changedforce,force0,ucell,iteration,memory,s,y,rho); + std::cout<<'s'< q=DotInVAndFloat(changedforce,-1); + int loopmax=std::min(memory,iteration); + std::vector a(loopmax); + for(int i=loopmax-1;i>=0;i--) + { + a[i]=rho[i]*this->DotInVAndV(s[i],q); + auto temp=this->DotInVAndFloat(y[i],a[i]); + q=this->VSubV(q,temp); + } + std::cout<<'q'< z=this->DotInVAndFloat(q,H0); + for(int i=0;iDotInVAndV(y[i],z); + auto temp=DotInVAndFloat(s[i],a[i]-b); + z=VAddV(z,temp); + } + auto temp0=DotInVAndFloat(z,-1); + dpos=ReshapeVToM(temp0); + auto temp1=DotInVAndFloat(changedforce,-1); + std::vector> g=ReshapeVToM(temp1); + energy=etot; + + //alpha_k=LineSearch(ucell,pos,g,energy); + + + //auto temp2=DotInVAndFloat(temp0,alpha_k); + auto temp2=DotInVAndFloat(temp0,1); + dpos=ReshapeVToM(temp2); + std::cout<<"dpos"<ReshapeMToV(pos); + pos_taud0=this->ReshapeMToV(pos_taud); + force0 = changedforce; +} +void LBFGS::Update(std::vector>& pos_taud, + std::vector& pos_taud0, + std::vector& force, + std::vector& force0, + UnitCell& ucell, + int iteration, + int memory, + std::vector>& s, + std::vector>& y, + std::vector& rho) +{ + if(iteration>0) + { + //std::vector dpos=this->VSubV(pos,pos0); + auto term=this->ReshapeMToV(pos_taud); + std::vector dpos = this->VSubV(term, pos_taud0); + for(int i=0;i<3*size;i++) + { + double shortest_move = dpos[i]; + //dpos[i]/=ModuleBase::BOHR_TO_A; + //dpos[i]/=ucell.lat0; + for (int cell = -1; cell <= 1; ++cell) + { + const double now_move = dpos[i] + cell; + if (std::abs(now_move) < std::abs(shortest_move)) + { + shortest_move = now_move; + } + } + //shortest_move=shortest_move*ModuleBase::BOHR_TO_A*ucell.lat0; + dpos[i]=shortest_move; + } + std::vector> c=ReshapeVToM(dpos); + for(int iat=0; iat move_ion_cart; + move_ion_cart.x = c[iat][0] *ModuleBase::BOHR_TO_A * ucell.lat0; + move_ion_cart.y = c[iat][1] * ModuleBase::BOHR_TO_A * ucell.lat0; + move_ion_cart.z = c[iat][2] * ModuleBase::BOHR_TO_A * ucell.lat0; + + //convert pos + ModuleBase::Vector3 move_ion_dr = move_ion_cart* ucell.latvec; + int it = ucell.iat2it[iat]; + int ia = ucell.iat2ia[iat]; + Atom* atom = &ucell.atoms[it]; + if(atom->mbl[ia].x == 1) + { + dpos[iat * 3] = move_ion_dr.x; + } + if(atom->mbl[ia].y == 1) + { + dpos[iat * 3 + 1] = move_ion_dr.y ; + } + if(atom->mbl[ia].z == 1) + { + dpos[iat * 3 + 2] = move_ion_dr.z ; + } + } + std::vector dforce = this->VSubV(force0, force); + double rho0=1.0/this->DotInVAndV(dpos,dforce); + s.push_back(dpos); + y.push_back(dforce); + rho.push_back(rho0); + } + + if(iteration>memory) + { + s.erase(s.begin()); + y.erase(y.begin()); + rho.erase(rho.begin()); + } +} +double LBFGS::LineSearch(UnitCell& ucell,std::vector>& r,std::vector>& g,double e) +{ + std::vector tempp=ReshapeMToV(dpos); + double p_size=0; + for(int i=0;i tempr=ReshapeMToV(r); + std::vector tempg=ReshapeMToV(g); + double c1=0.23; + double c2=0.46; + stpmax=50; + stpmin=1e-8; + xtrapl=1.1; + xtrapu=4; + no_update=false; + double phi0=e; + double derphi0=DotInVAndV(tempg,tempp); + double gms=maxstep*sqrt(3*size); + double alpha1=1; + bool no_update=false; + double fval=e; + double stp=0; + std::vector gval=tempg; + while(true) + { + stp=Step(alpha1,phi0,derphi0,c1,c2,xtol,isave,dsave); + if(task==1) + { + alpha1=stp; + fval=GetEnergy(ucell,stp); + gval=GetForce(ucell,stp); + phi0=fval; + auto temp=ReshapeMToV(dpos); + derphi0=DotInVAndV(gval,temp); + old_stp=alpha1; + if(no_update) + { + break; + } + } + else + { + break; + } + } + return stp; + +} +double LBFGS::Step(double stp,double f,double g,double c1,double c2,double xtol,std::vector& isave,std::vector& dsave) +{ + if(task==0) + { + if(stpstpmax) + { + std::cout<<"ERROR: STP GT maxstep"<=0) + { + std::cout<<"ERROR: INITIAL G >= 0"<=0) + { + stage=2; + } + if(bracket && (stp<=stmin || stp>=stmax)) + { + std::cout<<"WARNING: ROUNDING ERRORS PREVENT PROGRESS"<ftest||g>=gtest)) + { + std::cout<<"WARNING::STP=minstep"<=0.66*width1) + { + stp=stx+0.5*(sty-stx); + } + width1=width; + width=std::abs(sty-stx); + stmin=std::min(stx,sty); + stmax=std::max(stx,sty); + } + else + { + stmin=stp+xtrapl*(stp-stx); + stmax=stp+xtrapu*(stp-stx); + } + stp=std::max(stp,stpmin); + stp=std::min(stp,stpmax); + if(stx==stp&&stp==stpmax&&stmin>stpmax) + { + no_update=true; + } + if((bracket&&stp=stmax||(bracket&&(stmax-stmin)fx) + { + theta=3*(fx-fp)/(stp-stx)+gx+gp; + s=std::max({std::abs(theta),std::abs(gx),std::abs(gp)}); + gamma=s*sqrt((theta*theta-gx*gp)/(s*s)); + if(stp stx) + { + gamma = -gamma; + } + p = (gamma - gp) + theta; + q = ((gamma - gp) + gamma) + gx; + r = p / q; + stpc = stp + r * (stx - stp); + stpq = stp + (gp / (gp - gx)) * (stx - stp); + if (abs(stpc - stp) > abs(stpq - stp)) + { + stpf = stpc; + } + else + { + stpf = stpq; + } + bracket = true; + } + else if (abs(gp) < abs(gx)) + { + theta = 3.0 * (fx - fp) / (stp - stx) + gx + gp; + s = std::max({std::abs(theta), std::abs(gx), std::abs(gp)}); + gamma = s * sqrt(std::max(0.0, std::pow(theta / s, 2.0) - (gx / s) * (gp / s))); + if (stp > stx) + { + gamma = -gamma; + } + p = (gamma - gp) + theta; + q = (gamma + (gx - gp)) + gamma; + r = p / q; + if (r < 0.0 && gamma != 0.0) + { + stpc = stp + r * (stx - stp); + } + else if (stp > stx) + { + stpc = stpmax; + } + else + { + stpc = stpmin; + } + stpq = stp + (gp / (gp - gx)) * (stx - stp); + if (bracket) + { + if (abs(stpc - stp) < abs(stpq - stp)) + { + stpf = stpc; + } + else + { + stpf = stpq; + } + if (stp > stx) + { + stpf = std::min(stp + 0.66 * (sty - stp), stpf); + } + else + { + stpf = std::max(stp + 0.66 * (sty - stp), stpf); + } + } + else + { + if (abs(stpc - stp) > abs(stpq - stp)) + { + stpf = stpc; + } + else + { + stpf = stpq; + } + stpf = std::min(stpmax, stpf); + stpf = std::max(stpmin, stpf); + } + } + else + { + if(bracket) + { + theta = 3.0 * (fp - fy) / (sty - stp) + gy + gp; + s = std::max({std::abs(theta), std::abs(gy), std::abs(gp)}); + gamma = s * sqrt(pow(theta / s, 2.0) - (gy / s) * (gp / s)); + if (stp > sty) + { + gamma = -gamma; + } + p = (gamma - gp) + theta; + q = ((gamma - gp) + gamma) + gy; + r = p / q; + stpc = stp + r * (sty - stp); + stpf = stpc; + } + else if (stp > stx) + { + stpf = stpmax; + } + else + { + stpf = stpmin; + } + } + if (fp > fx) + { + sty = stp; + fy = fp; + gy = gp; + } + else + { + if (sign < 0) + { + sty = stx; + fy = fx; + gy = gx; + } + stx = stp; + fx = fp; + gx = gp; + } + stp = DetermineStep(stpf); +} + +double LBFGS::DetermineStep(double stp) +{ + double dr=stp-old_stp; + for(int i = 0; i < size; i++) + { + double k = 0; + for(int j = 0; j < 3; j++) + { + k += dpos[i][j] * dpos[i][j]*dr*dr; + } + steplength[i] = sqrt(k); + } + auto maxsteplength = max_element(steplength.begin(), steplength.end()); + double a = *maxsteplength; + if(a >= maxstep) + { + double scale = maxstep / a; + dr*=scale; + } + double k=old_stp+dr; + return k; +} +void LBFGS::Save(int a,double b,double c,double d,double e,double f,double g,double h,double i,double j,double k,double l,double m,double n) +{ + if(bracket) + { + isave[0]=1; + } + else + { + isave[0]=0; + } + isave[1]=a; + dsave[0]=b; + dsave[1]=c; + dsave[2]=d; + dsave[3]=e; + dsave[4]=f; + dsave[5]=g; + dsave[6]=h; + dsave[7]=i; + dsave[8]=j; + dsave[9]=k; + dsave[10]=l; + dsave[11]=m; + dsave[12]=n; +} + + + +double LBFGS::GetEnergy(UnitCell& ucell,double stp) +{ + double a[3*size]; + for(int i=0;ical_energy(); +} +std::vector LBFGS::GetForce(UnitCell& ucell,double stp) +{ + double a[3*size]; + for(int i=0;ical_force(ucell,b); + std::vector c=std::vector(3*size, 0.0); + for(int i = 0; i < b.nr; i++) + { + for(int j=0;j move_ion_cart; + move_ion_cart.x = dpos[iat][0] / ModuleBase::BOHR_TO_A / ucell.lat0; + move_ion_cart.y = dpos[iat][1] / ModuleBase::BOHR_TO_A / ucell.lat0; + move_ion_cart.z = dpos[iat][2] / ModuleBase::BOHR_TO_A / ucell.lat0; + + //convert to Direct coordinate + //note here the old GT is used + + //convert pos + ModuleBase::Vector3 move_ion_dr = move_ion_cart* ucell.GT; + + + int it = ucell.iat2it[iat]; + int ia = ucell.iat2ia[iat]; + Atom* atom = &ucell.atoms[it]; + + if(atom->mbl[ia].x == 1) + { + move_ion[iat * 3] = move_ion_dr.x; + } + if(atom->mbl[ia].y == 1) + { + move_ion[iat * 3 + 1] = move_ion_dr.y ; + } + if(atom->mbl[ia].z == 1) + { + move_ion[iat * 3 + 2] = move_ion_dr.z ; + } + } + ucell.update_pos_taud(move_ion); + pos = this->MAddM(pos, dpos);*/ +} + +void LBFGS::IsRestrain(std::vector>& dpos) +{ + Ions_Move_Basic::converged = Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177 grad= std::vector(3*size, 0.0); + int iat = 0; + for (int it = 0; it < ucell.ntype; it++) + { + Atom *atom = &ucell.atoms[it]; + for (int ia = 0; ia < ucell.atoms[it].na; ia++) + { + for (int ik = 0; ik < 3; ++ik) + { + if (atom->mbl[ia][ik]) + { + grad[3 * iat + ik] = -_force(iat, ik) * ucell.lat0; + } + } + ++iat; + } + } + Ions_Move_Basic::largest_grad = 0.0; + for (int i = 0; i < 3*size; i++) + { + if (Ions_Move_Basic::largest_grad < std::abs(grad[i])) + { + Ions_Move_Basic::largest_grad = std::abs(grad[i]); + } + } + Ions_Move_Basic::largest_grad /= ucell.lat0; + if (PARAM.inp.out_level == "ie") + { + std::cout << " LARGEST GRAD (eV/A) : " << Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.5291772109 + << std::endl; + } + +} + + +// matrix methods +std::vector LBFGS::ReshapeMToV(std::vector>& matrix) +{ + int size = matrix.size(); + std::vector result; + result.reserve(3*size); + for (const auto& row : matrix) { + result.insert(result.end(), row.begin(), row.end()); + } + return result; +} + +std::vector> LBFGS::MAddM(std::vector>& a, std::vector>& b) +{ + std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); + for(int i = 0; i < a.size(); i++) + { + for(int j = 0; j < a[0].size(); j++) + { + result[i][j] = a[i][j] + b[i][j]; + } + } + return result; +} + +std::vector LBFGS::VSubV(std::vector& a, std::vector& b) +{ + std::vector result = std::vector(a.size(), 0.0); + for(int i = 0; i < a.size(); i++) + { + result[i] = a[i] - b[i]; + } + return result; +} + +std::vector LBFGS::VAddV(std::vector& a, std::vector& b) +{ + std::vector result = std::vector(a.size(), 0.0); + for(int i = 0; i < a.size(); i++) + { + result[i] = a[i] + b[i]; + } + return result; +} + +std::vector> LBFGS::ReshapeVToM(std::vector& matrix) +{ + std::vector> result = std::vector>(matrix.size() / 3, std::vector(3)); + for(int i = 0; i < result.size(); i++) + { + for(int j = 0; j < 3; j++) + { + result[i][j] = matrix[i*3 + j]; + } + } + return result; +} + +std::vector LBFGS::DotInMAndV1(std::vector>& matrix, std::vector& vec) +{ + std::vector result(matrix.size(), 0.0); + for(int i = 0; i < result.size(); i++) + { + for(int j = 0; j < vec.size(); j++) + { + result[i] += matrix[i][j] * vec[j]; + } + } + return result; +} +std::vector LBFGS::DotInMAndV2(std::vector>& matrix, std::vector& vec) +{ + std::vector result(matrix.size(), 0.0); + for(int i = 0; i < result.size(); i++) + { + for(int j = 0; j < vec.size(); j++) + { + result[i] += matrix[j][i] * vec[j]; + } + } + return result; +} + +double LBFGS::DotInVAndV(std::vector& vec1, std::vector& vec2) +{ + double result = 0.0; + for(int i = 0; i < vec1.size(); i++) + { + result += vec1[i] * vec2[i]; + } + return result; +} + +std::vector> LBFGS::OuterVAndV(std::vector& a, std::vector& b) +{ + std::vector> result = std::vector>(a.size(), std::vector(b.size(), 0.0)); + for(int i = 0; i < a.size(); i++) + { + for(int j = 0; j < b.size(); j++) + { + result[i][j] = a[i] * b[j]; + } + } + return result; +} + +std::vector> LBFGS::MPlus(std::vector>& a, double b) +{ + std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); + for(int i = 0; i < a.size(); i++) + { + for(int j = 0; j < a[0].size(); j++) + { + result[i][j] = a[i][j] / b; + } + } + return result; +} + +std::vector> LBFGS::MSubM(std::vector>& a, std::vector>& b) +{ + std::vector> result = std::vector>(a.size(), std::vector(a[0].size(), 0.0)); + for(int i = 0; i < a.size(); i++) + { + for(int j = 0; j < a[0].size(); j++) + { + result[i][j] = a[i][j] - b[i][j]; + } + } + return result; +} + +std::vector LBFGS::DotInVAndFloat(std::vector& vec, double b) +{ + std::vector result(vec.size(), 0.0); + for(int i = 0; i < vec.size(); i++) + { + result[i] = vec[i] * b; + } + return result; +} \ No newline at end of file diff --git a/source/module_relax/relax_old/lbfgs.h b/source/module_relax/relax_old/lbfgs.h new file mode 100644 index 0000000000..8167290bd1 --- /dev/null +++ b/source/module_relax/relax_old/lbfgs.h @@ -0,0 +1,121 @@ +#ifndef LBFGS_H +#define LBFGS_H + +/** + * @file bfgs.h + * @author your name (you@domain.com) + * @brief + * @version 0.1 + * @date 2024-11-28 + * + * @copyright Copyright (c) 2024 + * + */ + +#include +#include +#include +#include +#include"module_base/lapack_connector.h" + +#include "module_base/matrix.h" +#include "module_base/matrix3.h" +#include "module_cell/unitcell.h" +#include "module_esolver/esolver.h" +#include "module_esolver/esolver_ks.h" + + +class LBFGS +{ +public: + + double alpha;//initialize H,diagonal element is alpha + double maxstep;//every movement smaller than maxstep + int size;//number of etoms + int memory; + double damping; + double H0; + int iteration; + double energy; + double alpha_k; + double xtol; + double stpmin; + double stpmax; + bool bracket; + double xtrapl; + double xtrapu; + double old_stp; + bool no_update; + int task; + ModuleESolver::ESolver* solver; + + + std::vector> H; + std::vector force0; + std::vector> force; + std::vector pos0; + std::vector> pos; + std::vector pos_taud0; + std::vector> pos_taud; + std::vector> dpos; + std::vector> s; + std::vector> y; + std::vector rho; + std::vector isave; + std::vector dsave; + std::vector steplength; + + /** + * @brief + * + * @param _size + */ + void allocate(const int _size);//initialize parameters + void relax_step(const ModuleBase::matrix _force,UnitCell& ucell,const double &etot,ModuleESolver::ESolver* p_esolver);// + void PrepareStep(std::vector>& force,std::vector>& pos,std::vector>& H,std::vector& pos0,std::vector& force0,std::vector>& dpos,UnitCell& ucell,const double &etot); + void IsRestrain(std::vector>& dpos); + + +private: + bool sign; + + void CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell); + void GetPos(UnitCell& ucell,std::vector>& pos); + void GetPostaud(UnitCell& ucell,std::vector>& pos_taud); + + void Update(std::vector>& pos_taud, + std::vector& pos_taud0, + std::vector& force, + std::vector& force0, + UnitCell& ucell, + int iteration, + int memory, + std::vector>& s, + std::vector>& y, + std::vector& rho); + double LineSearch(UnitCell& ucell,std::vector>& r,std::vector>& g,double e); + double Step(double stp,double f,double g,double c1,double c2,double xtol,std::vector& isave,std::vector& dsave); + void UpdateLineSearch(double& stx,double& fx,double& gx,double& sty,double& fy,double& gy,double& stp,double& fp,double& gp,double& stpmin,double& stpmax); + double DetermineStep(double stp); + void UpdatePos(UnitCell& ucell); + void Save(int a,double b,double c,double d,double e,double f,double g,double h,double i,double j,double k,double l,double m,double n); + double GetEnergy(UnitCell& ucell,double stp); + std::vector GetForce(UnitCell& ucell,double stp); + + + // matrix method + std::vector ReshapeMToV(std::vector>& matrix); + std::vector> MAddM(std::vector>& a, std::vector>& b); + std::vector VSubV(std::vector& a, std::vector& b); + std::vector VAddV(std::vector& a, std::vector& b); + std::vector> ReshapeVToM(std::vector& matrix); + std::vector DotInMAndV1(std::vector>& matrix, std::vector& vec); + std::vector DotInMAndV2(std::vector>& matrix, std::vector& vec); + double DotInVAndV(std::vector& vec1, std::vector& vec2); + std::vector> OuterVAndV(std::vector& a, std::vector& b); + std::vector> MPlus(std::vector>& a, double b); + std::vector> MSubM(std::vector>& a, std::vector>& b); + std::vector DotInVAndFloat(std::vector& vec, double b); +}; + +#endif // BFGS_H diff --git a/source/module_relax/relax_old/relax_old.cpp b/source/module_relax/relax_old/relax_old.cpp index b18219c6a7..a90d71f263 100644 --- a/source/module_relax/relax_old/relax_old.cpp +++ b/source/module_relax/relax_old/relax_old.cpp @@ -28,7 +28,8 @@ bool Relax_old::relax_step(const int& istep, ModuleBase::matrix force, ModuleBase::matrix stress, int& force_step, - int& stress_step) + int& stress_step, + ModuleESolver::ESolver* p_esolver) { ModuleBase::TITLE("Relax_old", "relax_step"); @@ -49,7 +50,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,p_esolver); if (!converged) { ucell.ionic_position_updated = true; @@ -129,10 +130,11 @@ bool Relax_old::do_relax(const int& istep, const ModuleBase::matrix& ionic_force, const double& total_energy, UnitCell& ucell, - int& jstep) + int& jstep, + ModuleESolver::ESolver* p_esolver) { 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,p_esolver); ++jstep; return IMM.get_converged(); } diff --git a/source/module_relax/relax_old/relax_old.h b/source/module_relax/relax_old/relax_old.h index 0f10e0ca0b..e5f4aab942 100644 --- a/source/module_relax/relax_old/relax_old.h +++ b/source/module_relax/relax_old/relax_old.h @@ -4,6 +4,8 @@ #include "ions_move_methods.h" #include "lattice_change_methods.h" #include "module_cell/unitcell.h" +#include "module_esolver/esolver.h" +#include "module_esolver/esolver_ks.h" class Relax_old { @@ -15,7 +17,8 @@ class Relax_old ModuleBase::matrix force, ModuleBase::matrix stress, int& force_step, - int& stress_step); + int& stress_step, + ModuleESolver::ESolver* p_esolver); private: Ions_Move_Methods IMM; @@ -28,7 +31,8 @@ class Relax_old const ModuleBase::matrix& ionic_force, const double& total_energy, UnitCell& ucell, - int& jstep); + int& jstep, + ModuleESolver::ESolver* p_esolver); bool do_cellrelax(const int& istep, const int& stress_step, const ModuleBase::matrix& stress, diff --git a/source/module_relax/relax_old/test/bfgs_test.cpp b/source/module_relax/relax_old/test/bfgs_test.cpp index ff766de494..e16663bbb6 100644 --- a/source/module_relax/relax_old/test/bfgs_test.cpp +++ b/source/module_relax/relax_old/test/bfgs_test.cpp @@ -5,8 +5,59 @@ #include "module_base/matrix.h" #include "module_relax/relax_old/ions_move_basic.h" +class BFGSTest : public ::testing::Test { +protected: + BFGS bfgs; + UnitCell ucell; + std::vector> force; + + void SetUp() override { + int size = 10; + bfgs.allocate(size); + + ucell.ntype = 2; + ucell.lat0 = 1.0; + ucell.atoms = new Atom[ucell.ntype]; + for (int i = 0; i < ucell.ntype; i++) { + ucell.atoms[i].na = 5; + ucell.atoms[i].tau = std::vector>(5); + ucell.atoms[i].taud = std::vector>(5); + ucell.atoms[i].mbl = std::vector>(5, {1, 1, 1}); + } + + force = std::vector>(size, std::vector(3, 0.0)); + for (int i = 0; i < force.size(); ++i) { + for (int j = 0; j < 3; ++j) { + force[i][j] = -0.1 * (i + 1); + } + } + } +}; + +TEST_F(BFGSTest, PrepareStep) { + bfgs.PrepareStep(force, bfgs.pos, bfgs.H, bfgs.pos0, bfgs.force0, bfgs.steplength, bfgs.dpos, ucell); + EXPECT_EQ(bfgs.steplength.size(), 10); + for (int i = 0; i < 10; ++i) { + EXPECT_GT(bfgs.steplength[i], 0); + } +} + +TEST_F(BFGSTest, RelaxStep) { + ModuleBase::matrix force; + force = ModuleBase::matrix(10, 3); + for (int i = 0; i < force.nr; ++i) + { + for (int j = 0; j < force.nc; ++j) + { + force(i, j) = -0.1 * (i + 1); + } + } + bfgs.relax_step(force, ucell); + EXPECT_TRUE(ucell.ionic_position_updated); + EXPECT_GT(Ions_Move_Basic::largest_grad, 0); +} -TEST(BFGSTest, AllocateTest) { +TEST_F(BFGSTest, AllocateTest) { BFGS bfgs; int size = 5; bfgs.allocate(size); @@ -20,47 +71,7 @@ TEST(BFGSTest, AllocateTest) { } } - -/*TEST(BFGSTest, RelaxStepTest) { - BFGS bfgs; - UnitCell ucell; - ModuleBase::matrix force(3, 3,true); - int size = 3; - - bfgs.allocate(size); - - force(0, 0)=0.1; - force(1, 1)=0.2; - force(2, 2)=0.3; - - ASSERT_NO_THROW(bfgs.relax_step(force, ucell)); - - - EXPECT_EQ(bfgs.pos.size(), size); -} - -TEST(BFGSTest, PrepareStepIntegrationTest) { - BFGS bfgs; - int size = 3; - bfgs.allocate(size); - - std::vector> force = {{1.0, 2.0, 3.0}, {4.0, 5.0, 6.0}, {7.0, 8.0, 9.0}}; - std::vector> pos = {{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}, {2.0, 2.0, 2.0}}; - std::vector> H = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; - std::vector pos0(size, 0.0); - std::vector force0(size, 0.0); - std::vector steplength(size, 0.0); - std::vector> dpos(size, std::vector(size, 0.0)); - - bfgs.PrepareStep(force, pos, H, pos0, force0, steplength, dpos); - - for (double step : steplength) { - EXPECT_GT(step, 0.0); - } -}*/ - - -TEST(BFGSTest, FullStepTest) +TEST_F(BFGSTest, FullStepTest) { BFGS bfgs; UnitCell ucell; From e2ada07517b7931338f9615268ca1b57d24af5e3 Mon Sep 17 00:00:00 2001 From: feiyang <2100011095@stu.pku.edu.cn> Date: Fri, 7 Mar 2025 09:36:35 +0000 Subject: [PATCH 7/7] add lbfgs method --- source/module_relax/relax_old/bfgs.cpp | 27 +++++++++--- source/module_relax/relax_old/lbfgs.cpp | 57 +++++++++++++++++++++---- source/module_relax/relax_old/lbfgs.h | 1 + 3 files changed, 71 insertions(+), 14 deletions(-) diff --git a/source/module_relax/relax_old/bfgs.cpp b/source/module_relax/relax_old/bfgs.cpp index 39bb3c991b..6cf0cb41f2 100644 --- a/source/module_relax/relax_old/bfgs.cpp +++ b/source/module_relax/relax_old/bfgs.cpp @@ -43,6 +43,15 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell) force[i][j]=_force(i,j)*ModuleBase::Ry_to_eV/ModuleBase::BOHR_TO_A; } } + /*std::cout<<"force"<PrepareStep(force,pos,H,pos0,force0,steplength,dpos,ucell); this->DetermineStep(steplength,dpos,maxstep); - - /*std::cout<<"force"<PrepareStep(force,pos,H,pos0,force0,dpos,ucell,etot); + this->DetermineStep(steplength,dpos,maxstep); //std::cout<<"enterstep2"<>& dpos, UnitCell& ucell, const double &etot) { - std::cout<<"force"< changedforce = this->ReshapeMToV(force); std::vector changedpos = this->ReshapeMToV(pos); this->Update(pos_taud,pos_taud0,changedforce,force0,ucell,iteration,memory,s,y,rho); - std::cout<<'s'< q=DotInVAndFloat(changedforce,-1); int loopmax=std::min(memory,iteration); std::vector a(loopmax); @@ -179,12 +189,12 @@ const double &etot) auto temp=this->DotInVAndFloat(y[i],a[i]); q=this->VSubV(q,temp); } - std::cout<<'q'< z=this->DotInVAndFloat(q,H0); for(int i=0;iReshapeMToV(pos); @@ -294,6 +313,24 @@ void LBFGS::Update(std::vector>& pos_taud, rho.erase(rho.begin()); } } +void LBFGS::DetermineStep(std::vector& steplength, + std::vector>& dpos, + double& maxstep) +{ + auto maxsteplength = max_element(steplength.begin(), steplength.end()); + double a = *maxsteplength; + if(a >= maxstep) + { + double scale = maxstep / a; + for(int i = 0; i < size; i++) + { + for(int j=0;j<3;j++) + { + dpos[i][j]*=scale; + } + } + } +} double LBFGS::LineSearch(UnitCell& ucell,std::vector>& r,std::vector>& g,double e) { std::vector tempp=ReshapeMToV(dpos); @@ -778,7 +815,11 @@ void LBFGS::UpdatePos(UnitCell& ucell) a[i*3+j]/=ModuleBase::BOHR_TO_A; } } - int k=0; + /*std::cout<<"a"<& isave,std::vector& dsave); void UpdateLineSearch(double& stx,double& fx,double& gx,double& sty,double& fy,double& gy,double& stp,double& fp,double& gp,double& stpmin,double& stpmax); double DetermineStep(double stp); + void DetermineStep(std::vector& steplength,std::vector>& dpos,double& maxstep); void UpdatePos(UnitCell& ucell); void Save(int a,double b,double c,double d,double e,double f,double g,double h,double i,double j,double k,double l,double m,double n); double GetEnergy(UnitCell& ucell,double stp);