update bfgs_trad method

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 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.

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.
 
+We have alse implemented an alternative 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.
+
 ### 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.

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<double>(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<std::vector<double>>& dpos)
     Ions_Move_Basic::converged = Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177<PARAM.inp.force_thr_ev;
 }
 
-void BFGS::CalculateLargestGrad(ModuleBase::matrix& _force,UnitCell& ucell)
+void BFGS::CalculateLargestGrad(const ModuleBase::matrix& _force,UnitCell& ucell)
 {
     std::vector<double> grad= std::vector<double>(3*size, 0.0);
     int iat = 0;

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<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);
 
 private:
     bool sign;
 
-    void CalculateLargestGrad(ModuleBase::matrix& _force,UnitCell& ucell);
+    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);