diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index 665c9f9da6..47670d1735 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -1487,15 +1487,20 @@ These variables are used to control the geometry relaxation.
 
 ### relax_method
 
-- **Type**: String
+- **Type**: Vector of string
 - **Description**: The methods to do geometry optimization.
+  the first element:
   - 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 traditional Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm.
+  - bfgs : using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm.
+  - lbfgs: using the Limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) 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.
-- **Default**: cg
+
+  the second element:
+  when the first element is bfgs, if the second parameter is 1, it indicates the use of the new BFGS algorithm; if the second parameter is not 1, it indicates the use of the old BFGS algorithm.
+- **Default**: cg 1
+- **Note**:In the 3.10-LTS version, the type of this parameter is std::string. It can be set to "cg","bfgs","cg_bfgs","bfgs_trad","lbfgs","sd","fire".
 
 ### relax_new
 
diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h
index 15ab2d3a3c..a8dd6d63bb 100644
--- a/source/source_io/module_parameter/input_parameter.h
+++ b/source/source_io/module_parameter/input_parameter.h
@@ -6,6 +6,7 @@
 #include <string>
 #include <vector>
 
+
 // It stores all input parameters both defined in INPUT file and not defined in
 // INPUT file
 struct Input_para
@@ -150,7 +151,7 @@ struct Input_para
     //  int		bessel_nao_lmax;		///< lmax used in descriptor
 
     // ==============   #Parameters (4.Relaxation) ===========================
-    std::string relax_method = "cg"; ///< methods to move_ion: sd, bfgs, cg...
+    std::vector<std::string> relax_method = {"cg","1"}; ///< methods to move_ion: sd, bfgs, cg...
     bool relax_new = true;
     bool relax = false; ///< allow relaxation along the specific direction
     double relax_scale_force = 0.5;
