diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index 9f992294779..8d516560420 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -964,15 +964,21 @@ calculations.
 - **Type**: Real
 - **Description**: In general, the formula of charge mixing can be written as $\rho_{new} = \rho_{old} + \beta * \rho_{update}$, where $\rho_{new}$ represents the new charge density after charge mixing, $\rho_{old}$ represents the charge density in previous step, $\rho_{update}$ is obtained through various mixing methods, and $\beta$ is set by the parameter `mixing_beta`. A lower value of 'mixing_beta' results in less influence of $\rho_{update}$ on $\rho_{new}$, making the self-consistent field (SCF) calculation more stable. However, it may require more steps to achieve convergence.
 We recommend the following options:
-  - **-10.0**: Program will auto set `mixing_beta` and `mixing_gg0` before charge mixing method starts.
-    - Default values of metal system (bandgap <= 1.0 eV) are `mixing_beta=0.2` and `mixing_gg0=1.0`;
-    - Default values of other systems (bandgap > 1.0eV) are `mixing_beta=0.7` and `mixing_gg0=1.0`.
+  - **0.8**: `nspin=1`
+  - **0.4**: `nspin=2`
+  - **0.2**: `nspin=4`
   - **0**: keep charge density unchanged, usually used for restarting with `init_chg=file` or testing.
   - **0.1 or less**: if convergence of SCF calculation is difficult to reach, please try `0 < mixing_beta < 0.1`.
   
-  Note: For low-dimensional large systems, the setup of `mixing_beta=0.1`, `mixing_ndim=20`, and `mixing_gg0=1.5` usually works well.
+  Note: For low-dimensional large systems, the setup of `mixing_beta=0.1`, `mixing_ndim=20`, and `mixing_gg0=1.0` usually works well.
 
-- **Default**: -10.0
+- **Default**: 0.8 for `nspin=1`, 0.4 for `nspin=2`, 0.2 for `nspin=4`.
+
+### mixing_beta_mag
+
+- **Type**: Real
+- **Description**: Mixing parameter of magnetic density.
+- **Default**: `4*mixing_beta`
 
 ### mixing_ndim
 
@@ -985,11 +991,18 @@ We recommend the following options:
 ### mixing_gg0
 
 - **Type**: Real
-- **Description**: Whether to perfom Kerker scaling.
+- **Description**: Whether to perfom Kerker scaling for charge density.
   -  **>0**: The high frequency wave vectors will be suppressed by multiplying a scaling factor $\frac{k^2}{k^2+gg0^2}$. Setting `mixing_gg0 = 1.0` is normally a good starting point. Kerker preconditioner will be automatically turned off if `mixing_beta <= 0.1`.
   -  **0**: No Kerker scaling is performed.
   
   For systems that are difficult to converge, particularly metallic systems, enabling Kerker scaling may aid in achieving convergence.
+- **Default**: 1.0
+
+### mixing_gg0_mag
+
+- **Type**: Real
+- **Description**: Whether to perfom Kerker preconditioner of magnetic density. 
+  Note: we do not recommand to open Kerker preconditioner of magnetic density unless the system is too hard to converge.
 - **Default**: 0.0
 
 ### mixing_tau
diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp
index e86b2958f6f..f8363b7fc15 100644
--- a/source/module_base/global_variable.cpp
+++ b/source/module_base/global_variable.cpp
@@ -248,7 +248,9 @@ std::string of_kernel_file = "WTkernel.txt";
 std::string MIXING_MODE = "broyden";
 double MIXING_BETA = 0.7;
 int MIXING_NDIM = 8;
-double MIXING_GG0 = 0.0;
+double MIXING_GG0 = 1.00;
+double MIXING_BETA_MAG = 1.6;
+double MIXING_GG0_MAG = 1.00;
 bool MIXING_TAU = 0;
 
 //==========================================================
diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h
index d18f927700d..7c69d4cccb7 100644
--- a/source/module_base/global_variable.h
+++ b/source/module_base/global_variable.h
@@ -278,6 +278,8 @@ extern double MIXING_BETA;
 extern int MIXING_NDIM;
 extern double MIXING_GG0;
 extern bool MIXING_TAU;
+extern double MIXING_BETA_MAG;
+extern double MIXING_GG0_MAG;
 
 //==========================================================
 // device flags added by denghui
diff --git a/source/module_base/module_mixing/broyden_mixing.h b/source/module_base/module_mixing/broyden_mixing.h
index 4f3db44f1d8..1b3ee66ff90 100755
--- a/source/module_base/module_mixing/broyden_mixing.h
+++ b/source/module_base/module_mixing/broyden_mixing.h
@@ -6,6 +6,7 @@
 #include "module_base/memory.h"
 #include "module_base/timer.h"
 #include "module_base/tool_title.h"
+#include "module_base/global_variable.h"
 
 namespace Base_Mixing
 {
@@ -36,6 +37,14 @@ class Broyden_Mixing : public Mixing
         this->mixing_beta = mixing_beta;
         this->coef = std::vector<double>(mixing_ndim + 1);
         this->beta = ModuleBase::matrix(mixing_ndim, mixing_ndim, true);
+        if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
+        {
+            this->two_beta = 0;
+        }
+        else if (GlobalV::NSPIN == 2)
+        {
+            this->two_beta = 1;
+        }
     }
     virtual ~Broyden_Mixing() override
     {
@@ -125,21 +134,43 @@ class Broyden_Mixing : public Mixing
         {
             F_tmp[i] = data_out[i] - data_in[i];
         }
-
         // get screened F
         if (screen != nullptr)
             screen(F_tmp.data());
-
         // container::Tensor data = data_in + mixing_beta * F;
         std::vector<FPTYPE> data(length);
+        // mix density and magnetic density sperately
+        if (this->two_beta == 0)
+        {
+            // rho_tot
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
 #endif
-        for (int i = 0; i < length; ++i)
-        {
-            data[i] = data_in[i] + this->mixing_beta * F_tmp[i];
+            for (int i = 0; i < length; ++i)
+            {
+                data[i] = data_in[i] + this->mixing_beta * F_tmp[i];
+            }
         }
+        else if (this->two_beta == 1)
+        {
+            // rho_tot
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
+#endif
+            for (int i = 0; i < length / 2; ++i)
+            {
+                data[i] = data_in[i] + this->mixing_beta * F_tmp[i];
+            }
+            // magnetism
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
+#endif
+            for (int i = length / 2; i < length; ++i)
+            {
+                data[i] = data_in[i] + GlobalV::MIXING_BETA_MAG * F_tmp[i];
+            }
 
+        }
         mdata.push(data.data());
 
         if (!need_calcoef)
@@ -312,6 +343,8 @@ class Broyden_Mixing : public Mixing
     }
     // the number of calculated dF
     int ndim_cal_dF = 0;
+    // if we should considere two beta
+    int two_beta = 0;
 };
 } // namespace Base_Mixing
 #endif
\ No newline at end of file
diff --git a/source/module_base/module_mixing/plain_mixing.h b/source/module_base/module_mixing/plain_mixing.h
index ba18b6f5546..c87e2ee4a2d 100644
--- a/source/module_base/module_mixing/plain_mixing.h
+++ b/source/module_base/module_mixing/plain_mixing.h
@@ -4,6 +4,7 @@
 #include "module_base/memory.h"
 #include "module_base/timer.h"
 #include "module_base/tool_title.h"
+#include "module_base/global_variable.h"
 
 namespace Base_Mixing
 {
@@ -19,6 +20,14 @@ class Plain_Mixing : public Mixing
         this->mixing_beta = mixing_beta;
         this->data_ndim = 1;
         this->coef = std::vector<double>(1, 1.0);
+        if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
+        {
+            this->two_beta = 0;
+        }
+        else if (GlobalV::NSPIN == 2)
+        {
+            this->two_beta = 1;
+        }
     }
     virtual ~Plain_Mixing() override{};
 
@@ -108,16 +117,43 @@ class Plain_Mixing : public Mixing
 
         // container::Tensor data = data_in + mixing_beta * F;
         std::vector<FPTYPE> data(length);
+        if (this->two_beta == 0)
+        {
+            // rho_tot
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
 #endif
-        for (int i = 0; i < length; ++i)
+            for (int i = 0; i < length; ++i)
+            {
+                data[i] = data_in[i] + this->mixing_beta * F_tmp[i];
+            }
+        }
+        else if (this->two_beta == 1)
         {
-            data[i] = data_in[i] + this->mixing_beta * F_tmp[i];
+            // rho_tot
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
+#endif
+            for (int i = 0; i < length / 2; ++i)
+            {
+                data[i] = data_in[i] + this->mixing_beta * F_tmp[i];
+            }
+            // magnetism
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
+#endif
+            for (int i = length / 2; i < length; ++i)
+            {
+                data[i] = data_in[i] + GlobalV::MIXING_BETA_MAG * F_tmp[i];
+            }
+
         }
 
         mdata.push(data.data());
     };
+
+    // if we should considere two beta
+    int two_beta = 0;
 };
 } // namespace Base_Mixing
 #endif
\ No newline at end of file
diff --git a/source/module_base/module_mixing/pulay_mixing.h b/source/module_base/module_mixing/pulay_mixing.h
index b5ca6544abb..04bb0e47130 100644
--- a/source/module_base/module_mixing/pulay_mixing.h
+++ b/source/module_base/module_mixing/pulay_mixing.h
@@ -6,6 +6,7 @@
 #include "module_base/memory.h"
 #include "module_base/timer.h"
 #include "module_base/tool_title.h"
+#include "module_base/global_variable.h"
 
 namespace Base_Mixing
 {
@@ -30,6 +31,14 @@ class Pulay_Mixing : public Mixing
         this->mixing_beta = mixing_beta;
         this->coef = std::vector<double>(mixing_ndim);
         this->beta = ModuleBase::matrix(mixing_ndim, mixing_ndim, true);
+        if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
+        {
+            this->two_beta = 0;
+        }
+        else if (GlobalV::NSPIN == 2)
+        {
+            this->two_beta = 1;
+        }
     }
     virtual ~Pulay_Mixing() override
     {
@@ -124,12 +133,36 @@ class Pulay_Mixing : public Mixing
 
         // container::Tensor data = data_in + mixing_beta * F;
         std::vector<FPTYPE> data(length);
+        if (this->two_beta == 0)
+        {
+            // rho_tot
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
 #endif
-        for (int i = 0; i < length; ++i)
+            for (int i = 0; i < length; ++i)
+            {
+                data[i] = data_in[i] + this->mixing_beta * F_tmp[i];
+            }
+        }
+        else if (this->two_beta == 1)
         {
-            data[i] = data_in[i] + this->mixing_beta * F_tmp[i];
+            // rho_tot
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
+#endif
+            for (int i = 0; i < length / 2; ++i)
+            {
+                data[i] = data_in[i] + this->mixing_beta * F_tmp[i];
+            }
+            // magnetism
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
+#endif
+            for (int i = length / 2; i < length; ++i)
+            {
+                data[i] = data_in[i] + GlobalV::MIXING_BETA_MAG * F_tmp[i];
+            }
+
         }
 
         mdata.push(data.data());
@@ -275,6 +308,8 @@ class Pulay_Mixing : public Mixing
     int mixing_ndim = -1;
     // start index for F
     int start_F = 0;
+    // if we should considere two beta
+    int two_beta = 0;
 };
 } // namespace Base_Mixing
 #endif
