Feature: mixing density and magnetism {mx,my,mz} separately for non-colinear calculation (#3330)

* new mixing for nspin=4

* revert get_drho

* change default value and open mixing_beta_mag for nspin=4

* support new mixing of nspin=4 also for PW basis

* remove seting nspin in front of seting mixing in Default_2()

* update ref of 104_PW_NC_magnetic

* update INPUT of 204_NO_KP_NC

Author: WHUweiqingzhou

source/module_elecstate/module_charge/charge_mixing.cpp

@@ -38,9 +38,9 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in,
     GlobalV::ofs_running<<"mixing_beta: "<< this->mixing_beta <<std::endl;
     GlobalV::ofs_running<<"mixing_gg0: "<< this->mixing_gg0 <<std::endl;
     GlobalV::ofs_running<<"mixing_gg0_min: "<< GlobalV::MIXING_GG0_MIN <<std::endl;
-    if (GlobalV::NSPIN==2)
+    if (GlobalV::NSPIN==2 || GlobalV::NSPIN==4)
     {
-        GlobalV::ofs_running<<"mixing_beta_mag: "<< GlobalV::MIXING_BETA_MAG <<std::endl;
+        GlobalV::ofs_running<<"mixing_beta_mag: "<< this->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;
@@ -271,7 +271,6 @@ void Charge_Mixing::mix_rho_recip(Charge* chr)
 
 void Charge_Mixing::mix_rho_recip_new(Charge* chr)
 {
-    // not support nspin=4 yet 2023/11/17
     // old support see mix_rho_recip()
     if (GlobalV::double_grid)
     {
@@ -281,7 +280,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
     std::complex<double>* rhog_in = nullptr;
     std::complex<double>* rhog_out = nullptr;
 
-    if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
+    if (GlobalV::NSPIN == 1)
     {
         rhog_in = chr->rhog_save[0];
         rhog_out = chr->rhog[0];    
@@ -316,8 +315,38 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
               };
         this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, twobeta_mix, true);
     }
+    else if (GlobalV::NSPIN == 4)
+    {
+        // I do not merge this part with nspin=2 because the method of nspin=2 is almost done,
+        // while nspin=4 is not finished yet. I will try more methods for nspin=4 in the future. 
+        rhog_in = chr->rhog_save[0];
+        rhog_out = chr->rhog[0];
+        const int npw = this->rhopw->npw;
+        auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1); // use old one
+        auto twobeta_mix
+            = [this, npw](std::complex<double>* out, const std::complex<double>* in, const std::complex<double>* sres) {
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 256)
+#endif
+                  for (int i = 0; i < npw; ++i)
+                  {
+                      out[i] = in[i] + this->mixing_beta * sres[i];
+                  }
+            // magnetism, mx, my, mz
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 256)
+#endif
+                  for (int i = npw; i < 4 * npw; ++i)
+                  {
+                      out[i] = in[i] + this->mixing_beta_mag * sres[i];
+                  }
+              };
+        this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, twobeta_mix, true);
+    }
 
     //  can choose inner_product_recip_new1 or inner_product_recip_new2
+    //  inner_product_recip_new1 is a simple sum
+    //  inner_product_recip_new2 is a hartree-like sum, unit is Ry
     auto inner_product_new
         = std::bind(&Charge_Mixing::inner_product_recip_new2, this, std::placeholders::_1, std::placeholders::_2);
     auto inner_product_old
@@ -326,7 +355,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
     {
         this->mixing->cal_coef(this->rho_mdata, inner_product_new);
     }
-    else if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
+    else if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) // nspin=4 can use old inner_product
     {
         this->mixing->cal_coef(this->rho_mdata, inner_product_old);
     }
@@ -344,6 +373,8 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
     {
         chr->rhopw->recip2real(chr->rhog[is], chr->rho[is]);
     }
+
+    // renormalize rho in R-space would induce a error in K-space
     //chr->renormalize_rho();
 
     // For kinetic energy density
@@ -400,7 +431,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
 {
     double* rhor_in;
     double* rhor_out;
-    if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
+    if (GlobalV::NSPIN == 1)
     {
         rhor_in = chr->rho_save[0];
         rhor_out = chr->rho[0];
@@ -420,19 +451,47 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 256)
 #endif
-                  for (int i = 0; i < nrxx; ++i)
-                  {
-                      out[i] = in[i] + this->mixing_beta * sres[i];
-                  }
+            for (int i = 0; i < nrxx; ++i)
+            {
+                out[i] = in[i] + this->mixing_beta * sres[i];
+            }
             // magnetism
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 256)
 #endif