diff --git a/source/source_io/read_input_item_relax.cpp b/source/source_io/read_input_item_relax.cpp
index 036252f23b..59c6ab187a 100644
--- a/source/source_io/read_input_item_relax.cpp
+++ b/source/source_io/read_input_item_relax.cpp
@@ -5,28 +5,53 @@
 
 namespace ModuleIO
 {
+
+
 void ReadInput::item_relax()
 {
     {
         Input_Item item("relax_method");
         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","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");
-                ModuleBase::WARNING_QUIT("ReadInput", warningstr);
-            }
+        item.read_value = [](const Input_Item& item, Parameter& para) {
+        if(item.get_size()==1)
+        {
+            para.input.relax_method[0] = item.str_values[0];
+            para.input.relax_method[1] = "1"; 
+        }
+        else if(item.get_size()>=2)
+        {
+            para.input.relax_method[0] = item.str_values[0];
+            para.input.relax_method[1] = item.str_values[1];
+        }
+        };
+        item.check_value = [](const Input_Item& item, const Parameter& para) {          
+        const std::vector<std::string> relax_methods = {"cg", "sd", "cg_bfgs","lbfgs","bfgs"};
+        if (std::find(relax_methods.begin(), relax_methods.end(), para.input.relax_method[0]) == relax_methods.end()) {
+            const std::string warningstr = nofound_str(relax_methods, "relax_method");
+            ModuleBase::WARNING_QUIT("ReadInput", warningstr);
+        }
         };
         this->add_item(item);
+        
+        // Input_Item item("relax_method");
+        // 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_old", "sd", "cg_bfgs","bfgs","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");
+        //         ModuleBase::WARNING_QUIT("ReadInput", warningstr);
+        //     }
+        // };
+        // this->add_item(item);
     }
     {
         Input_Item item("relax_new");
         item.annotation = "whether to use the new relaxation method";
         read_sync_bool(input.relax_new);
         item.reset_value = [](const Input_Item& item, Parameter& para) {
-            if (para.input.relax_new && para.input.relax_method != "cg")
+            if (para.input.relax_new && para.input.relax_method[0] != "cg")
             {
                 para.input.relax_new = false;
             }
diff --git a/source/source_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp
index 7e0287d90b..a8cdfde867 100644
--- a/source/source_io/test/read_input_ptest.cpp
+++ b/source/source_io/test/read_input_ptest.cpp
@@ -114,7 +114,7 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_EQ(param.inp.fixed_axes, "None");
     EXPECT_FALSE(param.inp.fixed_ibrav);
     EXPECT_FALSE(param.inp.fixed_atoms);
-    EXPECT_EQ(param.inp.relax_method, "cg");
+    EXPECT_EQ(param.inp.relax_method[0], "cg");
     EXPECT_DOUBLE_EQ(param.inp.relax_cg_thr, 0.5);
     EXPECT_EQ(param.inp.out_level, "ie");
     EXPECT_TRUE(param.globalv.out_md_control);
diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp
index 8e1eaba646..fb539164f4 100644
--- a/source/source_io/test_serial/read_input_item_test.cpp
+++ b/source/source_io/test_serial/read_input_item_test.cpp
@@ -766,7 +766,7 @@ TEST_F(InputTest, Item_test)
     }
     { // relax_method
         auto it = find_label("relax_method", readinput.input_lists);
-        param.input.relax_method = "none";
+        param.input.relax_method[0] = "none";
         testing::internal::CaptureStdout();
         EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), "");
         output = testing::internal::GetCapturedStdout();
@@ -775,12 +775,12 @@ TEST_F(InputTest, Item_test)
     { //relax_new
         auto it = find_label("relax_new", readinput.input_lists);
         param.input.relax_new = true;
-        param.input.relax_method = "cg";
+        param.input.relax_method[0] = "cg";
         it->second.reset_value(it->second, param);
         EXPECT_EQ(param.input.relax_new, true);
 
         param.input.relax_new = true;
-        param.input.relax_method = "none";
+        param.input.relax_method[0] = "none";
         it->second.reset_value(it->second, param);
         EXPECT_EQ(param.input.relax_new, false);
     }
diff --git a/source/source_relax/ions_move_basic.cpp b/source/source_relax/ions_move_basic.cpp
index 5b3e1b58c6..1bf1e025a6 100644
--- a/source/source_relax/ions_move_basic.cpp
+++ b/source/source_relax/ions_move_basic.cpp
@@ -23,7 +23,7 @@ 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::string Ions_Move_Basic::relax_method = "bfgs";
+std::vector<std::string> Ions_Move_Basic::relax_method = {"bfgs","2"};
 
 void Ions_Move_Basic::setup_gradient(const UnitCell &ucell, const ModuleBase::matrix &force, double *pos, double *grad)
 {
diff --git a/source/source_relax/ions_move_basic.h b/source/source_relax/ions_move_basic.h
index f6966d2d17..8f03085028 100644
--- a/source/source_relax/ions_move_basic.h
+++ b/source/source_relax/ions_move_basic.h
@@ -21,7 +21,7 @@ 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::string relax_method; // relaxation method,
+extern std::vector<std::string> relax_method; // relaxation method,
 extern int out_stru; // output the structure or not
 // funny way to pass this parameter, but nevertheless
 
diff --git a/source/source_relax/ions_move_cg.cpp b/source/source_relax/ions_move_cg.cpp
index 7da54c3831..9b8aa4b5c0 100644
--- a/source/source_relax/ions_move_cg.cpp
+++ b/source/source_relax/ions_move_cg.cpp
@@ -1,5 +1,4 @@
 #include "ions_move_cg.h"
-
 #include "ions_move_basic.h"
 #include "source_base/global_function.h"
 #include "source_base/global_variable.h"
@@ -77,20 +76,16 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
     // 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;
     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;
-
     
     // some arrays
     double *pos = new double[dim];
@@ -100,7 +95,6 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
     double *cg_grad = new double[dim];
     double best_x = 0.0;
     double fmin = 0.0;
-
     int flag = 0;
 
     ModuleBase::GlobalFunc::ZEROS(pos, dim);
@@ -137,7 +131,6 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
     {
         Ions_Move_Basic::check_converged(ucell, grad);
     }
-
     if (Ions_Move_Basic::converged)
     {
         Ions_Move_Basic::terminate(ucell);
@@ -155,21 +148,7 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
                           flag); // we use the last direction ,the last grad and the grad now to get the direction now
             ncggrad++;
 
-            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;
-                }
-            }
-
+            normalize(cg_gradn, cg_grad, dim);
             setup_move(move0, cg_gradn, steplength); // move the atom position
             Ions_Move_Basic::move_atoms(ucell, move0, pos);
 
@@ -181,16 +160,16 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
 
             f_cal(move0, move0, dim, xb); // xb = trial steplength
             f_cal(move0, grad, dim, fa);  // fa is the projection force in this direction
