Fix: fix nk_total in class RDMFT depend on the Symmetry::symm_flag

Author: JGHan7

source/module_esolver/esolver_ks.cpp

@@ -401,10 +401,10 @@ namespace ModuleESolver
                     }
 
                     // do rdmft calculation
-                    ModuleBase::matrix E_gradient_wg(this->pelec->wg.nr, this->pelec->wg.nc, true);
+                    ModuleBase::matrix E_gradient_occNum(this->pelec->wg.nr, this->pelec->wg.nc, true);
                     psi::Psi<T> E_gradient_wfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis()); 
-                    double Etotal = this->Run_rdmft(E_gradient_wg, E_gradient_wfc);   // add by jghan 2024-03-16
-
+                    double Etotal = this->Run_rdmft(E_gradient_occNum, E_gradient_wfc);   // add by jghan 2024-03-16
+    
                     // continue;
                 }
                 ModuleBase::timer::tick("RDMFT", "E & Egradient");

source/module_esolver/esolver_ks.h

@@ -51,7 +51,7 @@ namespace ModuleESolver
         // get current step of Ionic simulation
         virtual int getniter() override;
 
-        virtual double Run_rdmft(ModuleBase::matrix& E_gradient_wg, psi::Psi<T>& E_gradient_wfc) { return 0.0; };   // add by jghan, 2024-03-16
+        virtual double Run_rdmft(ModuleBase::matrix& E_gradient_occNum, psi::Psi<T>& E_gradient_wfc) { return 0.0; };   // add by jghan, 2024-03-16
         virtual void update_elec_rdmft(const ModuleBase::matrix& occ_number_in, const psi::Psi<T>& wfc_in) {};
 
     protected:

source/module_esolver/esolver_ks_lcao.cpp

@@ -939,7 +939,7 @@ namespace ModuleESolver
 
 
     //initialize the gradients of Etotal on wg and wfc, and set all elements to 0. 
-    ModuleBase::matrix E_gradient_wg(this->pelec->wg.nr, this->pelec->wg.nc, true);
+    ModuleBase::matrix E_gradient_occNum(this->pelec->wg.nr, this->pelec->wg.nc, true);
     psi::Psi<TK> E_gradient_wfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis()); 
     E_gradient_wfc.zero_out();
     double Etotal_RDMFT = 0.0;
@@ -976,7 +976,7 @@ namespace ModuleESolver
             LM.ParaV,
             occ_number,
             wfc_rdmft,
-            E_gradient_wg,
+            E_gradient_occNum,
             E_gradient_wfc,
             this->kv,
             this->UHM.GG,
@@ -997,7 +997,7 @@ namespace ModuleESolver
             LM.ParaV,
             occ_number,
             *(this->psi),
-            E_gradient_wg,
+            E_gradient_occNum,
             E_gradient_wfc,
             this->kv,
             this->UHM.GK,
@@ -1016,7 +1016,7 @@ namespace ModuleESolver
 
     // // test class rdmft
     // rdmft_solver.update_elec(occ_number, *(this->psi));
-    // this->Run_rdmft(E_gradient_wg, E_gradient_wfc);
+    // this->Run_rdmft(E_gradient_occNum, E_gradient_wfc);
     // std::cout << "\nrdmft_solver: " << "0.0000" << std::endl;
 
     // ModuleBase::timer::tick("RDMFT", "E & Egradient");
@@ -1135,9 +1135,9 @@ namespace ModuleESolver
 
 
 template <typename TK, typename TR>
-double ESolver_KS_LCAO<TK, TR>::Run_rdmft(ModuleBase::matrix& E_gradient_wg, psi::Psi<TK>& E_gradient_wfc)
+double ESolver_KS_LCAO<TK, TR>::Run_rdmft(ModuleBase::matrix& E_gradient_occNum, psi::Psi<TK>& E_gradient_wfc)
 {
-    return rdmft_solver.Run(E_gradient_wg, E_gradient_wfc);
+    return rdmft_solver.Run(E_gradient_occNum, E_gradient_wfc);
 }
 
 

source/module_esolver/esolver_ks_lcao.h

@@ -41,7 +41,7 @@ namespace ModuleESolver
         void nscf() override;
         void get_S();
 