\ No newline at end of file
diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp
index 709c49de470..d94a9abb9f8 100644
--- a/source/module_elecstate/module_charge/charge.cpp
+++ b/source/module_elecstate/module_charge/charge.cpp
@@ -258,6 +258,113 @@ void Charge::renormalize_rho(void)
     return;
 }
 
+// for real space magnetic density
+void Charge::get_rho_mag(void)
+{
+    ModuleBase::TITLE("Charge", "get_rho_tot_mag");
+
+    for (int ir = 0; ir < nrxx; ir++)
+    {
+        rho_mag[ir] = rho[0][ir] + rho[1][ir];
+        rho_mag_save[ir] = rho_save[0][ir] + rho_save[1][ir];
+    }
+    for (int ir = 0; ir < nrxx; ir++)
+    {
+        rho_mag[ir + nrxx] = rho[0][ir] - rho[1][ir];
+        rho_mag_save[ir + nrxx] = rho_save[0][ir] - rho_save[1][ir];
+    }
+    return;
+}
+
+void Charge::get_rho_from_mag(void)
+{
+    ModuleBase::TITLE("Charge", "get_rho_from_mag");
+    for (int is = 0; is < nspin; is++)
+    {
+        ModuleBase::GlobalFunc::ZEROS(rho[is], nrxx);
+        //ModuleBase::GlobalFunc::ZEROS(rho_save[is], nrxx);
+    }
+    for (int ir = 0; ir < nrxx; ir++)
+    {
+        rho[0][ir] = 0.5 * (rho_mag[ir] + rho_mag[ir+nrxx]);
+        rho[1][ir] = 0.5 * (rho_mag[ir] - rho_mag[ir+nrxx]);
+    }
+
+    return;
+}
+
+void Charge::allocate_rho_mag(void)
+{
+    rho_mag = new double[nrxx * nspin];
+    rho_mag_save = new double[nrxx * nspin];
+
+    ModuleBase::GlobalFunc::ZEROS(rho_mag, nrxx * nspin);
+    ModuleBase::GlobalFunc::ZEROS(rho_mag_save, nrxx * nspin);
+
+    return;
+}
+
+void Charge::destroy_rho_mag(void)
+{
+    delete[] rho_mag;
+    delete[] rho_mag_save;
+
+    return;
+}
+
+// for reciprocal space magnetic density
+void Charge::get_rhog_mag(void)
+{
+    ModuleBase::TITLE("Charge", "get_rhog_tot_mag");
+
+    for (int ig = 0; ig < ngmc; ig++)
+    {
+        rhog_mag[ig] = rhog[0][ig] + rhog[1][ig];
+        rhog_mag_save[ig] = rhog_save[0][ig] + rhog_save[1][ig];
+    }
+    for (int ig = 0; ig < ngmc; ig++)
+    {
+        rhog_mag[ig + ngmc] = rhog[0][ig] - rhog[1][ig];
+        rhog_mag_save[ig + ngmc] = rhog_save[0][ig] - rhog_save[1][ig];
+    }
+    return;
+}
+
+void Charge::get_rhog_from_mag(void)
+{
+    ModuleBase::TITLE("Charge", "get_rhog_from_mag");
+    for (int is = 0; is < nspin; is++)
+    {
+        ModuleBase::GlobalFunc::ZEROS(rhog[is], ngmc);
+    }
+    for (int ig = 0; ig < ngmc; ig++)
+    {
+        rhog[0][ig] = 0.5 * (rhog_mag[ig] + rhog_mag[ig+ngmc]);
+        rhog[1][ig] = 0.5 * (rhog_mag[ig] - rhog_mag[ig+ngmc]);
+    }
+
+    return;
+}
+
+void Charge::allocate_rhog_mag(void)
+{
+    rhog_mag = new std::complex<double>[ngmc * nspin];
+    rhog_mag_save = new std::complex<double>[ngmc * nspin];
+
+    ModuleBase::GlobalFunc::ZEROS(rhog_mag, ngmc * nspin);
+    ModuleBase::GlobalFunc::ZEROS(rhog_mag_save, ngmc * nspin);
+
+    return;
+}
+
+void Charge::destroy_rhog_mag(void)
+{
+    delete[] rhog_mag;
+    delete[] rhog_mag_save;
+
+    return;
+}
+
 //-------------------------------------------------------
 // superposition of atomic charges contained in the array
 // rho_at (read from pseudopotential files)
diff --git a/source/module_elecstate/module_charge/charge.h b/source/module_elecstate/module_charge/charge.h
index 5abcb64530d..b0ca0b0cdd2 100644
--- a/source/module_elecstate/module_charge/charge.h
+++ b/source/module_elecstate/module_charge/charge.h
@@ -36,6 +36,12 @@ class Charge
     double **rho = nullptr;
     double **rho_save = nullptr;
 
+    // for magnetic density, onyl support nspin=2 now
+    double *rho_mag = nullptr;
+    double *rho_mag_save = nullptr; 
+    std::complex<double> *rhog_mag = nullptr;
+    std::complex<double> *rhog_mag_save = nullptr;
+
     std::complex<double> **rhog = nullptr;
     std::complex<double> **rhog_save = nullptr;
 
@@ -83,6 +89,49 @@ class Charge
 
     void renormalize_rho(void);
 
+    // for magnetic density
+    /**
+     * @brief allocate rho_mag[is*nnrx] and rho_mag_save[is*nnrx]
+    */
+    void allocate_rho_mag(void);
+
+    /**
+     * @brief destroy rho_mag[is*nnrx] and rho_mag_save[is*nnrx]
+    */
+    void destroy_rho_mag(void);
+
+    /**
+     * @brief get rho_mag[is*nnrx]
+    */
+    void get_rho_mag(void);
+
+    /**
+     * @brief get rho[is][nnrx] from rho_mag[is*nnrx]
+    */
+    void get_rho_from_mag(void);
+
+    /**
+     * @brief // allocate rhog_mag[is*ngmc] and rhog_mag_save[is*ngmc]
+    */
+    void allocate_rhog_mag(void);
+
+    /**
+     * @brief destroy rhog_mag[is*ngmc] and rhog_mag_save[is*ngmc]
+    */
+    void destroy_rhog_mag(void);
+
+    /**
+     * @brief get rhog_mag[is*ngmc]
+    */
+    void get_rhog_mag(void);
+
+    /**
+     * @brief get rhog[is][nnrx] from rhog_mag[is*ngmc]
+    */
+    void get_rhog_from_mag(void);
+
+    double sum_rho(void) const;
+
     void save_rho_before_sum_band(void);
 
 	// for non-linear core correction
@@ -131,7 +180,6 @@ class Charge
     int nspin; // number of spins
     ModulePW::PW_Basis* rhopw = nullptr;
   private:
-    double sum_rho(void) const;
 
     void destroy();    // free arrays  liuyu 2023-03-12
 
diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp
index f4efb193366..ff175c85f36 100755
--- a/source/module_elecstate/module_charge/charge_mixing.cpp
+++ b/source/module_elecstate/module_charge/charge_mixing.cpp
@@ -8,6 +8,7 @@
 #include "module_base/parallel_reduce.h"
 #include "module_base/timer.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
+
 Charge_Mixing::Charge_Mixing()
 {
 }
@@ -26,9 +27,21 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in,
     this->mixing_mode = mixing_mode_in;
     this->mixing_beta = mixing_beta_in;
     this->mixing_ndim = mixing_ndim_in;
-    this->mixing_gg0 = mixing_gg0_in; // mohan add 2014-09-27
+    this->mixing_gg0 = mixing_gg0_in;
     this->mixing_tau = mixing_tau_in;
 
+    GlobalV::ofs_running<<"\n----------- Double Check Mixing Parameters Begin ------------"<<std::endl;
+    GlobalV::ofs_running<<"mixing_type: "<< this->mixing_mode <<std::endl;
+    GlobalV::ofs_running<<"mixing_beta: "<< this->mixing_beta <<std::endl;
+    GlobalV::ofs_running<<"mixing_gg0: "<< this->mixing_gg0 <<std::endl;
+    if (GlobalV::NSPIN==2)
+    {
+        GlobalV::ofs_running<<"mixing_beta_mag: "<< GlobalV::MIXING_BETA_MAG <<std::endl;
+        GlobalV::ofs_running<<"mixing_gg0_mag: "<< GlobalV::MIXING_GG0_MAG <<std::endl;
+    }
+    GlobalV::ofs_running<<"mixing_ndim: "<< this->mixing_ndim <<std::endl;
+    GlobalV::ofs_running<<"----------- Double Check Mixing Parameters End ------------"<<std::endl;
+
     if (this->mixing_mode == "broyden")
     {
         delete this->mixing;
@@ -41,6 +54,7 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in,
     }
     else if (this->mixing_mode == "pulay")
     {
+        delete this->mixing;
         this->mixing = new Base_Mixing::Pulay_Mixing(this->mixing_ndim, this->mixing_beta);
     }
     else
@@ -90,36 +104,31 @@ void Charge_Mixing::auto_set(const double& bandgap_in, const UnitCell& ucell_)
         this->autoset = false;
     }
     GlobalV::ofs_running << "--------------AUTO-SET---------------" << std::endl;
-    // 0.2 for metal and 0.7 for others
-    if (bandgap_in * ModuleBase::Ry_to_eV < 1.0)
+    // 0.8 for nspin=1 and 0.4 for others
+    if (GlobalV::NSPIN == 1)
     {
-        this->mixing->mixing_beta = this->mixing_beta = 0.2;
+        this->mixing->mixing_beta = this->mixing_beta = 0.8;
+        GlobalV::ofs_running << "      Autoset mixing_beta to " << this->mixing_beta << std::endl;
     }
     else
     {
-        this->mixing->mixing_beta = this->mixing_beta = 0.7;
+        this->mixing->mixing_beta = this->mixing_beta = 0.4;
+        GlobalV::MIXING_BETA_MAG = 4 * this->mixing_beta;
+        GlobalV::ofs_running << "      Autoset mixing_beta to " << this->mixing_beta << std::endl;
+        GlobalV::ofs_running << "      Autoset mixing_beta_mag to " << GlobalV::MIXING_BETA_MAG << std::endl;
     }
     GlobalV::MIXING_BETA = mixing_beta;