-
             fmax = fa;
             sd = false;
 
-            if (Ions_Move_Basic::relax_method == "cg_bfgs")
+            if (Ions_Move_Basic::relax_method[0] == "cg_bfgs")
             {
                 if (Ions_Move_Basic::largest_grad * ModuleBase::Ry_to_eV / 0.529177
                     < RELAX_CG_THR) // cg to bfgs  by pengfei 13-8-8
                 {
-                    Ions_Move_Basic::relax_method = "bfgs";
+                    Ions_Move_Basic::relax_method[0] = "bfgs";
+                    Ions_Move_Basic::relax_method[1] = "2";
                 }
                 Ions_Move_Basic::best_xxx = steplength;
             }
@@ -204,7 +183,6 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
                 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;
@@ -214,21 +192,7 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
                     goto CG_begin;
                 }
 
-                double norm = 0.0;
-                for (int i = 0; i < dim; ++i)
-                {
-                    norm += pow(cg_grad0[i], 2);
-                }
-                norm = sqrt(norm);
-
-                if (norm != 0.0)
-                {
-                    for (int i = 0; i < dim; ++i)
-                    {
-                        cg_gradn[i] = cg_grad0[i] / norm;
-                    }
-                }
-
+                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))
@@ -238,7 +202,6 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
 
                 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);
@@ -250,7 +213,6 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
             {
                 double xtemp, ftemp;
                 f_cal(move0, grad, dim, fc);
-
                 fmin = std::abs(fc);
                 nbrent++;
 
@@ -275,24 +237,9 @@ void Ions_Move_CG::start(UnitCell &ucell, const ModuleBase::matrix &force, const
                         goto CG_begin;
                     }
 
-                    double norm = 0.0;
-                    for (int i = 0; i < dim; ++i)
-                    {
-                        norm += pow(cg_grad0[i], 2);
-                    }
-                    norm = sqrt(norm);
-
-                    if (norm != 0.0)
-                    {
-                        for (int i = 0; i < dim; ++i)
-                        {
-                            cg_gradn[i] = cg_grad0[i] / norm;
-                        }
-                    }
-
+                    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;
                 }
             }
@@ -342,7 +289,6 @@ void Ions_Move_CG::setup_cg_grad(double *grad,
             cgp_gp += cg_grad0[i] * grad0[i];
             cgp_g += cg_grad0[i] * grad[i];
         }
-
         assert(g_gp != 0.0);
         const double gamma1 = gg / gp_gp; // FR
         // const double gamma2 = -(gg - g_gp)/(cgp_g - cgp_gp);  //CW
@@ -508,3 +454,22 @@ void Ions_Move_CG::setup_move(double *move, double *cg_gradn, const double &trus
     }
     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;
+}
diff --git a/source/source_relax/ions_move_cg.h b/source/source_relax/ions_move_cg.h
index d143902486..662ef0c14f 100644
--- a/source/source_relax/ions_move_cg.h
+++ b/source/source_relax/ions_move_cg.h
@@ -38,6 +38,7 @@ class Ions_Move_CG
                      const double &fb,
                      const double x,
                      double &best_x);
+    void normalize(double *cg_gradn, const double *cg_grad, int dim);
 };
 
 #endif
