From 09ef2fce428577e32a9daa8faaffb6bc04312c79 Mon Sep 17 00:00:00 2001
From: wenfei-li <liwenfei@gmail.com>
Date: Wed, 8 Nov 2023 19:37:54 +0800
Subject: [PATCH 1/2] Feature : add mixing of compensation charge

---
 .../module_charge/charge_mixing.cpp           | 53 ++++++++++++++++++-
 .../module_charge/charge_mixing.h             |  1 +
 2 files changed, 53 insertions(+), 1 deletion(-)

diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp
index 1ed6df65fe..fefdc128ca 100755
--- a/source/module_elecstate/module_charge/charge_mixing.cpp
+++ b/source/module_elecstate/module_charge/charge_mixing.cpp
@@ -58,6 +58,9 @@ 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));
     }
+
+    if(GlobalV::use_paw) this->mixing->init_mixing_data(this->nhat_mdata, this->rhopw->nrxx * GlobalV::NSPIN, sizeof(double));
+
     // 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 +304,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 +365,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 +404,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 +460,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");
diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h
index 7cfcab4947..124a828a6e 100644
--- a/source/module_elecstate/module_charge/charge_mixing.h
+++ b/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

From 92e83edef30e912050f2f920ef0580d92f8393cd Mon Sep 17 00:00:00 2001
From: wenfei-li <liwenfei@gmail.com>
Date: Thu, 9 Nov 2023 14:10:37 +0800
Subject: [PATCH 2/2] add use_paw

---
 source/module_elecstate/module_charge/charge_mixing.cpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp
index fefdc128ca..8a0912d580 100755
--- a/source/module_elecstate/module_charge/charge_mixing.cpp
+++ b/source/module_elecstate/module_charge/charge_mixing.cpp
@@ -59,7 +59,9 @@ 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));