-        virtual double Run_rdmft(ModuleBase::matrix& E_gradient_wg, psi::Psi<TK>& E_gradient_wfc) override;         // add by jghan for rdmft calculation, 2024-03-16
+        virtual double Run_rdmft(ModuleBase::matrix& E_gradient_occNum, psi::Psi<TK>& E_gradient_wfc) override;         // add by jghan for rdmft calculation, 2024-03-16
         virtual void update_elec_rdmft(const ModuleBase::matrix& occ_number_in, const psi::Psi<TK>& wfc_in) override;
 
     protected:

source/module_rdmft/CMakeLists.txt

@@ -1,9 +1,9 @@
 add_library(
     rdmft
     OBJECT
-    rdmft_test.cpp
     rdmft.cpp
     rdmft_tools.cpp
+    rdmft_test.cpp
 )
 
 # add_library(

source/module_rdmft/rdmft.cpp

@@ -13,6 +13,7 @@
 
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h"
 #include "module_elecstate/module_dm/cal_dm_psi.h"
+#include "module_cell/module_symmetry/symmetry.h"
 
 #include <iostream>
 #include <cmath>
@@ -79,10 +80,12 @@ void RDMFT<TK, TR>::init(Gint_Gamma& GG_in, Gint_k& GK_in, Parallel_Orbitals& Pa
     ucell = &ucell_in;
     kv = &kv_in;
     charge = &charge_in;
-    nk_total = kv->nkstot_full;
     XC_func_rdmft = XC_func_rdmft_in;
     alpha_power = alpha_power_in;
 
+    if (ModuleSymmetry::Symmetry::symm_flag == -1) nk_total = kv->nkstot_full;
+    else nk_total = kv->nks;
+
     // XC_func_rdmft = "power"; // just for test
 
     // create desc[] and something about MPI to Eij(nbands*nbands)
@@ -175,11 +178,15 @@ void RDMFT<TK, TR>::update_elec(const ModuleBase::matrix& occ_number_in, const p
         }
     }
 
+    std::cout << "\nrdmft_solver: " << "0.1" << std::endl;
+
     // update wfc
     TK* pwfc_in = &wfc_in(0, 0, 0);
     TK* pwfc = &wfc(0, 0, 0);
     for(int i=0; i<wfc.size(); ++i) pwfc[i] = pwfc_in[i];
 