diff --git a/source/source_relax/ions_move_methods.cpp b/source/source_relax/ions_move_methods.cpp
index 50ab3286c1..15f27d67d3 100644
--- a/source/source_relax/ions_move_methods.cpp
+++ b/source/source_relax/ions_move_methods.cpp
@@ -16,28 +16,28 @@ void Ions_Move_Methods::allocate(const int &natom)
 {
     Ions_Move_Basic::dim = natom * 3;
 
-    if (Ions_Move_Basic::relax_method == "bfgs")
+    if (Ions_Move_Basic::relax_method[0] == "bfgs"&&Ions_Move_Basic::relax_method[1] != "1")
     {
         this->bfgs.allocate();
     }
-    else if (Ions_Move_Basic::relax_method == "sd")
+    else if (Ions_Move_Basic::relax_method[0] == "sd")
     {
         this->sd.allocate();
     }
-    else if (Ions_Move_Basic::relax_method == "cg")
+    else if (Ions_Move_Basic::relax_method[0] == "cg")
     {
         this->cg.allocate();
     }
-    else if (Ions_Move_Basic::relax_method == "cg_bfgs")
+    else if (Ions_Move_Basic::relax_method[0] == "cg_bfgs")
     {
         this->cg.allocate();
         this->bfgs.allocate(); // added by pengfei  13-8-8
     }
-    else if(Ions_Move_Basic::relax_method == "bfgs_trad")
+    else if(Ions_Move_Basic::relax_method[0] == "bfgs"&&Ions_Move_Basic::relax_method[1] == "1")
     {
         this->bfgs_trad.allocate(natom);       
     }
-    else if(Ions_Move_Basic::relax_method == "lbfgs")
+    else if(Ions_Move_Basic::relax_method[0] == "lbfgs")
     {
         this->lbfgs.allocate(natom);       
     }
@@ -56,34 +56,32 @@ void Ions_Move_Methods::cal_movement(const int &istep,
                                      UnitCell &ucell)
 {
     ModuleBase::TITLE("Ions_Move_Methods", "init");
-
     // Ions_Move_Basic::istep = istep;
     Ions_Move_Basic::istep = force_step;
-
-    if (Ions_Move_Basic::relax_method == "bfgs")
+    if (Ions_Move_Basic::relax_method[0] == "bfgs"&&Ions_Move_Basic::relax_method[1] != "1")
     {
         // move_ions
         // output tau
         // check all symmery
         bfgs.start(ucell, f, etot);
     }
-    else if (Ions_Move_Basic::relax_method == "sd")
+    else if (Ions_Move_Basic::relax_method[0] == "sd")
     {
         sd.start(ucell, f, etot);
     }
-    else if (Ions_Move_Basic::relax_method == "cg")
+    else if (Ions_Move_Basic::relax_method[0] == "cg")
     {
         cg.start(ucell, f, etot);
     }
-    else if (Ions_Move_Basic::relax_method == "cg_bfgs")
+    else if (Ions_Move_Basic::relax_method[0] == "cg_bfgs")
     {
         cg.start(ucell, f, etot); // added by pengfei 13-8-10
     }
-    else if(Ions_Move_Basic::relax_method == "bfgs_trad")
+    else if(Ions_Move_Basic::relax_method[0] == "bfgs"&&Ions_Move_Basic::relax_method[1] == "1")
     {
         bfgs_trad.relax_step(f,ucell);        
     }
-    else if(Ions_Move_Basic::relax_method == "lbfgs")
+    else if(Ions_Move_Basic::relax_method[0] == "lbfgs")
     {
         lbfgs.relax_step(f,ucell,etot);        
     }
diff --git a/source/source_relax/lattice_change_cg.cpp b/source/source_relax/lattice_change_cg.cpp
index 94a1bf581d..442b37e671 100644
--- a/source/source_relax/lattice_change_cg.cpp
+++ b/source/source_relax/lattice_change_cg.cpp
@@ -165,21 +165,7 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_
                           flag); // we use the last direction ,the last grad and the grad now to get the direction now
             ncggrad++;
 
-            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;
-                }
-            }
-
+            normalize(cg_gradn, cg_grad, dim);
             setup_move(move0, cg_gradn, steplength); // move the atom position
             Lattice_Change_Basic::change_lattice(ucell, move0, lat);
 
@@ -214,21 +200,7 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_
                     goto CG_begin;
                 }
 
-                double norm = 0.0;
-                for (int i = 0; i < dim; ++i)
-                {
-                    norm += pow(cg_grad0[i], 2);
-                }
-                norm = sqrt(norm);
-
-                if (norm != 0.0)
-                {
-                    for (int i = 0; i < dim; ++i)
-                    {
-                        cg_gradn[i] = cg_grad0[i] / norm;
-                    }
-                }
-
+                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))
@@ -279,21 +251,7 @@ void Lattice_Change_CG::start(UnitCell &ucell, const ModuleBase::matrix &stress_
                         goto CG_begin;
                     }
 
-                    double norm = 0.0;
-                    for (int i = 0; i < dim; ++i)
-                    {
-                        norm += pow(cg_grad0[i], 2);
-                    }
-                    norm = sqrt(norm);
-
-                    if (norm != 0.0)
-                    {
-                        for (int i = 0; i < dim; ++i)
-                        {
-                            cg_gradn[i] = cg_grad0[i] / norm;
-                        }
-                    }
-
+                    normalize(cg_gradn, cg_grad0, dim);
                     setup_move(move, cg_gradn, best_x);
                     Lattice_Change_Basic::change_lattice(ucell, move, lat);
 
