update bfgs method and modify the parameter force in relax_step to be…

… passed by reference

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.

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