+    std::cout << "\nrdmft_solver: " << "0.2" << std::endl;
+
     // update charge
     if( GlobalV::GAMMA_ONLY_LOCAL )
     {
@@ -206,22 +213,35 @@ void RDMFT<TK, TR>::update_elec(const ModuleBase::matrix& occ_number_in, const p
     {
         // calculate DMK and DMR
         elecstate::DensityMatrix<TK, double> DM(kv, ParaV, GlobalV::NSPIN);
+        std::cout << "\nrdmft_solver: " << "0.21" << std::endl;
+        std::cout << "\nwfc.get_nk(): " << wfc.get_nk() << std::endl;
+        std::cout << "\nkv->nks: " << kv->nks << std::endl;
+        std::cout << "\nwg.nr: " << wg.nr << std::endl;
         elecstate::cal_dm_psi(ParaV, wg, wfc, DM);
+        std::cout << "\nrdmft_solver: " << "0.22" << std::endl;
         DM.init_DMR(&GlobalC::GridD, &GlobalC::ucell);
+        std::cout << "\nrdmft_solver: " << "0.23" << std::endl;
         DM.cal_DMR();
 
+        std::cout << "\nrdmft_solver: " << "0.3" << std::endl;
+
         // this code is copying from function ElecStateLCAO<TK>::psiToRho(), in elecstate_lcao.cpp
         for (int is = 0; is < GlobalV::NSPIN; is++)
         {
             ModuleBase::GlobalFunc::ZEROS(charge->rho[is], charge->nrxx);
         }
 
+        std::cout << "\nrdmft_solver: " << "0.4" << std::endl;
+
         GK->transfer_DM2DtoGrid(DM.get_DMR_vector());
+        std::cout << "\nrdmft_solver: " << "0.5" << std::endl;
         //double** invaild_ptr = nullptr;   // use invaild_ptr replace loc.DM_R in the future
         Gint_inout inout(loc->DM_R, charge->rho, Gint_Tools::job_type::rho);  // what is Local_Orbital_Charge& loc_in? ///////////////
+        std::cout << "\nrdmft_solver: " << "0.6" << std::endl;
         GK->cal_gint(&inout);
-
+        std::cout << "\nrdmft_solver: " << "0.7" << std::endl;
         charge->renormalize_rho();
+        std::cout << "\nrdmft_solver: " << "0.8" << std::endl;
     }
 }
 
@@ -465,7 +485,8 @@ double RDMFT<TK, TR>::cal_rdmft()
     add_occNum(occ_number, wfcHwfc_TV, wfcHwfc_hartree, wfcHwfc_XC, occNum_wfcHamiltWfc, XC_func_rdmft, alpha_power);
 
     // get the total energy
-    add_wfcHwfc(kv->wk, occ_number, wfcHwfc_TV, wfcHwfc_hartree, wfcHwfc_XC, Etotal_n_k, XC_func_rdmft, alpha_power);
+    // add_wfcHwfc(kv->wk, occ_number, wfcHwfc_TV, wfcHwfc_hartree, wfcHwfc_XC, Etotal_n_k, XC_func_rdmft, alpha_power);
+    add_wfcHwfc(wg, wk_fun_occNum, wfcHwfc_TV, wfcHwfc_hartree, wfcHwfc_XC, Etotal_n_k, XC_func_rdmft, alpha_power);
 
     E_RDMFT[3] = getEnergy(Etotal_n_k);
     Parallel_Reduce::reduce_all(E_RDMFT[3]);
@@ -511,19 +532,29 @@ void RDMFT<TK, TR>::cal_Energy()
 
 
 template <typename TK, typename TR>
-double RDMFT<TK, TR>::Run(ModuleBase::matrix& E_gradient_wg, psi::Psi<TK>& E_gradient_wfc)
+double RDMFT<TK, TR>::Run(ModuleBase::matrix& E_gradient_occNum, psi::Psi<TK>& E_gradient_wfc)
 {
     this->get_V_hartree();
     this->get_V_XC();
     this->cal_rdmft();
     this->cal_Energy();
 
-    E_gradient_wg = (occNum_wfcHamiltWfc);
+    E_gradient_occNum = (occNum_wfcHamiltWfc);
 
     TK* pwfc = &occNum_HamiltWfc(0, 0, 0);
     TK* pwfc_out = &E_gradient_wfc(0, 0, 0);
     for(int i=0; i<wfc.size(); ++i) pwfc_out[i] = pwfc[i];
 
+    // test
+    rdmft::printMatrix_pointer(E_gradient_occNum.nr, E_gradient_occNum.nc, &E_gradient_occNum(0, 0), "E_gradient_occNum");
+    rdmft::printMatrix_pointer(occ_number.nr, occ_number.nc, &occ_number(0, 0), "occ_number");
+    rdmft::printMatrix_pointer(wfcHwfc_TV.nr, wfcHwfc_TV.nc, &wfcHwfc_TV(0, 0), "wfcHwfc_TV");
+    rdmft::printMatrix_pointer(wfcHwfc_hartree.nr, wfcHwfc_hartree.nc, &wfcHwfc_hartree(0, 0), "wfcHwfc_hartree");
+    rdmft::printMatrix_pointer(wfcHwfc_XC.nr, wfcHwfc_XC.nc, &wfcHwfc_XC(0, 0), "wfcHwfc_XC");
+    rdmft::printMatrix_pointer(E_gradient_wfc.get_nbands(), E_gradient_wfc.get_nbasis(), &E_gradient_wfc(0, 0, 0), "E_gradient_wfc(ik=0)");
+    rdmft::printMatrix_pointer(E_gradient_wfc.get_nbands(), E_gradient_wfc.get_nbasis(), &E_gradient_wfc(2, 0, 0), "E_gradient_wfc(ik=2)");
+    // test
+
     return E_RDMFT[3];
 }
 

source/module_rdmft/rdmft.h

@@ -156,7 +156,7 @@ class RDMFT
     void update_elec(const ModuleBase::matrix& occ_number_in, const psi::Psi<TK>& wfc_in);
 
     // do all calculation after update occNum&wfc, get Etotal and the gradient of energy with respect to the occNum&wfc
-    double Run(ModuleBase::matrix& E_gradient_wg, psi::Psi<TK>&E_gradient_wfc);
+    double Run(ModuleBase::matrix& E_gradient_occNum, psi::Psi<TK>&E_gradient_wfc);
 
 
 

source/module_rdmft/rdmft_tools.cpp

@@ -117,47 +117,21 @@ void occNum_Mul_wfcHwfc(const ModuleBase::matrix& occ_number, const ModuleBase::
 }
 
 
-// Default symbol = 0 for the gradient of Etotal with respect to occupation numbers
-// symbol = 1 for the relevant calculation of Etotal
+// for the gradient of Etotal with respect to occupation numbers
 void add_occNum(const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
-                const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft, const double alpha, int symbol)
-{
-    
+                const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft, const double alpha)
+{ 
     occNum_wfcHwfc.zero_out();
-    
-    if( symbol==0 )
-    {
-        occNum_Mul_wfcHwfc(occ_number, wfcHwfc_XC_in, occNum_wfcHwfc, 4, XC_func_rdmft, alpha);
-        occNum_wfcHwfc+=(wfcHwfc_TV_in);
-        occNum_wfcHwfc+=(wfcHwfc_hartree_in);
-    }
-    // else if( symbol==1 )
-    // {
-    //     occNum_Mul_wfcHwfc(wg, wfcHwfc_TV_in, occNum_wfcHwfc);
-    //     occNum_Mul_wfcHwfc(wg, wfcHwfc_hartree_in, occNum_wfcHwfc, 1);
-    //     occNum_Mul_wfcHwfc(wg, wfcHwfc_XC_in, occNum_wfcHwfc, 3, XC_func_rdmft, alpha);
-    // }
-    else std::cout << "\n\n\n!!!!!!\nthere are something wrong when calling rdmft_test() and calculation add_occNum()\n!!!!!!\n\n\n"; 
+    occNum_Mul_wfcHwfc(occ_number, wfcHwfc_XC_in, occNum_wfcHwfc, 4, XC_func_rdmft, alpha);
+    occNum_wfcHwfc+=(wfcHwfc_TV_in);
+    occNum_wfcHwfc+=(wfcHwfc_hartree_in);
 }
 
 
 // do wk*g(occNum)*wfcHwfc and add for TV, hartree, XC. This function just use once, so it can be replace and delete
-void add_wfcHwfc(const std::vector<double>& wk_in, const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
+void add_wfcHwfc(const ModuleBase::matrix& wg, const ModuleBase::matrix& wk_fun_occNum, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
                 const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft, const double alpha)
 {
-    /****** delete when function rdmft_cal() -> class RDMFT ******/
-    ModuleBase::matrix wg(occ_number);
-    for(int ik=0; ik < wg.nr; ++ik)
-    {
-        for(int inb=0; inb < wg.nc; ++inb) wg(ik, inb) *= wk_in[ik];
-    }
-    ModuleBase::matrix wk_fun_occNum(occ_number.nr, occ_number.nc, true);
-    for(int ik=0; ik<wg.nr; ++ik)
-    {
-        for(int inb=0; inb<wg.nc; ++inb) wk_fun_occNum(ik, inb) = wk_in[ik] * occNum_func(occ_number(ik, inb), 2, XC_func_rdmft, alpha);
-    }
-    /****** delete when function rdmft_cal() -> class RDMFT ******/
-
     occNum_wfcHwfc.zero_out();
     occNum_Mul_wfcHwfc(wg, wfcHwfc_TV_in, occNum_wfcHwfc);
     occNum_Mul_wfcHwfc(wg, wfcHwfc_hartree_in, occNum_wfcHwfc, 1);

source/module_rdmft/rdmft_tools.h

@@ -261,11 +261,13 @@ void occNum_Mul_wfcHwfc(const ModuleBase::matrix& occ_number, const ModuleBase::
 // Default symbol = 0 for the gradient of Etotal with respect to occupancy
 // symbol = 1 for the relevant calculation of Etotal
 void add_occNum(const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
-            const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft = "HF", const double alpha = 0.656, int symbol = 0);
+            const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft = "HF", const double alpha = 0.656);
 
 
-// do wk*g(occNum)*wfcHwfc and add for TV, hartree, XC. This function just use once, so it can be replace and delete
-void add_wfcHwfc(const std::vector<double>& wk_in, const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
+// // do wk*g(occNum)*wfcHwfc and add for TV, hartree, XC. This function just use once, so it can be replace and delete
+// void add_wfcHwfc(const std::vector<double>& wk_in, const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
+//                 const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft, const double alpha);
+void add_wfcHwfc(const ModuleBase::matrix& wg, const ModuleBase::matrix& wk_fun_occNum, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
                 const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft, const double alpha);