@@ -512,3 +470,22 @@ void Lattice_Change_CG::setup_move(double *move, double *cg_gradn, const double
     }
     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;
+}
diff --git a/source/source_relax/lattice_change_cg.h b/source/source_relax/lattice_change_cg.h
index 3a6e724637..b93f22c81b 100644
--- a/source/source_relax/lattice_change_cg.h
+++ b/source/source_relax/lattice_change_cg.h
@@ -40,6 +40,8 @@ class Lattice_Change_CG
                      const double &fb,
                      const double x,
                      double &best_x);
+                  
+    void normalize(double *cg_gradn, const double *cg_grad, int dim);
 };
 
 #endif
diff --git a/source/source_relax/test/ions_move_cg_test.cpp b/source/source_relax/test/ions_move_cg_test.cpp
index f3ce632f0f..c0aa34dcec 100644
--- a/source/source_relax/test/ions_move_cg_test.cpp
+++ b/source/source_relax/test/ions_move_cg_test.cpp
@@ -138,7 +138,7 @@ TEST_F(IonsMoveCGTest, TestStartSd)
     // setup data
     Ions_Move_Basic::istep = 1;
     Ions_Move_Basic::converged = false;
-    Ions_Move_Basic::relax_method = "cg_bfgs";
+    Ions_Move_Basic::relax_method[0] = "cg_bfgs";
     Ions_Move_CG::RELAX_CG_THR = 100.0;
     UnitCell ucell;
     setupucell(ucell);
@@ -162,7 +162,7 @@ TEST_F(IonsMoveCGTest, TestStartSd)
     EXPECT_THAT(output, testing::HasSubstr(expected_output));
     EXPECT_EQ(Ions_Move_Basic::converged, false);
     EXPECT_EQ(Ions_Move_Basic::update_iter, 5);
-    EXPECT_EQ(Ions_Move_Basic::relax_method, "bfgs");
+    EXPECT_EQ(Ions_Move_Basic::relax_method[0], "bfgs");
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::largest_grad, 0.01);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::best_xxx, -1.0);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::relax_bfgs_init, 1.0);
@@ -201,7 +201,7 @@ TEST_F(IonsMoveCGTest, TestStartTrialGoto)
     EXPECT_THAT(output, testing::HasSubstr(expected_output));
     EXPECT_EQ(Ions_Move_Basic::converged, false);
     EXPECT_EQ(Ions_Move_Basic::update_iter, 5);
-    EXPECT_EQ(Ions_Move_Basic::relax_method, "bfgs");
+    EXPECT_EQ(Ions_Move_Basic::relax_method[0], "bfgs");
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::largest_grad, 0.001);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::best_xxx, -1.0);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::relax_bfgs_init, 10.0);
@@ -239,7 +239,7 @@ TEST_F(IonsMoveCGTest, TestStartTrial)
     EXPECT_THAT(output, testing::HasSubstr(expected_output));
     EXPECT_EQ(Ions_Move_Basic::converged, false);
     EXPECT_EQ(Ions_Move_Basic::update_iter, 5);
-    EXPECT_EQ(Ions_Move_Basic::relax_method, "bfgs");
+    EXPECT_EQ(Ions_Move_Basic::relax_method[0], "bfgs");
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::largest_grad, 0.01);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::best_xxx, -1.0);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::relax_bfgs_init, 70.0);
@@ -279,7 +279,7 @@ TEST_F(IonsMoveCGTest, TestStartNoTrialGotoCase1)
     EXPECT_THAT(output, testing::HasSubstr(expected_output));
     EXPECT_EQ(Ions_Move_Basic::converged, false);
     EXPECT_EQ(Ions_Move_Basic::update_iter, 5);
-    EXPECT_EQ(Ions_Move_Basic::relax_method, "bfgs");
+    EXPECT_EQ(Ions_Move_Basic::relax_method[0], "bfgs");
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::largest_grad, 0.001);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::best_xxx, -1.0);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::relax_bfgs_init, 490.0);
@@ -318,7 +318,7 @@ TEST_F(IonsMoveCGTest, TestStartNoTrialGotoCase2)
     EXPECT_THAT(output, testing::HasSubstr(expected_output));
     EXPECT_EQ(Ions_Move_Basic::converged, false);
     EXPECT_EQ(Ions_Move_Basic::update_iter, 5);
