add LBFGS

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<std::string> relax_methods = {"cg", "bfgs", "sd", "cg_bfgs","bfgs_trad"};
+            const std::vector<std::string> 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");

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

source/module_relax/relax_driver.cpp

@@ -12,10 +12,9 @@
 #include "module_cell/print_cell.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)
@@ -27,7 +26,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;

source/module_relax/relax_old/bfgs.cpp

@@ -44,6 +44,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"<<std::endl;
+    for(int i=0;i<size;i++)
+    {
+        for(int j=0;j<3;j++)
+        {
+            std::cout<<force[i][j]<<' ';
+        }
+        std::cout<<std::endl;
+    }*/
 
     int k=0;
     for(int i=0;i<ucell.ntype;i++)
@@ -68,17 +77,25 @@ void BFGS::relax_step(const ModuleBase::matrix& _force,UnitCell& ucell)
 
     this->PrepareStep(force,pos,H,pos0,force0,steplength,dpos,ucell);
     this->DetermineStep(steplength,dpos,maxstep);
-
-    /*std::cout<<"force"<<std::endl;
+    /*std::cout<<"dpos"<<std::endl;
+    for(int i=0;i<force.size();i++)
+    {
+        for(int j=0;j<3;j++)
+        {
+            std::cout<<dpos[i][j]<<' ';
+        }
+        std::cout<<std::endl;
+    }
+    std::cout<<"force"<<std::endl;
     for(int i=0;i<size;i++)
     {
         for(int j=0;j<3;j++)
         {
             std::cout<<force[i][j]<<' ';
         }
         std::cout<<std::endl;
-    }
-    std::cout<<"dpos"<<std::endl;
+    }*/
+    /*std::cout<<"dpos"<<std::endl;
     for(int i=0;i<size;i++)
     {
         for(int j=0;j<3;j++)
@@ -332,9 +349,7 @@ void BFGS::UpdatePos(UnitCell& ucell)
             a[i*3+j]/=ModuleBase::BOHR_TO_A;
         }
     }
-    std::cout<<std::endl;
-    int k=0;
-    unitcell::update_pos_tau(ucell.lat,a,ucell.ntype,ucell.nat,ucell.atoms);
+    ucell.update_pos_tau(a);
     /*double move_ion[3*size];
     ModuleBase::zeros(move_ion, size*3);
 

source/module_relax/relax_old/bfgs.h

@@ -3,7 +3,7 @@
 
 /**
  * @file bfgs.h
- * @author your name ([email protected])
+ * @author 19hello ([email protected])
  * @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<double> steplength;
-    std::vector<std::vector<double>> H;
-    std::vector<double> force0;
+    std::vector<double> steplength;//the length of atoms displacement 
+    std::vector<std::vector<double>> H;//Hessian matrix
+    std::vector<double> force0;//force in previous step
     std::vector<std::vector<double>> force;
-    std::vector<double> pos0;
+    std::vector<double> pos0;//atom pos in previous step(cartesian coordinates)
     std::vector<std::vector<double>> pos;
-    std::vector<double> pos_taud0;
+    std::vector<double> pos_taud0;//atom pos in previous step(relative coordinates)
     std::vector<std::vector<double>> pos_taud;
     std::vector<std::vector<double>> 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<std::vector<double>>& force,std::vector<std::vector<double>>& pos,std::vector<std::vector<double>>& H,std::vector<double>& pos0,std::vector<double>& force0,std::vector<double>& steplength,std::vector<std::vector<double>>& dpos,UnitCell& ucell);
-    void IsRestrain(std::vector<std::vector<double>>& dpos);
+    void relax_step(const ModuleBase::matrix& _force,UnitCell& ucell);//a full iteration step
+    void PrepareStep(std::vector<std::vector<double>>& force,std::vector<std::vector<double>>& pos,std::vector<std::vector<double>>& H,std::vector<double>& pos0,std::vector<double>& force0,std::vector<double>& steplength,std::vector<std::vector<double>>& dpos,UnitCell& ucell);//calculate the atomic displacement in one iteration step
+    void IsRestrain(std::vector<std::vector<double>>& 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<std::vector<double>>& pos);
     void GetPostaud(UnitCell& ucell,std::vector<std::vector<double>>& pos_taud);
-    void Update(std::vector<double>& pos, std::vector<double>& force,std::vector<std::vector<double>>& H,UnitCell& ucell);
-    void DetermineStep(std::vector<double>& steplength,std::vector<std::vector<double>>& dpos,double& maxstep);
-    void UpdatePos(UnitCell& ucell);
+    void Update(std::vector<double>& pos, std::vector<double>& force,std::vector<std::vector<double>>& H,UnitCell& ucell);//update hessian matrix
+    void DetermineStep(std::vector<double>& steplength,std::vector<std::vector<double>>& dpos,double& maxstep);//normalize large atomic displacements based on maxstep
+    void UpdatePos(UnitCell& ucell);//update ucell with the new coordinates 
 
 
     // matrix method

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.");

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