-    GlobalV::ofs_running << "      Autoset mixing_beta to " << this->mixing_beta << std::endl;
 
-    bool has_trans_metal = false;
-    // find elements of cell
-    for (int it = 0; it < ucell_.ntype; it++)
-    {
-        if (ModuleBase::IsTransMetal.find(ucell_.atoms[it].ncpp.psd) != ModuleBase::IsTransMetal.end())
-        {
-            if (ModuleBase::IsTransMetal.at(ucell_.atoms[it].ncpp.psd))
-            {
-                has_trans_metal = true;
-            }
-        }
-    }
     // auto set kerker mixing_gg0 = 1.0 as default
     this->mixing_gg0 = 1.0;
-
     GlobalV::ofs_running << "      Autoset mixing_gg0 to " << this->mixing_gg0 << std::endl;
+    if (GlobalV::NSPIN == 1)
+    {
+        GlobalV::MIXING_GG0_MAG = 0.0;
+        GlobalV::ofs_running << "      Autoset mixing_gg0_mag to " << GlobalV::MIXING_GG0_MAG << std::endl;
+    }
+
     GlobalV::ofs_running << "-------------------------------------" << std::endl;
-    // auto set for inhomogeneous system
 }
 
 double Charge_Mixing::get_drho(Charge* chr, const double nelec)
@@ -197,6 +206,7 @@ void Charge_Mixing::mix_rho_recip(Charge* chr)
     // NOTE: chr->rhopw is dense, while this->rhopw is smooth
     std::complex<double>* rhog_in = nullptr;
     std::complex<double>* rhog_out = nullptr;
+
     if (GlobalV::double_grid)
     {
         rhog_in = new std::complex<double>[GlobalV::NSPIN * this->rhopw->npw];
@@ -323,21 +333,140 @@ void Charge_Mixing::mix_rho_recip(Charge* chr)
     return;
 }
 
-void Charge_Mixing::mix_rho_real(Charge* chr)
+void Charge_Mixing::mix_rho_recip_new(Charge* chr)
 {
-    // electronic density
-    double* rhor_in = chr->rho_save[0];
-    double* rhor_out = chr->rho[0];
+    // not support nspin=4 yet 2023/11/17
+    // old support see mix_rho_recip()
+    if (GlobalV::double_grid)
+    {
+        ModuleBase::WARNING_QUIT("Charge_Mixing", "double_grid is not supported for new mixing method yet.");
+    }
+
+    std::complex<double>* rhog_in = nullptr;
+    std::complex<double>* rhog_out = nullptr;
+    
+    if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
+    {
+        rhog_in = chr->rhog_save[0];
+        rhog_out = chr->rhog[0];    
+    }
+    else if (GlobalV::NSPIN == 2)
+    {
+        chr->allocate_rhog_mag();
+        chr->get_rhog_mag();
+        rhog_in = chr->rhog_mag_save;
+        rhog_out = chr->rhog_mag;
+    }
+
+    auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1);
+    this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, true);
+    //  can choose inner_product_recip_new1 or inner_product_recip_new2
+    auto inner_product_new
+        = std::bind(&Charge_Mixing::inner_product_recip_new1, this, std::placeholders::_1, std::placeholders::_2);
+    auto inner_product_old
+        = std::bind(&Charge_Mixing::inner_product_recip, this, std::placeholders::_1, std::placeholders::_2);
+    if (GlobalV::NSPIN == 2)
+    {
+        this->mixing->cal_coef(this->rho_mdata, inner_product_new);
+    }
+    else if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
+    {
+        this->mixing->cal_coef(this->rho_mdata, inner_product_old);
+    }
+
+    this->mixing->mix_data(this->rho_mdata, rhog_out);
 
+    if (GlobalV::NSPIN == 2)
+    {
+        chr->get_rhog_from_mag();
+        chr->destroy_rhog_mag();
+    }
+
+    // rhog to rho
+    for (int is = 0; is < GlobalV::NSPIN; is++)
+    {
+        chr->rhopw->recip2real(chr->rhog[is], chr->rho[is]);
+    }
+    //chr->renormalize_rho();
+
+    // For kinetic energy density
+    if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau)
+    {
+        std::vector<std::complex<double>> kin_g(GlobalV::NSPIN * chr->rhopw->npw);
+        std::vector<std::complex<double>> kin_g_save(GlobalV::NSPIN * chr->rhopw->npw);
+        std::complex<double>* taug_in = nullptr;
+        std::complex<double>* taug_out = nullptr;
+        if (GlobalV::double_grid)
+        {
+            taug_in = new std::complex<double>[GlobalV::NSPIN * this->rhopw->npw];
+            taug_out = new std::complex<double>[GlobalV::NSPIN * this->rhopw->npw];
+        }
+        else
+        {
+            taug_in = kin_g_save.data();
+            taug_out = kin_g.data();
+        }
+
+        for (int is = 0; is < GlobalV::NSPIN; ++is)
+        {
+            chr->rhopw->real2recip(chr->kin_r[is], &kin_g[is * chr->rhopw->npw]);
+            chr->rhopw->real2recip(chr->kin_r_save[is], &kin_g_save[is * chr->rhopw->npw]);
+
+            // similar to rhog and rhog_save
+            if (GlobalV::double_grid)
+            {
+                std::memcpy(&taug_in[is * rhopw->npw],
+                            &kin_g_save[is * chr->rhopw->npw],
+                            rhopw->npw * sizeof(std::complex<double>));
+                std::memcpy(&taug_out[is * rhopw->npw],
+                            &kin_g[is * chr->rhopw->npw],
+                            rhopw->npw * sizeof(std::complex<double>));
+            }
+        }
+        // Note: there is no kerker modification for tau because I'm not sure
+        // if we should have it. If necessary we can try it in the future.
+        this->mixing->push_data(this->tau_mdata, taug_in, taug_out, nullptr, false);
+
+        this->mixing->mix_data(this->tau_mdata, taug_out);
+
+        // kin_g to kin_r
+        for (int is = 0; is < GlobalV::NSPIN; is++)
+        {
+            chr->rhopw->recip2real(&kin_g[is * chr->rhopw->npw], chr->kin_r[is]);
+        }
+    }
+
+    return;
+}
+
+void Charge_Mixing::mix_rho_real(Charge* chr)
+{
+    double* rhor_in;
+    double* rhor_out;
+    if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
+    {
+        rhor_in = chr->rho_save[0];
+        rhor_out = chr->rho[0];    
+    }
+    else if (GlobalV::NSPIN == 2)
+    {
+        chr->allocate_rho_mag();
+        chr->get_rho_mag();
+        rhor_in = chr->rho_mag_save;
+        rhor_out = chr->rho_mag;
+    }
     auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1);
     this->mixing->push_data(this->rho_mdata, rhor_in, rhor_out, screen, true);
-
     auto inner_product
         = std::bind(&Charge_Mixing::inner_product_real, this, std::placeholders::_1, std::placeholders::_2);
     this->mixing->cal_coef(this->rho_mdata, inner_product);
-
     this->mixing->mix_data(this->rho_mdata, rhor_out);
-
+    if (GlobalV::NSPIN == 2)
+    {
+        chr->get_rho_from_mag();
+        chr->destroy_rho_mag();
+    }
+    chr->renormalize_rho();
     double *taur_out, *taur_in;
     if ((XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) && mixing_tau)
     {
@@ -349,6 +478,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
 
         this->mixing->mix_data(this->tau_mdata, taur_out);
     }
+
 }
 
 void Charge_Mixing::mix_reset()
@@ -427,9 +557,16 @@ void Charge_Mixing::mix_rho(Charge* chr)
     // --------------------Mixing Body--------------------
     if (GlobalV::SCF_THR_TYPE == 1)
     {
-        mix_rho_recip(chr);
+        if (GlobalV::double_grid)
+        {
+            mix_rho_recip(chr);
+        }
+        else // new mixing method do not support double_grid yet
+        {
+            mix_rho_recip_new(chr);
+        }
     }
-    else
+    else if (GlobalV::SCF_THR_TYPE == 2)
     {
         mix_rho_real(chr);
     }
@@ -506,10 +643,53 @@ void Charge_Mixing::Kerker_screen_recip(std::complex<double>* drhog)
     return;
 }
 
-void Charge_Mixing::Kerker_screen_real(double* drhor)
+void Charge_Mixing::Kerker_screen_recip_new(std::complex<double>* drhog)
 {
     if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1)
         return;
+    double fac, gg0, amin;
+
+    // implement Kerker for density and magnetization separately
+    for (int is = 0; is < GlobalV::NSPIN; ++is)
+    {
+        // new mixing method only support nspin=2 not nspin=4
+        if (is == 1 && GlobalV::NSPIN == 2)
+        {
+            if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1)
+            {
+                for (int ig = 0; ig < this->rhopw->npw; ig++)
+                {
+                    drhog[is * this->rhopw->npw + ig] *= 1;
+                }
+                break;
+            }
+            fac = GlobalV::MIXING_GG0_MAG;
+            amin = GlobalV::MIXING_BETA_MAG;
+        }
+        else
+        {
+            fac = this->mixing_gg0;
+            amin = this->mixing_beta;
+        }
+
+        gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 512)
+#endif
+        for (int ig = 0; ig < this->rhopw->npw; ++ig)
+        {
+            double gg = this->rhopw->gg[ig];
+            double filter_g = std::max(gg / (gg + gg0), 0.1 / amin);
+            drhog[is * this->rhopw->npw + ig] *= filter_g;
+        }
+    }
+    return;
+}
+
+void Charge_Mixing::Kerker_screen_real(double* drhor)
+{
+    if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1)
+        return;
     std::vector<std::complex<double>> drhog(this->rhopw->npw * GlobalV::NSPIN);
     std::vector<double> drhor_filter(this->rhopw->nrxx * GlobalV::NSPIN);
     for (int is = 0; is < GlobalV::NSPIN; ++is)
@@ -518,23 +698,48 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
         // Thus we cannot use Kerker_screen_recip(drhog.data()) directly after it.
         this->rhopw->real2recip(drhor + is * this->rhopw->nrxx, drhog.data() + is * this->rhopw->npw);
     }