-    EXPECT_EQ(Ions_Move_Basic::relax_method, "bfgs");
+    EXPECT_EQ(Ions_Move_Basic::relax_method[0], "bfgs");
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::largest_grad, 0.01);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::best_xxx, -1.0);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::relax_bfgs_init, 70.0);
@@ -358,7 +358,7 @@ TEST_F(IonsMoveCGTest, TestStartNoTrial)
     EXPECT_THAT(output, testing::HasSubstr(expected_output));
     EXPECT_EQ(Ions_Move_Basic::converged, false);
     EXPECT_EQ(Ions_Move_Basic::update_iter, 5);
-    EXPECT_EQ(Ions_Move_Basic::relax_method, "bfgs");
+    EXPECT_EQ(Ions_Move_Basic::relax_method[0], "bfgs");
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::largest_grad, 0.001);
     EXPECT_DOUBLE_EQ(Ions_Move_Basic::best_xxx, -1.0);
     EXPECT_NEAR(Ions_Move_Basic::relax_bfgs_init, 1.2345679012345678, 1e-12);
diff --git a/source/source_relax/test/ions_move_methods_test.cpp b/source/source_relax/test/ions_move_methods_test.cpp
index 27e7315d68..f0027be236 100644
--- a/source/source_relax/test/ions_move_methods_test.cpp
+++ b/source/source_relax/test/ions_move_methods_test.cpp
@@ -40,19 +40,19 @@ class IonsMoveMethodsTest : public ::testing::Test
 // Test the allocate() function
 TEST_F(IonsMoveMethodsTest, Allocate)
 {
-    Ions_Move_Basic::relax_method = "bfgs";
+    Ions_Move_Basic::relax_method[0] = "bfgs";
     imm.allocate(natom);
     EXPECT_EQ(Ions_Move_Basic::dim, 6);
 
-    Ions_Move_Basic::relax_method = "sd";
+    Ions_Move_Basic::relax_method[0] = "sd";
     imm.allocate(natom);
     EXPECT_EQ(Ions_Move_Basic::dim, 6);
 
-    Ions_Move_Basic::relax_method = "cg";
+    Ions_Move_Basic::relax_method[0] = "cg";
     imm.allocate(natom);
     EXPECT_EQ(Ions_Move_Basic::dim, 6);
 
-    Ions_Move_Basic::relax_method = "cg_bfgs";
+    Ions_Move_Basic::relax_method[0] = "cg_bfgs";
     imm.allocate(natom);
     EXPECT_EQ(Ions_Move_Basic::dim, 6);
 }
