Feature : add mixing of compensation charge

source/module_elecstate/module_charge/charge_mixing.cpp

@@ -58,6 +58,11 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in,
     {
         this->mixing->init_mixing_data(this->rho_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double));
     }
+
+#ifdef USE_PAW
+    if(GlobalV::use_paw) this->mixing->init_mixing_data(this->nhat_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double));
+#endif
+
     // Note: we can not init tau_mdata here temporarily, since set_xc_type() is after it.
     // this->mixing->init_mixing_data(this->tau_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double));
     return;
@@ -301,6 +306,19 @@ void Charge_Mixing::mix_rho_recip(Charge* chr)
         }
     }
 
+#ifdef USE_PAW
+    if(GlobalV::use_paw)
+    {
+        double *nhat_out, *nhat_in;
+        nhat_in = chr->nhat_save[0];
+        nhat_out = chr->nhat[0];
+        // 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->nhat_mdata, nhat_in, nhat_out, nullptr, false);
+
+        this->mixing->mix_data(this->nhat_mdata, nhat_out);
+    }
+#endif
     return;
 }
 
@@ -349,6 +367,9 @@ void Charge_Mixing::mix_reset()
             this->mixing->init_mixing_data(this->tau_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double));
         }
     }
+#ifdef USE_PAW
+    if(GlobalV::use_paw) this->mixing->init_mixing_data(this->nhat_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double));
+#endif
 }
 
 void Charge_Mixing::mix_rho(Charge* chr)
@@ -385,7 +406,23 @@ void Charge_Mixing::mix_rho(Charge* chr)
             kin_r123[ir] = chr->kin_r[0][ir];
         }
     }
-
+#ifdef USE_PAW
+    std::vector<double> nhat_r123;
+    if(GlobalV::use_paw)
+    {
+        nhat_r123.resize(GlobalV::NSPIN * nrxx);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 512)
+#endif
+        for(int ir = 0 ; ir < nrxx ; ++ir)
+        {
+            for(int is = 0; is < GlobalV::NSPIN; ++is)
+            {
+                nhat_r123[ir+is*nrxx] = chr->nhat[0][ir];
+            }
+        }
+    }        
+#endif
     // --------------------Mixing Body--------------------
     if (GlobalV::SCF_THR_TYPE == 1)
     {
@@ -425,6 +462,22 @@ void Charge_Mixing::mix_rho(Charge* chr)
         }
     }
 
+#ifdef USE_PAW
+    if(GlobalV::use_paw)
+    {
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static, 512)
+#endif
+        for(int ir = 0 ; ir < nrxx ; ++ir)
+        {
+            for(int is = 0; is < GlobalV::NSPIN; ++is)
+            {
+                chr->nhat_save[is][ir] = nhat_r123[ir+is*nrxx];
+            }
+        }
+    }
+#endif
+
     if (new_e_iteration)
         new_e_iteration = false;
     ModuleBase::timer::tick("Charge", "mix_rho");

source/module_elecstate/module_charge/charge_mixing.h

@@ -14,6 +14,7 @@ class Charge_Mixing
     Base_Mixing::Mixing* mixing = nullptr;
     Base_Mixing::Mixing_Data rho_mdata;
     Base_Mixing::Mixing_Data tau_mdata;
+    Base_Mixing::Mixing_Data nhat_mdata;
 
     /**
      * @brief reset mixing