-    // Kerker_screen_recip(drhog.data()); Note that we can not use it.
-    // However, we should rewrite it:
-    const double fac = this->mixing_gg0;
-    const double gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2);
-#ifdef _OPENMP
-#pragma omp parallel for collapse(2) schedule(static, 512)
-#endif
+    // implement Kerker for density and magnetization separately
+    double fac, gg0, amin;
     for (int is = 0; is < GlobalV::NSPIN; is++)
     {
+
+        if (is == 1 && GlobalV::NSPIN == 2)
+        {
+            if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1)
+            {
+                for (int ig = 0; ig < this->rhopw->npw; ig++)
+                {
+                    drhog[is * this->rhopw->npw + ig] *= 0;
+                }
+                break;
+            }
+            fac = GlobalV::MIXING_GG0_MAG;
+            amin = GlobalV::MIXING_BETA_MAG;
+        }
+        else
+        {
+            fac = this->mixing_gg0;
+            amin = this->mixing_beta;
+        }
+        
+        gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 512)
+#endif
         for (int ig = 0; ig < this->rhopw->npw; ig++)
         {
             double gg = this->rhopw->gg[ig];
-            double filter_g = std::max(gg / (gg + gg0), 0.1 / this->mixing_beta);
+            // I have not decided how to handle gg=0 part, will be changed in future
+            //if (gg == 0)
+            //{
+            //    drhog[is * this->rhopw->npw + ig] *= 0;
+            //    continue;
+            //}
+            double filter_g = std::max(gg / (gg + gg0), 0.1 / amin);
             drhog[is * this->rhopw->npw + ig] *= (1 - filter_g);
         }
     }
-
+    // inverse FT
     for (int is = 0; is < GlobalV::NSPIN; ++is)
     {
         this->rhopw->recip2real(drhog.data() + is * this->rhopw->npw, drhor_filter.data() + is * this->rhopw->nrxx);
@@ -549,6 +754,73 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
     }
 }
 