@@ -60,7 +60,7 @@ TEST_F(IonsMoveMethodsTest, Allocate)
 // Test the allocate() function warning quit
 TEST_F(IonsMoveMethodsTest, AllocateWarningQuit)
 {
-    Ions_Move_Basic::relax_method = "none";
+    Ions_Move_Basic::relax_method[0] = "none";
     GlobalV::ofs_warning.open("log");
     imm.allocate(natom);
     GlobalV::ofs_warning.close();
@@ -81,22 +81,22 @@ TEST_F(IonsMoveMethodsTest, CalMovement)
     const double etot = 0.0;
     UnitCell ucell;
 
-    Ions_Move_Basic::relax_method = "bfgs";
+    Ions_Move_Basic::relax_method[0] = "bfgs";
     imm.allocate(natom);
     imm.cal_movement(istep, force_step, f, etot, ucell);
     EXPECT_EQ(Ions_Move_Basic::istep, force_step);
 
-    Ions_Move_Basic::relax_method = "sd";
+    Ions_Move_Basic::relax_method[0] = "sd";
     imm.allocate(natom);
     imm.cal_movement(istep, force_step, f, etot, ucell);
     EXPECT_EQ(Ions_Move_Basic::istep, force_step);
 
-    Ions_Move_Basic::relax_method = "cg";
+    Ions_Move_Basic::relax_method[0] = "cg";
     imm.allocate(natom);
     imm.cal_movement(istep, force_step, f, etot, ucell);
     EXPECT_EQ(Ions_Move_Basic::istep, force_step);
 
-    Ions_Move_Basic::relax_method = "cg_bfgs";
+    Ions_Move_Basic::relax_method[0] = "cg_bfgs";
     imm.allocate(natom);
     imm.cal_movement(istep, force_step, f, etot, ucell);
     EXPECT_EQ(Ions_Move_Basic::istep, force_step);
@@ -110,7 +110,7 @@ TEST_F(IonsMoveMethodsTest, CalMovementWarningQuit)
     const ModuleBase::matrix f(3, 3);
     const double etot = 0.0;
     UnitCell ucell;
-    Ions_Move_Basic::relax_method = "none";
+    Ions_Move_Basic::relax_method[0] = "none";
     imm.allocate(natom);
 
     GlobalV::ofs_warning.open("log");
diff --git a/tests/01_PW/058_PW_RE_MB/INPUT b/tests/01_PW/058_PW_RE_MB/INPUT
index a501fb7922..cb65f32dae 100644
--- a/tests/01_PW/058_PW_RE_MB/INPUT
+++ b/tests/01_PW/058_PW_RE_MB/INPUT
@@ -8,7 +8,7 @@ calculation      relax
 relax_nmax       2
 cal_force        1
 force_thr_ev     0.01
-relax_method     bfgs
+relax_method     bfgs 2
 relax_new        0
 
 # Self-Consistent Field
diff --git a/tests/01_PW/059_PW_RE_MB_traj/INPUT b/tests/01_PW/059_PW_RE_MB_traj/INPUT
index 2d3058510e..59254d64f2 100644
--- a/tests/01_PW/059_PW_RE_MB_traj/INPUT
+++ b/tests/01_PW/059_PW_RE_MB_traj/INPUT
@@ -16,7 +16,7 @@ scf_nmax  100
 relax_nmax      2
 cal_force       1
 force_thr_ev    0.01
-relax_method    bfgs_trad
+relax_method    bfgs
 
 #Parameters (4.Basis)
 basis_type  pw
diff --git a/tests/02_NAO_Gamma/010_NO_GO_RE_MB/INPUT b/tests/02_NAO_Gamma/010_NO_GO_RE_MB/INPUT
index 71cb9f7036..f0d05d15e8 100644
--- a/tests/02_NAO_Gamma/010_NO_GO_RE_MB/INPUT
+++ b/tests/02_NAO_Gamma/010_NO_GO_RE_MB/INPUT
@@ -17,7 +17,7 @@ scf_nmax			100
 relax_nmax           2
 cal_force           1
 force_thr_ev    0.01
-relax_method     bfgs
+relax_method     bfgs 2
 
 #Parameters (4.Basis)
 basis_type		lcao
diff --git a/tests/03_NAO_multik/26_NO_KP_RE_MB/INPUT b/tests/03_NAO_multik/26_NO_KP_RE_MB/INPUT
index 8478291249..fd216bc63c 100644
--- a/tests/03_NAO_multik/26_NO_KP_RE_MB/INPUT
+++ b/tests/03_NAO_multik/26_NO_KP_RE_MB/INPUT
@@ -17,7 +17,7 @@ scf_nmax			100
 relax_nmax           2
 cal_force           1
 force_thr_ev    0.01
-relax_method     bfgs
+relax_method     bfgs 2
 
 #Parameters (4.Basis)
 basis_type		lcao
diff --git a/tests/09_DeePKS/102_NO_GO_deepks_relax/INPUT b/tests/09_DeePKS/102_NO_GO_deepks_relax/INPUT
index c111a61dc5..d1cce921c3 100644
--- a/tests/09_DeePKS/102_NO_GO_deepks_relax/INPUT
+++ b/tests/09_DeePKS/102_NO_GO_deepks_relax/INPUT
@@ -2,7 +2,7 @@ INPUT_PARAMETERS
 #Parameters (1.General)
 suffix			    autotest
 calculation     	relax
-relax_method        bfgs
+relax_method        bfgs 2
 relax_nmax		    10
 dft_functional		pbe
 
diff --git a/tests/09_DeePKS/102_NO_KP_deepks_relax/INPUT b/tests/09_DeePKS/102_NO_KP_deepks_relax/INPUT
index 329564682c..64337f5fdd 100644
--- a/tests/09_DeePKS/102_NO_KP_deepks_relax/INPUT
+++ b/tests/09_DeePKS/102_NO_KP_deepks_relax/INPUT
@@ -2,7 +2,7 @@ INPUT_PARAMETERS
 #Parameters (1.General)
 suffix			    autotest
 calculation     	relax
-relax_method        bfgs
+relax_method        bfgs 2
 relax_nmax		    10
 dft_functional		pbe