Feature: a new mixing method for magnetic calculations (#3226)

* improved implementation of kerker preconditioner

* fix bug in charge_mixing_test

* make corresponding changes in Documentation after PR#3133

* extra add

* new mixing method for nspin=2

* add some comment for gg=0

* nspin2 mixing method for pw

* remove autoset and set default

* fix charge_mixing_test & add some documentation of mixing_beta_mag/mixing_gg0_mag

* use old mixing to support double_grid

* fix fatal error in CI

* recover Autotest.sh

* 108

* 108

* revert 101_PW_15_lowz

* revert 107_PW_OK

* revert 186_PW_KG_100

* revert 186_PW_SKG_10D10S

* revert 201_NO_15_f_pseudopots

* revert 201_NO_KP_DJ_CF_CS_GaAs

* revert 204_NO_KP_AFM

* revert 207_NO_KP_OHS_SPIN4

* revert 207_NO_KP_OTdH

* revert 207_NO_OK

* revert 220_NO_KP_MD_wfc_out

* revert 101_PW_15_pseudopots_LibxcLDA

* Autotest.sh

* fix bug in deepks_tests

* fix bugs in io_input_test

* update ref of 260 and 360

* do not use renormalize_rho in pw calculations

* set mixing_gg0=0.0 for 601_TD_DFT tests

* change mixing_type=plain to broyden

* add two_beta to replace NSPIN in mixing class

* add UnitTests of inner_product_recip_new1/2 & kerker_recip_new

Author: WHUweiqingzhou

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

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;
 
 //==========================================================

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

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

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

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