+// this test will be removed once new mixing method is finished
+void Charge_Mixing::Kerker_screen_real_test(double* drhor)
+{
+    // for total charge density
+    if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1)
+        return;
+    std::vector<double> drhor_filter(this->rhopw->nrxx);
+    std::vector<std::complex<double>> drhog(this->rhopw->npw);
+    // Note after this process some G which is higher than Gmax will be filtered.
+    // FT
+    this->rhopw->real2recip(drhor, drhog.data());
+    const double fac = this->mixing_gg0;
+    const double gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 512)
+#endif
+    for (int ig = 0; ig < this->rhopw->npw; ig++)
+    {
+        double gg = this->rhopw->gg[ig];
+        double filter_g = std::max(gg / (gg + gg0), 0.1 / this->mixing_beta);
+        drhog[ig] *= (1 - filter_g);
+    }
+    // inverse FT
+    this->rhopw->recip2real(drhog.data(), drhor_filter.data());
+
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 512)
+#endif
+    for (int ir = 0; ir < this->rhopw->nrxx; ir++)
+    {
+        drhor[ir] -= drhor_filter[ir];
+    }
+
+    // for magnetic density
+    if (GlobalV::NSPIN == 2)
+    {
+        if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1)
+            return;
+        std::vector<double> drhor_mag_filter(this->rhopw->nrxx);
+        std::vector<std::complex<double>> drhog_mag(this->rhopw->npw);
+        // Note after this process some G which is higher than Gmax will be filtered.
+        // FT
+        this->rhopw->real2recip(drhor + this->rhopw->nrxx, drhog_mag.data());
+        const double fac_mag = GlobalV::MIXING_GG0_MAG;
+        const double gg0_mag = std::pow(fac_mag * 0.529177 / GlobalC::ucell.tpiba, 2);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 512)
+#endif
+        for (int ig = 0; ig < this->rhopw->npw; ig++)
+        {
+            double gg = this->rhopw->gg[ig];
+            double filter_g = std::max(gg / (gg + gg0_mag), 0.1 / GlobalV::MIXING_BETA_MAG);
+            drhog_mag[ig] *= (1 - filter_g);
+        }
+        // inverse FT
+        this->rhopw->recip2real(drhog_mag.data(), drhor_mag_filter.data());
+        // value magnetism
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 512)
+#endif
+        for (int ir = 0; ir < this->rhopw->nrxx; ir++)
+        {
+            drhor[ir + this->rhopw->nrxx] -= drhor_mag_filter[ir];
+        }
+    }
+}
+
 double Charge_Mixing::inner_product_recip(std::complex<double>* rho1, std::complex<double>* rho2)
 {
     std::complex<double>** rho1_2d = new std::complex<double>*[GlobalV::NSPIN];
@@ -564,6 +836,82 @@ double Charge_Mixing::inner_product_recip(std::complex<double>* rho1, std::compl
     return result;
 }
 
+// a simple inner product
+double Charge_Mixing::inner_product_recip_new1(std::complex<double>* rho1, std::complex<double>* rho2)
+{
+    double rnorm = 0.0;
+#ifdef _OPENMP
+#pragma omp parallel for reduction(+ : rnorm)
+#endif
+    for (int ig = 0; ig < this->rhopw->npw * GlobalV::NSPIN; ++ig)
+    {
+        rnorm += (conj(rho1[ig]) * rho2[ig]).real();
+    }
+#ifdef __MPI
+    Parallel_Reduce::reduce_pool(rnorm);
+#endif
+    return rnorm;
+}
+
+// a Hartree-like inner product
+double Charge_Mixing::inner_product_recip_new2(std::complex<double>* rhog1, std::complex<double>* rhog2)
+{
+    static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / GlobalC::ucell.tpiba2;
+    static const double fac2 = ModuleBase::e2 * ModuleBase::FOUR_PI / (ModuleBase::TWO_PI * ModuleBase::TWO_PI);
+
+    double sum = 0.0;
+
+#ifdef _OPENMP
+#pragma omp parallel for reduction(+ : sum)
+#endif
+        for (int ig = 0; ig < this->rhopw->npw; ++ig)
+        {
+            if (this->rhopw->gg[ig] < 1e-8)
+                continue;
+            sum += (conj(rhog1[ig]) * (rhog2[ig])).real() / this->rhopw->gg[ig];
+        }
+        sum *= fac;
+
+        if (GlobalV::GAMMA_ONLY_PW)
+        {
+            sum *= 2.0;
+        }
+
+        // (2) Second part of density error.
+        // including |G|=0 term.
+        double sum2 = 0.0;
+
+        sum2 += fac2 * (conj(rhog1[0 + this->rhopw->npw]) * rhog2[0 + this->rhopw->npw]).real();
+
+        double mag = 0.0;
+#ifdef _OPENMP
+#pragma omp parallel for reduction(+ : mag)
+#endif
+        for (int ig = 0; ig < this->rhopw->npw; ig++)
+        {
+            mag += (conj(rhog1[ig + this->rhopw->npw]) * rhog2[ig + this->rhopw->npw]).real();
+        }
+        mag *= fac2;
+
+        // if(GlobalV::GAMMA_ONLY_PW);
+        if (GlobalV::GAMMA_ONLY_PW) // Peize Lin delete ; 2020.01.31
+        {
+            mag *= 2.0;
+        }
+
+        // std::cout << " sum=" << sum << " mag=" << mag << std::endl;
+        sum2 += mag;
+        sum += sum2;
+
+#ifdef __MPI
+    Parallel_Reduce::reduce_pool(sum);
+#endif
+
+    sum *= GlobalC::ucell.omega * 0.5;
+
+    return sum;
+}
+
 double Charge_Mixing::inner_product_real(double* rho1, double* rho2)
 {
     double rnorm = 0.0;
diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h
index 124a828a6ef..f101c916b1b 100644
--- a/source/module_elecstate/module_charge/charge_mixing.h
+++ b/source/module_elecstate/module_charge/charge_mixing.h
@@ -33,6 +33,7 @@ class Charge_Mixing
      *
      */
     void mix_rho_recip(Charge* chr);
+    void mix_rho_recip_new(Charge* chr);
 
     /**
      * @brief charge mixing for real space
@@ -45,18 +46,22 @@ class Charge_Mixing
      *
      */
     void Kerker_screen_recip(std::complex<double>* rhog);
+    void Kerker_screen_recip_new(std::complex<double>* rhog);
 
     /**
      * @brief Kerker screen method for real space
      *
      */
     void Kerker_screen_real(double* rho);
+    void Kerker_screen_real_test(double* rho);
 
     /**
      * @brief Inner product of two complex vectors
      *
      */
     double inner_product_recip(std::complex<double>* rho1, std::complex<double>* rho2);
+    double inner_product_recip_new1(std::complex<double>* rho1, std::complex<double>* rho2);
+    double inner_product_recip_new2(std::complex<double>* rho1, std::complex<double>* rho2);
 
     /**
      * @brief Inner product of two double vectors
@@ -126,7 +131,7 @@ class Charge_Mixing
     // General parameters
     //======================================
     std::string mixing_mode = "broyden";
-    double mixing_beta = 0.7;
+    double mixing_beta = 0.8;
     int mixing_ndim = 8;
     double mixing_gg0 = 0.0; // mohan add 2014-09-27
     bool mixing_tau = false;
diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp
index 796246732bf..78e2f65edd6 100644
--- a/source/module_elecstate/test/charge_mixing_test.cpp
+++ b/source/module_elecstate/test/charge_mixing_test.cpp
@@ -22,6 +22,176 @@ Charge::~Charge()
 Charge::Charge()
 {
 }
+
+double Charge::sum_rho(void) const
+{
+    ModuleBase::TITLE("Charge", "sum_rho");
+
+    double sum_rho = 0.0;
+    int nspin0 = (nspin == 2) ? 2 : 1;
+
+    for (int is = 0; is < nspin0; is++)
+    {
+        for (int ir = 0; ir < nrxx; ir++)
+        {
+            if(GlobalV::use_paw)
+            {
+                sum_rho += this->rho[is][ir] + this->nhat[is][ir];
+            }
+            else
+            {
+                sum_rho += this->rho[is][ir];
+            }
+        }
+    }
+
+    // multiply the sum of charge density by a factor
+    sum_rho *= 500.0 / static_cast<double>(this->rhopw->nxyz);
+
+#ifdef __MPI
+    Parallel_Reduce::reduce_pool(sum_rho);
+#endif
+
+    // mohan fixed bug 2010-01-18,
+    // sum_rho may be smaller than 1, like Na bcc.
+    //if (sum_rho <= 0.1)
+    //{
+        //GlobalV::ofs_warning << " sum_rho=" << sum_rho << std::endl;
+        //ModuleBase::WARNING_QUIT("Charge::renormalize_rho", "Can't find even an electron!");
+    //}
+
+    return sum_rho;
+}
+
+void Charge::renormalize_rho(void)
+{
+    ModuleBase::TITLE("Charge", "renormalize_rho");
+
+    const double sr = this->sum_rho();
+    GlobalV::ofs_warning << std::setprecision(15);
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_warning, "charge before normalized", sr);
+    const double normalize_factor = GlobalV::nelec / sr;
+
+    for (int is = 0; is < nspin; is++)
+    {
+        for (int ir = 0; ir < nrxx; ir++)
+        {
+            rho[is][ir] *= normalize_factor;
+        }
+    }
+
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_warning, "charge after normalized", this->sum_rho());
+
+    GlobalV::ofs_running << std::setprecision(6);
+    return;
+}
+
+// for real space magnetic density
+void Charge::get_rho_mag(void)
+{
+    ModuleBase::TITLE("Charge", "get_rho_tot_mag");
+
+    for (int ir = 0; ir < nrxx; ir++)
+    {
+        rho_mag[ir] = rho[0][ir] + rho[1][ir];
+        rho_mag_save[ir] = rho_save[0][ir] + rho_save[1][ir];
+    }
+    for (int ir = 0; ir < nrxx; ir++)
+    {
+        rho_mag[ir + nrxx] = rho[0][ir] - rho[1][ir];
+        rho_mag_save[ir + nrxx] = rho_save[0][ir] - rho_save[1][ir];
+    }
+    return;
+}
+
+void Charge::get_rho_from_mag(void)
+{
+    ModuleBase::TITLE("Charge", "get_rho_from_mag");
+    for (int is = 0; is < nspin; is++)
+    {
+        ModuleBase::GlobalFunc::ZEROS(rho[is], nrxx);
+        //ModuleBase::GlobalFunc::ZEROS(rho_save[is], nrxx);
+    }
+    for (int ir = 0; ir < nrxx; ir++)
+    {
+        rho[0][ir] = 0.5 * (rho_mag[ir] + rho_mag[ir+nrxx]);
+        rho[1][ir] = 0.5 * (rho_mag[ir] - rho_mag[ir+nrxx]);
+    }
+
+    return;
+}
+
+void Charge::allocate_rho_mag(void)
+{
+    rho_mag = new double[nrxx * nspin];
+    rho_mag_save = new double[nrxx * nspin];
+
+    ModuleBase::GlobalFunc::ZEROS(rho_mag, nrxx * nspin);
+    ModuleBase::GlobalFunc::ZEROS(rho_mag_save, nrxx * nspin);
+
+    return;
+}
+
+void Charge::destroy_rho_mag(void)
+{
+    delete[] rho_mag;
+    delete[] rho_mag_save;
+
+    return;
+}
+
+// for reciprocal space magnetic density
+void Charge::get_rhog_mag(void)
+{
+    ModuleBase::TITLE("Charge", "get_rhog_tot_mag");
+
+    for (int ig = 0; ig < ngmc; ig++)
+    {
+        rhog_mag[ig] = rhog[0][ig] + rhog[1][ig];
+        rhog_mag_save[ig] = rhog_save[0][ig] + rhog_save[1][ig];
+    }
+    for (int ig = 0; ig < ngmc; ig++)
+    {
+        rhog_mag[ig + ngmc] = rhog[0][ig] - rhog[1][ig];
+        rhog_mag_save[ig + ngmc] = rhog_save[0][ig] - rhog_save[1][ig];
+    }
+    return;
+}
+
+void Charge::get_rhog_from_mag(void)
+{
+    ModuleBase::TITLE("Charge", "get_rhog_from_mag");
+    for (int is = 0; is < nspin; is++)
+    {
+        ModuleBase::GlobalFunc::ZEROS(rhog[is], ngmc);
+    }
+    for (int ig = 0; ig < ngmc; ig++)
+    {
+        rhog[0][ig] = 0.5 * (rhog_mag[ig] + rhog_mag[ig+ngmc]);
+        rhog[1][ig] = 0.5 * (rhog_mag[ig] - rhog_mag[ig+ngmc]);
+    }
+
+    return;
+}
+
+void Charge::allocate_rhog_mag(void)
+{
+    rhog_mag = new std::complex<double>[ngmc * nspin];
+    rhog_mag_save = new std::complex<double>[ngmc * nspin];
+
+    ModuleBase::GlobalFunc::ZEROS(rhog_mag, ngmc * nspin);
+    ModuleBase::GlobalFunc::ZEROS(rhog_mag_save, ngmc * nspin);
+
+    return;
+}
+
+void Charge::destroy_rhog_mag(void)
+{
+    delete[] rhog_mag;
+    delete[] rhog_mag_save;
+
+    return;
+}
 void Charge::set_rhopw(ModulePW::PW_Basis* rhopw_in)
 {
     this->rhopw = rhopw_in;
@@ -136,6 +306,7 @@ TEST_F(ChargeMixingTest, SetMixingTest)
     EXPECT_THAT(output, testing::HasSubstr("This Mixing mode is not implemended yet,coming soon."));
 }
 
+/**
 TEST_F(ChargeMixingTest, AutoSetTest)
 {
     Charge_Mixing CMtest;
@@ -167,6 +338,7 @@ TEST_F(ChargeMixingTest, AutoSetTest)
     CMtest.auto_set(1.0, GlobalC::ucell);
     EXPECT_EQ(CMtest.mixing_gg0, 1.0);
 }
+**/
 
 TEST_F(ChargeMixingTest, KerkerScreenRecipTest)
 {
@@ -224,6 +396,66 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipTest)
     delete[] drhor_ref;
 }
 
+TEST_F(ChargeMixingTest, KerkerScreenRecipNewTest)
+{
+    Charge_Mixing CMtest;
+    CMtest.set_rhopw(&pw_basis);
+    GlobalC::ucell.tpiba = 1.0;
+
+    // for new kerker
+    GlobalV::NSPIN = 2;
+    CMtest.mixing_gg0 = 0.0;
+    GlobalV::MIXING_GG0_MAG = 0.0;
+    std::complex<double>* drhog = new std::complex<double>[GlobalV::NSPIN*pw_basis.npw];
+    std::complex<double>* drhog_old = new std::complex<double>[GlobalV::NSPIN*pw_basis.npw];
+    double* drhor = new double[GlobalV::NSPIN*pw_basis.nrxx];
+    double* drhor_ref = new double[GlobalV::NSPIN*pw_basis.nrxx];
+    for (int i = 0; i < GlobalV::NSPIN*pw_basis.npw; ++i)
+    {
+        drhog_old[i] = drhog[i] = std::complex<double>(1.0, 1.0);
+    }
+    CMtest.Kerker_screen_recip_new(drhog);
+    for (int i = 0; i < GlobalV::NSPIN*pw_basis.npw; ++i)
+    {
+        EXPECT_EQ(drhog[i], drhog_old[i]);
+    }
+
+    // RECIPROCAL
+    CMtest.mixing_gg0 = 1.0;
+    GlobalV::MIXING_GG0_MAG = 0.0;
+    CMtest.Kerker_screen_recip_new(drhog);
+    const double gg0 = std::pow(0.529177, 2);
+    for (int i = 0; i < pw_basis.npw; ++i)
+    {
+        std::complex<double> ration = drhog[i] / drhog[i+pw_basis.npw];
+        double gg = this->pw_basis.gg[i];
+        double ration_ref = std::max(gg / (gg + gg0), 0.1 / CMtest.mixing_beta);
+        EXPECT_NEAR(ration.real(), ration_ref, 1e-10);
+        EXPECT_NEAR(ration.imag(), 0, 1e-10);
+    }
+
+    // REAL
+    pw_basis.recip2real(drhog, drhor_ref);
+    pw_basis.recip2real(drhog_old, drhor);
+
+    CMtest.mixing_gg0 = 0.0;
+    GlobalV::MIXING_GG0_MAG = 0.0;
+    // nothing happens
+    CMtest.Kerker_screen_real(drhor);
+
+    CMtest.mixing_gg0 = 1.0;
+    CMtest.Kerker_screen_real(drhor);
+    for (int i = 0; i < pw_basis.nrxx; ++i)
+    {
+        EXPECT_NEAR(drhor[i], drhor_ref[i], 1e-8);
+    }
+
+    delete[] drhog;
+    delete[] drhog_old;
+    delete[] drhor;
+    delete[] drhor_ref;
+}
+
 TEST_F(ChargeMixingTest, InnerDotTest)
 {
     // REAL
@@ -293,6 +525,50 @@ TEST_F(ChargeMixingTest, InnerDotTest)
     EXPECT_NEAR(inner, 110668.61166927818, 1e-8);
 }
 
+TEST_F(ChargeMixingTest, InnerDotNewTest)
+{
+    Charge_Mixing CMtest;
+    CMtest.set_rhopw(&pw_basis);
+    GlobalV::NSPIN = 1;
+
+    // inner_product_recip_new1
+    std::vector<std::complex<double>> drhog1(pw_basis.npw);
+    std::vector<std::complex<double>> drhog2(pw_basis.npw);
+    for (int i = 0; i < pw_basis.npw; ++i)
+    {
+        drhog1[i] = 1.0;
+        drhog2[i] = double(i);
+    }
+    double inner = CMtest.inner_product_recip_new1(drhog1.data(), drhog2.data());
+    EXPECT_NEAR(inner, 0.5 * pw_basis.npw * (pw_basis.npw - 1), 1e-8);
+
+    // inner_product_recip_new2
+    GlobalV::NSPIN = 2;
+    drhog1.resize(pw_basis.npw * GlobalV::NSPIN);
+    drhog2.resize(pw_basis.npw * GlobalV::NSPIN);
+    std::vector<std::complex<double>> drhog1_mag(pw_basis.npw * GlobalV::NSPIN);
+    std::vector<std::complex<double>> drhog2_mag(pw_basis.npw * GlobalV::NSPIN);
+    for (int i = 0; i < pw_basis.npw * GlobalV::NSPIN; ++i)
+    {
+        drhog1[i] = std::complex<double>(1.0, double(i));
+        drhog2[i] = std::complex<double>(1.0, 1.0);
+    }
+    // set mag
+    for (int i = 0; i < pw_basis.npw; ++i)
+    {
+        drhog1_mag[i] = drhog1[i] + drhog1[i+pw_basis.npw];
+        drhog1_mag[i+pw_basis.npw] = drhog1[i] - drhog1[i+pw_basis.npw];
+        drhog2_mag[i] = drhog2[i] + drhog2[i+pw_basis.npw];
+        drhog2_mag[i+pw_basis.npw] = drhog2[i] - drhog2[i+pw_basis.npw];
+    }
+    GlobalV::GAMMA_ONLY_PW = false;
+    inner = CMtest.inner_product_recip_new2(drhog1_mag.data(), drhog2_mag.data());
+    EXPECT_NEAR(inner, 236763.82650318215, 1e-8);
+    GlobalV::GAMMA_ONLY_PW = true;
+    inner = CMtest.inner_product_recip_new2(drhog1_mag.data(), drhog2_mag.data());
+    EXPECT_NEAR(inner, 236763.82650318215 * 2, 1e-8);
+}
+
 TEST_F(ChargeMixingTest, MixRhoTest)
 {
     GlobalV::double_grid = false;
diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
index 52351b38b71..a16b24065cd 100644
--- a/source/module_esolver/esolver_ks.cpp
+++ b/source/module_esolver/esolver_ks.cpp
@@ -79,14 +79,19 @@ namespace ModuleESolver
                              GlobalV::MIXING_NDIM,
                              GlobalV::MIXING_GG0,
                              GlobalV::MIXING_TAU);
+        // I use default value to replace autoset                     
         // using bandgap to auto set mixing_beta
-        if (std::abs(GlobalV::MIXING_BETA + 10.0) < 1e-6)
+        // if (std::abs(GlobalV::MIXING_BETA + 10.0) < 1e-6)
+        //{
+        //    p_chgmix->need_auto_set();
+        //}
+        if (GlobalV::MIXING_BETA > 1.0 || GlobalV::MIXING_BETA < 0.0)
         {
-            p_chgmix->need_auto_set();
+            ModuleBase::WARNING("INPUT", "You'd better set mixing_beta to [0.0, 1.0]!");
         }
-        else if (GlobalV::MIXING_BETA > 1.0 || GlobalV::MIXING_BETA < 0.0)
+        if (GlobalV::MIXING_BETA_MAG < 0.0)
         {
-            ModuleBase::WARNING("INPUT", "You'd better set mixing_beta to [0.0, 1.0]!");
+            ModuleBase::WARNING("INPUT", "You'd better set mixing_beta_mag >= 0.0!");
         }
         
 #ifdef USE_PAW
diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp
index a474d6f0761..de1d0003a86 100644
--- a/source/module_io/input.cpp
+++ b/source/module_io/input.cpp
@@ -297,9 +297,11 @@ void Input::Default(void)
     //  charge mixing
     //----------------------------------------------------------
     mixing_mode = "broyden";
-    mixing_beta = -10.0;
+    mixing_beta = -10;
     mixing_ndim = 8;
-    mixing_gg0 = 0.00; // used in kerker method. mohan add 2014-09-27
+    mixing_gg0 = 1.00; // use Kerker defaultly
+    mixing_beta_mag = -10.0; // only set when nspin == 2
+    mixing_gg0_mag = 0.0; // defaultly exclude Kerker from mixing magnetic density
     mixing_tau = false;
     mixing_dftu = false;
     //----------------------------------------------------------
@@ -1230,6 +1232,14 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, mixing_gg0);
         }