-                  for (int i = nrxx; i < 2 * nrxx; ++i)
-                  {
-                      out[i] = in[i] + this->mixing_beta_mag * sres[i];
-                  }
-              };
+            for (int i = nrxx; i < 2 * nrxx; ++i)
+            {
+                out[i] = in[i] + this->mixing_beta_mag * sres[i];
+            }
+        };
+        this->mixing->push_data(this->rho_mdata, rhor_in, rhor_out, screen, twobeta_mix, true);
+    }
+    else if (GlobalV::NSPIN == 4)
+    {
+        // I do not merge this part with nspin=2 because the method of nspin=2 is almost done,
+        // while nspin=4 is not finished yet. I will try more methods for nspin=4 in the future. 
+        rhor_in = chr->rho_save[0];
+        rhor_out = chr->rho[0];
+        const int nrxx = this->rhopw->nrxx;
+        auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1);
+        auto twobeta_mix
+            = [this, nrxx](double* out, const double* in, const double* sres) {
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 256)
+#endif
+            for (int i = 0; i < nrxx; ++i)
+            {
+                out[i] = in[i] + this->mixing_beta * sres[i];
+            }
+            // magnetism, mx, my, mz
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 256)
+#endif
+            for (int i = nrxx; i < 4 * nrxx; ++i)
+            {
+                out[i] = in[i] + this->mixing_beta_mag * sres[i];
+            }
+        };
         this->mixing->push_data(this->rho_mdata, rhor_in, rhor_out, screen, twobeta_mix, true);
     }
 
@@ -642,14 +701,18 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex<double>* drhog)
     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 (is >= 1)
         {
             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;
-                }
+#ifdef __DEBUG
+                assert(is == 1); // make sure break works
+#endif
+                double is_mag = GlobalV::NSPIN - 1;
+                //for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++)
+                //{
+                //    drhog[is * this->rhopw->npw + ig] *= 1;
+                //}
                 break;
             }
             fac = GlobalV::MIXING_GG0_MAG;
@@ -692,13 +755,17 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
     for (int is = 0; is < GlobalV::NSPIN; is++)
     {
 
-        if (is == 1 && GlobalV::NSPIN == 2)
+        if (is >= 1)
         {
             if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1)
             {
-                for (int ig = 0; ig < this->rhopw->npw; ig++)
+#ifdef __DEBUG
+                assert(is == 1); // make sure break works
+#endif
+                double is_mag = GlobalV::NSPIN - 1;
+                for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++)
                 {
-                    drhog[is * this->rhopw->npw + ig] *= 0;
+                    drhog[is * this->rhopw->npw + ig] = 0;
                 }
                 break;
             }

source/module_io/input.cpp

@@ -301,7 +301,7 @@ void Input::Default(void)
     mixing_beta = -10;
     mixing_ndim = 8;
     mixing_gg0 = 1.00; // use Kerker defaultly
-    mixing_beta_mag = -10.0; // only set when nspin == 2
+    mixing_beta_mag = -10.0; // only set when nspin == 2 || nspin == 4
     mixing_gg0_mag = 0.0; // defaultly exclude Kerker from mixing magnetic density
     mixing_gg0_min = 0.1; // defaultly minimum kerker coefficient
     mixing_tau = false;
@@ -2975,6 +2975,13 @@ void Input::Default_2(void) // jiyy add 2019-08-04
             scf_thr_type = 1;
         }
     }
+
+    // set nspin with noncolin
+    if (noncolin || lspinorb)
+    {
+        nspin = 4;
+    }
+
     // mixing parameters
     if (mixing_beta < 0.0)
     {
@@ -2990,28 +2997,25 @@ void Input::Default_2(void) // jiyy add 2019-08-04
         }
         else if (nspin == 4) // I will add this 
         {
-            mixing_beta = 0.2;
+            mixing_beta = 0.4;
+            mixing_beta_mag = 1.6;
+            mixing_gg0_mag = 0.0;
         }     
     }
     else
     {
-        if (nspin == 2 && mixing_beta_mag < 0.0)
+        if ((nspin == 2 || nspin == 4) && mixing_beta_mag < 0.0)
         {
             if (mixing_beta <= 0.4)
             {
                 mixing_beta_mag = 4 * mixing_beta;
             }
             else
             {
-                mixing_beta_mag = 1.6;
+                mixing_beta_mag = 1.6; // 1.6 can be discussed
             }
         }
     }
-    // set nspin with noncolin
-    if (noncolin || lspinorb)
-    {
-        nspin = 4;
-    }
 }
 #ifdef __MPI
 void Input::Bcast()

tests/integrate/104_PW_NC_magnetic/result.ref

@@ -1,6 +1,6 @@
-etotref -5866.174730432665
-etotperatomref -2933.0873652163
+etotref -5866.197307043378
+etotperatomref -2933.0986535217
 pointgroupref O_h
 spacegroupref O_h
 nksibzref 1
-totaltimeref 0.77367
+totaltimeref 2.23

tests/integrate/204_NO_KP_NC_deltaspin/INPUT

@@ -25,6 +25,7 @@ smearing_sigma     0.07
 
 mixing_type    pulay
 mixing_beta    0.3
+mixing_beta_mag 0.3
 mixing_gg0     0.0
 
 ks_solver    genelpa