+        else if (strcmp("mixing_beta_mag", word) == 0)
+        {
+            read_value(ifs, mixing_beta_mag);
+        }
+        else if (strcmp("mixing_gg0_mag", word) == 0)
+        {
+            read_value(ifs, mixing_gg0_mag);
+        }
         else if (strcmp("mixing_tau", word) == 0)
         {
             read_bool(ifs, mixing_tau);
@@ -2951,6 +2961,38 @@ void Input::Default_2(void) // jiyy add 2019-08-04
             scf_thr_type = 1;
         }
     }
+    // mixing parameters
+    if (mixing_beta < 0.0)
+    {
+        if (nspin == 1)
+        {
+            mixing_beta = 0.8;
+        }
+        else if (nspin == 2)
+        {
+            mixing_beta = 0.4;
+            mixing_beta_mag = 1.6;
+            mixing_gg0_mag = 0.0;
+        }
+        else if (nspin == 4) // I will add this 
+        {
+            mixing_beta = 0.2;
+        }     
+    }
+    else
+    {
+        if (nspin == 2 && mixing_beta_mag < 0.0)
+        {
+            if (mixing_beta <= 0.4)
+            {
+                mixing_beta_mag = 4 * mixing_beta;
+            }
+            else
+            {
+                mixing_beta_mag = 1.6;
+            }
+        }
+    }
 }
 #ifdef __MPI
 void Input::Bcast()
@@ -3106,6 +3148,8 @@ void Input::Bcast()
     Parallel_Common::bcast_double(mixing_beta);
     Parallel_Common::bcast_int(mixing_ndim);
     Parallel_Common::bcast_double(mixing_gg0); // mohan add 2014-09-27
+    Parallel_Common::bcast_double(mixing_beta_mag);
+    Parallel_Common::bcast_double(mixing_gg0_mag);
     Parallel_Common::bcast_bool(mixing_tau);
     Parallel_Common::bcast_bool(mixing_dftu);
 
diff --git a/source/module_io/input.h b/source/module_io/input.h
index 252d161bc18..1213c74889f 100644
--- a/source/module_io/input.h
+++ b/source/module_io/input.h
@@ -230,6 +230,9 @@ class Input
     double mixing_beta; // 0 : no_mixing
     int mixing_ndim; // used in Broyden method
     double mixing_gg0; // used in kerker method. mohan add 2014-09-27
+    double mixing_beta_mag;
+    double mixing_gg0_mag;
+
     bool mixing_tau; // whether to mix tau in mgga
     bool mixing_dftu; //whether to mix locale in DFT+U
 
diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp
index b8306372085..a65b214873d 100644
--- a/source/module_io/input_conv.cpp
+++ b/source/module_io/input_conv.cpp
@@ -741,6 +741,8 @@ void Input_Conv::Convert(void)
     GlobalV::MIXING_BETA = INPUT.mixing_beta;
     GlobalV::MIXING_NDIM = INPUT.mixing_ndim;
     GlobalV::MIXING_GG0 = INPUT.mixing_gg0;
+    GlobalV::MIXING_BETA_MAG = INPUT.mixing_beta_mag;
+    GlobalV::MIXING_GG0_MAG = INPUT.mixing_gg0_mag;
     GlobalV::MIXING_TAU = INPUT.mixing_tau;
     
     ModuleBase::timer::tick("Input_Conv", "Convert");
diff --git a/source/module_io/parameter_pool.cpp b/source/module_io/parameter_pool.cpp
index 8b9cc08132f..aa4a73e3aa9 100644
--- a/source/module_io/parameter_pool.cpp
+++ b/source/module_io/parameter_pool.cpp
@@ -827,6 +827,14 @@ bool input_parameters_set(std::map<std::string, InputParameter> input_parameters
     {
         INPUT.mixing_gg0 = *static_cast<double*>(input_parameters["mixing_gg0"].get());
     }
+    else if (input_parameters.count("mixing_beta_mag") != 0)
+    {
+        INPUT.mixing_beta_mag = *static_cast<double*>(input_parameters["mixing_beta_mag"].get());
+    }
+    else if (input_parameters.count("mixing_gg0_mag") != 0)
+    {
+        INPUT.mixing_gg0_mag = *static_cast<double*>(input_parameters["mixing_gg0_mag"].get());
+    }
     else if (input_parameters.count("mixing_tau") != 0)
     {
         INPUT.mixing_tau = *static_cast<bool*>(input_parameters["mixing_tau"].get());
diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp
index b36fe3dbf3a..d931d34ae9c 100644
--- a/source/module_io/test/input_test.cpp
+++ b/source/module_io/test/input_test.cpp
@@ -155,7 +155,7 @@ TEST_F(InputTest, Default)
         EXPECT_EQ(INPUT.mixing_mode,"broyden");
         EXPECT_DOUBLE_EQ(INPUT.mixing_beta,-10.0);
         EXPECT_EQ(INPUT.mixing_ndim,8);
-        EXPECT_DOUBLE_EQ(INPUT.mixing_gg0,0.00);
+        EXPECT_DOUBLE_EQ(INPUT.mixing_gg0,1.00);
         EXPECT_EQ(INPUT.init_wfc,"atomic");
         EXPECT_EQ(INPUT.mem_saver,0);
         EXPECT_EQ(INPUT.printe,100);
diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp
index 2f592331d46..14613e950b4 100644
--- a/source/module_io/test/input_test_para.cpp
+++ b/source/module_io/test/input_test_para.cpp
@@ -161,7 +161,7 @@ TEST_F(InputParaTest,Bcast)
         EXPECT_EQ(INPUT.mixing_mode,"broyden");
         EXPECT_DOUBLE_EQ(INPUT.mixing_beta,-10.0);
         EXPECT_EQ(INPUT.mixing_ndim,8);
-        EXPECT_DOUBLE_EQ(INPUT.mixing_gg0,0.00);
+        EXPECT_DOUBLE_EQ(INPUT.mixing_gg0,1.00);
         EXPECT_EQ(INPUT.init_wfc,"atomic");
         EXPECT_EQ(INPUT.mem_saver,0);
         EXPECT_EQ(INPUT.printe,100);
diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp
index 00b15c1b39d..e90f461218b 100644
--- a/source/module_io/write_input.cpp
+++ b/source/module_io/write_input.cpp
@@ -246,6 +246,8 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_beta", mixing_beta, "mixing parameter: 0 means no new charge");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_ndim", mixing_ndim, "mixing dimension in pulay or broyden");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_gg0", mixing_gg0, "mixing parameter in kerker");
+    ModuleBase::GlobalFunc::OUTP(ofs, "mixing_beta_mag", mixing_beta_mag, "mixing parameter for magnetic density");
+    ModuleBase::GlobalFunc::OUTP(ofs, "mixing_gg0_mag", mixing_gg0_mag, "mixing parameter in kerker");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_tau", mixing_tau, "whether to mix tau in mGGA calculation");
     ModuleBase::GlobalFunc::OUTP(ofs, "mixing_dftu", mixing_dftu, "whether to mix locale in DFT+U calculation");
 
diff --git a/tests/deepks/603_NO_deepks_ethanol/INPUT b/tests/deepks/603_NO_deepks_ethanol/INPUT
index 05a4fc3595a..4e3362ec4aa 100644
--- a/tests/deepks/603_NO_deepks_ethanol/INPUT
+++ b/tests/deepks/603_NO_deepks_ethanol/INPUT
@@ -25,6 +25,7 @@ smearing_sigma			0.02
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.4
+mixing_gg0      0.0
 
 #Parameters (6.File)
 deepks_out_labels          0 
diff --git a/tests/deepks/603_NO_deepks_watertrimer/INPUT b/tests/deepks/603_NO_deepks_watertrimer/INPUT
index 86e6cd48ea1..750f55a04fc 100644
--- a/tests/deepks/603_NO_deepks_watertrimer/INPUT
+++ b/tests/deepks/603_NO_deepks_watertrimer/INPUT
@@ -12,6 +12,7 @@ dft_functional lda
 gamma_only 1
 mixing_type broyden
 mixing_beta 0.400000
+mixing_gg0  0.0
 nbands 18
 smearing_method gaussian
 smearing_sigma 0.020000
diff --git a/tests/integrate/101_PW_15_lowz/INPUT b/tests/integrate/101_PW_15_lowz/INPUT
index fda37553da0..ffb71c46e9e 100644
--- a/tests/integrate/101_PW_15_lowz/INPUT
+++ b/tests/integrate/101_PW_15_lowz/INPUT
@@ -24,6 +24,7 @@ smearing_sigma		0.1
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 
 nx                      24
 ny                      24
diff --git a/tests/integrate/101_PW_15_paw/INPUT b/tests/integrate/101_PW_15_paw/INPUT
index e9b3ce53bb6..33a6b82d919 100644
--- a/tests/integrate/101_PW_15_paw/INPUT
+++ b/tests/integrate/101_PW_15_paw/INPUT
@@ -26,3 +26,4 @@ smearing_sigma			0.002
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
diff --git a/tests/integrate/107_PW_OK/INPUT b/tests/integrate/107_PW_OK/INPUT
index de663ebf153..83905206921 100644
--- a/tests/integrate/107_PW_OK/INPUT
+++ b/tests/integrate/107_PW_OK/INPUT
@@ -26,4 +26,5 @@ smearing_sigma			0.002
 #Parameters (5.Mixing)
 mixing_type		plain
 mixing_beta		0.7
+mixing_gg0      0.0
 
diff --git a/tests/integrate/186_PW_KG_100/INPUT b/tests/integrate/186_PW_KG_100/INPUT
index 3bba8e5f8a0..a51bb74d946 100644
--- a/tests/integrate/186_PW_KG_100/INPUT
+++ b/tests/integrate/186_PW_KG_100/INPUT
@@ -26,6 +26,7 @@ smearing_sigma			0.6
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.4
+mixing_gg0      0.0
 
 cal_cond      1
 cond_fwhm   8
diff --git a/tests/integrate/186_PW_NLKG_100/INPUT b/tests/integrate/186_PW_NLKG_100/INPUT
index a7f6344438d..bb86d547444 100644
--- a/tests/integrate/186_PW_NLKG_100/INPUT
+++ b/tests/integrate/186_PW_NLKG_100/INPUT
@@ -24,6 +24,7 @@ smearing_sigma			0.6
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.4
+mixing_gg0      0.0
 
 cal_cond      1
 cond_fwhm   8
diff --git a/tests/integrate/186_PW_NSCF_KG_100/INPUT b/tests/integrate/186_PW_NSCF_KG_100/INPUT
index 7057256c29a..7427b20472b 100644
--- a/tests/integrate/186_PW_NSCF_KG_100/INPUT
+++ b/tests/integrate/186_PW_NSCF_KG_100/INPUT
@@ -26,6 +26,7 @@ smearing_sigma			0.6
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.4
+mixing_gg0      0.0
 
 cal_cond      1
 cond_fwhm   8
diff --git a/tests/integrate/186_PW_SDOS_10D10S/INPUT b/tests/integrate/186_PW_SDOS_10D10S/INPUT
index 3d44dd1c770..838e2c143b1 100644
--- a/tests/integrate/186_PW_SDOS_10D10S/INPUT
+++ b/tests/integrate/186_PW_SDOS_10D10S/INPUT
@@ -33,6 +33,7 @@ smearing_sigma			0.6
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.4
+mixing_gg0      0.0
 
 out_dos          1
 dos_emin_ev      -20
diff --git a/tests/integrate/186_PW_SDOS_MALL/INPUT b/tests/integrate/186_PW_SDOS_MALL/INPUT
index cd0b58d85d8..791f415caf9 100644
--- a/tests/integrate/186_PW_SDOS_MALL/INPUT
+++ b/tests/integrate/186_PW_SDOS_MALL/INPUT
@@ -33,6 +33,7 @@ smearing_sigma			0.6
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.4
+mixing_gg0      0.0
 
 out_dos          1
 dos_emin_ev      -20
diff --git a/tests/integrate/186_PW_SKG_10D10S/INPUT b/tests/integrate/186_PW_SKG_10D10S/INPUT
index 4887c16356a..6c664a6fc55 100644
--- a/tests/integrate/186_PW_SKG_10D10S/INPUT
+++ b/tests/integrate/186_PW_SKG_10D10S/INPUT
@@ -33,6 +33,7 @@ smearing_sigma			0.6
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.4
+mixing_gg0      0.0
 
 cal_cond      1
 cond_fwhm   8
diff --git a/tests/integrate/186_PW_SKG_ALL/INPUT b/tests/integrate/186_PW_SKG_ALL/INPUT
index ca3ceb016d2..3d218ef2283 100644
--- a/tests/integrate/186_PW_SKG_ALL/INPUT
+++ b/tests/integrate/186_PW_SKG_ALL/INPUT
@@ -32,6 +32,7 @@ smearing_sigma			0.6
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.4
+mixing_gg0      0.0
 
 cal_cond      1
 cond_fwhm   16
diff --git a/tests/integrate/186_PW_SNLKG_10D10S/INPUT b/tests/integrate/186_PW_SNLKG_10D10S/INPUT
index 90acac914f7..440e0ad9e14 100644
--- a/tests/integrate/186_PW_SNLKG_10D10S/INPUT
+++ b/tests/integrate/186_PW_SNLKG_10D10S/INPUT
@@ -34,6 +34,7 @@ smearing_sigma			0.6
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.4
+mixing_gg0      0.0
 
 cal_cond      1
 cond_fwhm   8
diff --git a/tests/integrate/201_NO_15_f_pseudopots/INPUT b/tests/integrate/201_NO_15_f_pseudopots/INPUT
index 6b52938fe80..1e9185df942 100644
--- a/tests/integrate/201_NO_15_f_pseudopots/INPUT
+++ b/tests/integrate/201_NO_15_f_pseudopots/INPUT
@@ -23,6 +23,7 @@ smearing_sigma			0.002
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 
 out_pot 2
-out_dm1 1
\ No newline at end of file
+out_dm1 1
diff --git a/tests/integrate/201_NO_KP_DJ_CF_CS_GaAs/INPUT b/tests/integrate/201_NO_KP_DJ_CF_CS_GaAs/INPUT
index 228224e2c1c..7076d40ae49 100644
--- a/tests/integrate/201_NO_KP_DJ_CF_CS_GaAs/INPUT
+++ b/tests/integrate/201_NO_KP_DJ_CF_CS_GaAs/INPUT
@@ -32,6 +32,7 @@ pw_diag_thr                0.00001
 
 mixing_type         broyden 
 mixing_beta         0.4
+mixing_gg0          0.0
 mixing_ndim         5
 
 cal_force 1
diff --git a/tests/integrate/204_NO_KP_FM/INPUT b/tests/integrate/204_NO_KP_FM/INPUT
index 5eb89d942cc..d69629638ab 100644
--- a/tests/integrate/204_NO_KP_FM/INPUT
+++ b/tests/integrate/204_NO_KP_FM/INPUT
@@ -24,5 +24,8 @@ smearing_sigma			0.002
 
 #Parameters (5.Mixing)
 mixing_type		broyden
-mixing_beta		0.7
+mixing_beta		0.4
+mixing_gg0      0.0
+mixing_beta_mag 1.6
+mixing_gg0_mag  0.0
 
diff --git a/tests/integrate/204_NO_KP_NC_deltaspin/INPUT b/tests/integrate/204_NO_KP_NC_deltaspin/INPUT
index 5bf9796ed77..c533c35cce8 100644
--- a/tests/integrate/204_NO_KP_NC_deltaspin/INPUT
+++ b/tests/integrate/204_NO_KP_NC_deltaspin/INPUT
@@ -25,6 +25,7 @@ smearing_sigma     0.07
 
 mixing_type    pulay
 mixing_beta    0.3
+mixing_gg0     0.0
 
 ks_solver    genelpa
 basis_type   lcao
diff --git a/tests/integrate/207_NO_KP_OH2/INPUT b/tests/integrate/207_NO_KP_OH2/INPUT
index ac07aa5885d..6b2e88b08d3 100644
--- a/tests/integrate/207_NO_KP_OH2/INPUT
+++ b/tests/integrate/207_NO_KP_OH2/INPUT
@@ -24,6 +24,7 @@ smearing_sigma			0.002
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 
 out_mat_hs2  1
 ks_solver       scalapack_gvx
diff --git a/tests/integrate/207_NO_KP_OHS_SPIN4/INPUT b/tests/integrate/207_NO_KP_OHS_SPIN4/INPUT
index 1b8443a762f..637638e34e2 100644
--- a/tests/integrate/207_NO_KP_OHS_SPIN4/INPUT
+++ b/tests/integrate/207_NO_KP_OHS_SPIN4/INPUT
@@ -23,6 +23,7 @@ smearing_sigma			0.002
 #Parameters (5.Mixing)
 mixing_type		plain
 mixing_beta		0.7
+mixing_gg0      0.0
 
 out_mat_hs	1
 out_mat_hs2	1
diff --git a/tests/integrate/207_NO_OK/INPUT b/tests/integrate/207_NO_OK/INPUT
index 31d95f76956..820fbf7e9f6 100644
--- a/tests/integrate/207_NO_OK/INPUT
+++ b/tests/integrate/207_NO_OK/INPUT
@@ -27,3 +27,4 @@ smearing_sigma			0.002
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
diff --git a/tests/integrate/220_NO_KP_MD_wfc_out/INPUT b/tests/integrate/220_NO_KP_MD_wfc_out/INPUT
index 88eff5100b3..b83a0a782eb 100644
--- a/tests/integrate/220_NO_KP_MD_wfc_out/INPUT
+++ b/tests/integrate/220_NO_KP_MD_wfc_out/INPUT
@@ -23,6 +23,7 @@ force_thr_ev    1.0e-3
 ks_solver       scalapack_gvx
 mixing_type     broyden
 mixing_beta     0.7
+mixing_gg0      0.0
 
 out_wfc_lcao	1
 out_interval	2
diff --git a/tests/integrate/260_NO_15_PK_PU_AF/INPUT b/tests/integrate/260_NO_15_PK_PU_AF/INPUT
index b49841472d0..63b9f76acc3 100644
--- a/tests/integrate/260_NO_15_PK_PU_AF/INPUT
+++ b/tests/integrate/260_NO_15_PK_PU_AF/INPUT
@@ -24,7 +24,7 @@ smearing_sigma    0.07
 
 mixing_type    broyden
 mixing_beta    0.4
-#mixing_gg0    1.5
+#mixing_gg0    1.0
 
 ks_solver    scalapack_gvx
 basis_type    lcao
diff --git a/tests/integrate/260_NO_15_PK_PU_AF/result.ref b/tests/integrate/260_NO_15_PK_PU_AF/result.ref
index 2c8823dbe53..62f49c80fca 100644
--- a/tests/integrate/260_NO_15_PK_PU_AF/result.ref
+++ b/tests/integrate/260_NO_15_PK_PU_AF/result.ref
@@ -1,3 +1,4 @@
-etotref -7300.8388984694393002
-etotperatomref -1825.2097246174
-totaltimeref 10.12041
+etotref -7304.1083525045141869
+etotperatomref -1826.0270881261
+totaltimeref 29.0044
+28.78
diff --git a/tests/integrate/304_NO_GO_AF/INPUT b/tests/integrate/304_NO_GO_AF/INPUT
index 3922d2a1111..f60ff0064e2 100644
--- a/tests/integrate/304_NO_GO_AF/INPUT
+++ b/tests/integrate/304_NO_GO_AF/INPUT
@@ -23,6 +23,7 @@ smearing_sigma          0.02
 
 mixing_type    broyden
 mixing_beta    0.4
+mixing_gg0     0.0
 
 
 ks_solver     scalapack_gvx
diff --git a/tests/integrate/304_NO_GO_AF_atommag/INPUT b/tests/integrate/304_NO_GO_AF_atommag/INPUT
index a8e580be14e..d72d002acb2 100644
--- a/tests/integrate/304_NO_GO_AF_atommag/INPUT
+++ b/tests/integrate/304_NO_GO_AF_atommag/INPUT
@@ -24,6 +24,7 @@ smearing_sigma          0.02
 
 mixing_type    broyden
 mixing_beta    0.4
+mixing_gg0     0.0
 
 
 ks_solver     scalapack_gvx
diff --git a/tests/integrate/304_NO_GO_FM/INPUT b/tests/integrate/304_NO_GO_FM/INPUT
index 06ab7c782ce..e884b05ce17 100644
--- a/tests/integrate/304_NO_GO_FM/INPUT
+++ b/tests/integrate/304_NO_GO_FM/INPUT
@@ -23,6 +23,6 @@ smearing_method		gauss
 smearing_sigma			0.002
 
 #Parameters (5.Mixing)
-mixing_type		plain
+mixing_type		broyden
 mixing_beta		0.7
 
diff --git a/tests/integrate/304_NO_GO_ocp/INPUT b/tests/integrate/304_NO_GO_ocp/INPUT
index ce20aa0faf6..2c0a7a509ea 100644
--- a/tests/integrate/304_NO_GO_ocp/INPUT
+++ b/tests/integrate/304_NO_GO_ocp/INPUT
@@ -23,7 +23,7 @@ smearing_method		gauss
 smearing_sigma			0.002
 
 #Parameters (5.Mixing)
-mixing_type		plain
+mixing_type		broyden
 mixing_beta		0.7
 
 ocp                           1
diff --git a/tests/integrate/360_NO_15_GO_PU_AF/INPUT b/tests/integrate/360_NO_15_GO_PU_AF/INPUT
index 583b5c9f2de..b0fc0b85d84 100644
--- a/tests/integrate/360_NO_15_GO_PU_AF/INPUT
+++ b/tests/integrate/360_NO_15_GO_PU_AF/INPUT
@@ -23,8 +23,8 @@ cal_stress    0
 #relax_bfgs_init    0.5
 
 mixing_type    broyden
-mixing_beta    0.4
-#mixing_gg0    1.5
+mixing_beta    0.2
+mixing_gg0    1.0
 
 ks_solver    scalapack_gvx
 basis_type    lcao
diff --git a/tests/integrate/360_NO_15_GO_PU_AF/result.ref b/tests/integrate/360_NO_15_GO_PU_AF/result.ref
index b348b95b872..5f30fed3618 100644
--- a/tests/integrate/360_NO_15_GO_PU_AF/result.ref
+++ b/tests/integrate/360_NO_15_GO_PU_AF/result.ref
@@ -1,3 +1,4 @@
-etotref -7296.6620153242829474
-etotperatomref -1824.1655038311
-totaltimeref 8.44589
+etotref -7297.9758315540857438
+etotperatomref -1824.4939578885
+totaltimeref 29.0043
+24.90
diff --git a/tests/integrate/601_NO_TDDFT_CO/INPUT b/tests/integrate/601_NO_TDDFT_CO/INPUT
index 813d27bbdef..9b7c61b9fbb 100755
--- a/tests/integrate/601_NO_TDDFT_CO/INPUT
+++ b/tests/integrate/601_NO_TDDFT_CO/INPUT
@@ -19,6 +19,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          1
diff --git a/tests/integrate/601_NO_TDDFT_CO_occ/INPUT b/tests/integrate/601_NO_TDDFT_CO_occ/INPUT
index 3c5f187fc4b..2ba9cf3c5ec 100755
--- a/tests/integrate/601_NO_TDDFT_CO_occ/INPUT
+++ b/tests/integrate/601_NO_TDDFT_CO_occ/INPUT
@@ -20,6 +20,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 #out_alllog		1
 
diff --git a/tests/integrate/601_NO_TDDFT_H2/INPUT b/tests/integrate/601_NO_TDDFT_H2/INPUT
index 05712539087..4739384c7d5 100755
--- a/tests/integrate/601_NO_TDDFT_H2/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2/INPUT
@@ -19,6 +19,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          1
diff --git a/tests/integrate/601_NO_TDDFT_H2_etrs/INPUT b/tests/integrate/601_NO_TDDFT_H2_etrs/INPUT
index 42c092f82ec..1d9fd788090 100755
--- a/tests/integrate/601_NO_TDDFT_H2_etrs/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_etrs/INPUT
@@ -19,6 +19,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          1
diff --git a/tests/integrate/601_NO_TDDFT_H2_kpoint/INPUT b/tests/integrate/601_NO_TDDFT_H2_kpoint/INPUT
index 9b0879d4ee9..f18dbf13cf3 100755
--- a/tests/integrate/601_NO_TDDFT_H2_kpoint/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_kpoint/INPUT
@@ -19,6 +19,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          1
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_gauss/INPUT b/tests/integrate/601_NO_TDDFT_H2_len_gauss/INPUT
index 7e9cfac1ec9..d0a81a84a70 100755
--- a/tests/integrate/601_NO_TDDFT_H2_len_gauss/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_len_gauss/INPUT
@@ -18,6 +18,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          	1
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/INPUT b/tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/INPUT
index 29d554186cc..c1b53837b09 100755
--- a/tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/INPUT
@@ -19,6 +19,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          	1
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_heavi/INPUT b/tests/integrate/601_NO_TDDFT_H2_len_heavi/INPUT
index c40954ccd4e..a73bc56510e 100755
--- a/tests/integrate/601_NO_TDDFT_H2_len_heavi/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_len_heavi/INPUT
@@ -19,6 +19,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          	1
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_hhg/INPUT b/tests/integrate/601_NO_TDDFT_H2_len_hhg/INPUT
index d6e6360c0ab..34d929852c3 100755
--- a/tests/integrate/601_NO_TDDFT_H2_len_hhg/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_len_hhg/INPUT
@@ -18,6 +18,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          	1
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_trape/INPUT b/tests/integrate/601_NO_TDDFT_H2_len_trape/INPUT
index a57bb8a27b2..eadbf296b5a 100755
--- a/tests/integrate/601_NO_TDDFT_H2_len_trape/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_len_trape/INPUT
@@ -18,6 +18,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          	1
diff --git a/tests/integrate/601_NO_TDDFT_H2_len_trigo/INPUT b/tests/integrate/601_NO_TDDFT_H2_len_trigo/INPUT
index 97acff3ac0e..09a8744e190 100755
--- a/tests/integrate/601_NO_TDDFT_H2_len_trigo/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_len_trigo/INPUT
@@ -18,6 +18,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          	1
diff --git a/tests/integrate/601_NO_TDDFT_H2_oldedm/INPUT b/tests/integrate/601_NO_TDDFT_H2_oldedm/INPUT
index 46017b24e9c..42d446a3cf4 100755
--- a/tests/integrate/601_NO_TDDFT_H2_oldedm/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_oldedm/INPUT
@@ -19,6 +19,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          	1
diff --git a/tests/integrate/601_NO_TDDFT_H2_restart/INPUT b/tests/integrate/601_NO_TDDFT_H2_restart/INPUT
index 1d0af3ecf00..8a50b8a04df 100755
--- a/tests/integrate/601_NO_TDDFT_H2_restart/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_restart/INPUT
@@ -19,6 +19,7 @@ md_nstep		7
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 read_file_dir 		./restart
diff --git a/tests/integrate/601_NO_TDDFT_H2_taylor/INPUT b/tests/integrate/601_NO_TDDFT_H2_taylor/INPUT
index 4a9ad73c0a5..3821ab0ed9a 100755
--- a/tests/integrate/601_NO_TDDFT_H2_taylor/INPUT
+++ b/tests/integrate/601_NO_TDDFT_H2_taylor/INPUT
@@ -19,6 +19,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          1
diff --git a/tests/integrate/601_NO_TDDFT_O3/INPUT b/tests/integrate/601_NO_TDDFT_O3/INPUT
index 2e35e0da3b2..661e8960929 100755
--- a/tests/integrate/601_NO_TDDFT_O3/INPUT
+++ b/tests/integrate/601_NO_TDDFT_O3/INPUT
@@ -20,6 +20,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          1
diff --git a/tests/integrate/601_NO_TDDFT_graphene_kpoint/INPUT b/tests/integrate/601_NO_TDDFT_graphene_kpoint/INPUT
index 2cb0516a56c..3f6ec4298d9 100755
--- a/tests/integrate/601_NO_TDDFT_graphene_kpoint/INPUT
+++ b/tests/integrate/601_NO_TDDFT_graphene_kpoint/INPUT
@@ -22,6 +22,7 @@ md_nstep		2
 
 mixing_type		broyden
 mixing_beta		0.7
+mixing_gg0      0.0
 scf_thr			1.0e-6
 
 cal_stress          1