Refactor: using HSolverPW, Psi, ElecStatePW, HamiltPW in ESolver_KS_P, evc should be deleted later

Author: dyzheng

CMakeLists.txt

@@ -276,6 +276,10 @@ target_link_libraries(${ABACUS_BIN_NAME}
     ri
     driver
     xc
+    hsolver
+    elecstate
+    hamilt
+    psi
     esolver
     -lm
 )

source/module_elecstate/elecstate.h

@@ -11,6 +11,7 @@ namespace elecstate
 class ElecState
 {
   public:
+    ElecState(){};
     virtual void init(Charge *chg_in, // pointer for class Charge
                       const K_Vectors *klist_in,
                       int nk_in, // number of k points
@@ -68,6 +69,8 @@ class ElecState
     // occupation weight for each k-point and band
     ModuleBase::matrix wg;
 
+    std::string classname = "none";
+
   protected:
     // calculate ebands for all k points and all occupied bands
     void calEBand();

source/module_elecstate/elecstate_lcao.h

@@ -24,6 +24,7 @@ class ElecStateLCAO : public ElecState
         this->loc = loc_in;
         this->uhm = uhm_in;
         this->lowf = lowf_in;
+        this->classname = "ElecStateLCAO";
     }
     // void init(Charge* chg_in):charge(chg_in){} override;
 

source/module_elecstate/elecstate_pw.cpp

@@ -10,7 +10,18 @@ namespace elecstate
 void ElecStatePW::psiToRho(const psi::Psi<std::complex<double>>& psi)
 {
     this->calculate_weights();
+
     this->calEBand();
+
+    for(int is=0; is<GlobalV::NSPIN; is++)
+	{
+		ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], this->charge->nrxx);
+		if (XC_Functional::get_func_type() == 3)
+		{
+            ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
+        }
+	}
+
     for (int ik = 0; ik < psi.get_nk(); ++ik)
     {
         psi.fix_k(ik);
@@ -126,6 +137,7 @@ void ElecStatePW::rhoBandK(const psi::Psi<std::complex<double>>& psi)
         }
     }
     else
+    {
         for (int ibnd = 0; ibnd < nbands; ibnd++)
         {
             ///
@@ -171,6 +183,7 @@ void ElecStatePW::rhoBandK(const psi::Psi<std::complex<double>>& psi)
                 }
             }
         }
+    }
 
     return;
 }

source/module_elecstate/elecstate_pw.h

@@ -13,6 +13,7 @@ class ElecStatePW : public ElecState
     ElecStatePW(const PW_Basis* basis_in, Charge* chg_in, int nbands_in) : basis(basis_in)
     {
         init(chg_in, basis_in->Klist, basis_in->Klist->nks, nbands_in);
+        this->classname = "ElecStatePW";
     }
     // void init(Charge* chg_in):charge(chg_in){} override;
 

source/module_esolver/esolver_ks.h

@@ -3,6 +3,9 @@
 #include "./esolver_fp.h"
 #include "string.h"
 #include "fstream"
+#include "module_hsolver/hsolver.h"
+#include "module_hamilt/hamilt.h"
+#include "module_elecstate/elecstate.h"
 // #include "estates.h"
 // #include "h2e.h"
 namespace ModuleESolver
@@ -56,6 +59,9 @@ class ESolver_KS: public ESolver_FP
         void reset_diagethr(std::ofstream &ofs_running, const double hsover_error);
 
 
+    hsolver::HSolver* phsol = nullptr;
+    elecstate::ElecState* pelec = nullptr;
+    hamilt::Hamilt* phami = nullptr;
 
     protected:
         std::string basisname; //PW or LCAO

source/module_esolver/esolver_ks_pw.cpp

@@ -26,6 +26,10 @@
 //-----stress------------------
 #include "../src_pw/stress_pw.h"
 //---------------------------------------------------
+#include "module_hsolver/hsolver_pw.h"
+#include "module_elecstate/elecstate_pw.h"
+#include "module_hamilt/hamilt_pw.h"
+#include "module_hsolver/diago_iter_assist.h"
 
 namespace ModuleESolver
 {
@@ -190,6 +194,41 @@ void ESolver_KS_PW:: beforescf()
     {
         srho.begin(is, GlobalC::CHR,GlobalC::pw, GlobalC::Pgrid, GlobalC::symm);
     }
+    //init Psi, HSolver, ElecState, Hamilt
+    hsolver::DiagoIterAssist::PW_DIAG_NMAX = GlobalV::PW_DIAG_NMAX;
+    hsolver::DiagoIterAssist::PW_DIAG_THR = GlobalV::PW_DIAG_THR;
+    const PW_Basis* pbas = &(GlobalC::pw);
+    if(this->phsol == nullptr)
+    {
+        this->phsol = new hsolver::HSolverPW(pbas);
+    }
+    else if(this->phsol->classname != "HSolverPW")
+    {
+        delete[] this->phsol;
+        this->phsol = new hsolver::HSolverPW(pbas);
+    }
+    this->phsol->method = GlobalV::KS_SOLVER;
+    if(this->pelec == nullptr)
+    {
+        this->pelec = new elecstate::ElecStatePW( pbas, (Charge*)(&(GlobalC::CHR)), GlobalV::NBANDS);
+    }
+    else if(this->pelec->classname != "ElecStatePW")
+    {
+        delete[] this->pelec;
+        this->pelec = new elecstate::ElecStatePW( pbas, (Charge*)(&(GlobalC::CHR)), GlobalV::NBANDS);
+    }
+    Hamilt_PW* hpw = &(GlobalC::hm.hpw);
+    if(this->phami == nullptr)
+    {
+        this->phami = new hamilt::HamiltPW(hpw);
+    }
+    else if(this->phami->classname != "HamiltPW")
+    {
+        delete[] this->phami;
+        this->phami = new hamilt::HamiltPW(hpw);
+    }
+    //initial psi
+    GlobalC::wf.evc_transform_psi();
 } 
 
 void ESolver_KS_PW:: eachiterinit(const int iter)
@@ -217,22 +256,37 @@ void ESolver_KS_PW:: eachiterinit(const int iter)
 //Temporary, it should be replaced by hsolver later.
 void ESolver_KS_PW:: hamilt2density(const int istep, const int iter, const double ethr)
 {
-    GlobalV::PW_DIAG_THR = ethr;
-    this->c_bands(istep,iter);
-    
-    GlobalC::en.eband  = 0.0;
-    GlobalC::en.demet  = 0.0;
-    GlobalC::en.ef     = 0.0;
-    GlobalC::en.ef_up  = 0.0;
-    GlobalC::en.ef_dw  = 0.0;
+    if(this->phsol != nullptr)
+    {
+        // reset energy 
+        this->pelec->eband  = 0.0;
+        this->pelec->demet  = 0.0;
+        this->pelec->ef     = 0.0;
+        GlobalC::en.ef_up  = 0.0;
+        GlobalC::en.ef_dw  = 0.0;
+        // choose if psi should be diag in subspace
+        // be careful that istep start from 0 and iter start from 1
+        if((istep==0||istep==1)&&iter==1) 
+        {
+            hsolver::DiagoIterAssist::need_subspace = false;
+        }
+        else 
+        {
+            hsolver::DiagoIterAssist::need_subspace = true;
+        }
 
-    // calculate weights of each band.
-    Occupy::calculate_weights();
+        hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; 
+        this->phsol->solve(this->phami, GlobalC::wf.psi[0], this->pelec);
 
-    // calculate new charge density according to
-    // new wave functions.
-    // calculate the new eband here.
-    GlobalC::CHR.sum_band();
+        // transform energy for print
+        GlobalC::en.eband = this->pelec->eband;
+        GlobalC::en.demet = this->pelec->demet;
+        GlobalC::en.ef = this->pelec->ef;
+    }
+    else
+    {
+        ModuleBase::WARNING_QUIT("ESolver_KS_PW", "HSolver has not been initialed!");
+    }
 
 	// add exx
 #ifdef __LCAO
@@ -350,6 +404,17 @@ void ESolver_KS_PW:: eachiterfinish(const int iter, const bool conv_elec)
 
 void ESolver_KS_PW::afterscf(const bool conv_elec)
 {
+    //temporary transform psi to evc 
+    // psi back to evc
+    GlobalC::wf.psi_transform_evc();
+    for(int ik=0; ik<this->pelec->ekb.nr; ++ik)
+    {
+        for(int ib=0; ib<this->pelec->ekb.nc; ++ib)
+        {
+            GlobalC::wf.ekb[ik][ib] = this->pelec->ekb(ik, ib);
+            GlobalC::wf.wg(ik, ib) = this->pelec->wg(ik, ib);
+        }
+    }
 #ifdef __LCAO
     if(GlobalC::chi0_hilbert.epsilon)                 // pengfei 2016-11-23
     {

source/module_hamilt/hamilt.h

@@ -25,6 +25,8 @@ class Hamilt
     virtual void matrix(MatrixBlock<std::complex<double>> &hk_in, MatrixBlock<std::complex<double>> &sk_in){return;}
     virtual void matrix(MatrixBlock<double> &hk_in, MatrixBlock<double> &sk_in){return;}
 
+    std::string classname = "none";
+
   protected:
     // array, save operations from each operators
     // would be implemented later

source/module_hamilt/hamilt_lcao.h

@@ -44,6 +44,7 @@ template <typename T, typename T1> class HamiltLCAO : public Hamilt
         this->GK = GK_in;
         this->genH = genH_in;
         this->LM = LM_in;
+        this->classname = "HamiltLCAO";
     }
     //~HamiltLCAO();
 

source/module_hamilt/hamilt_pw.cpp

@@ -8,8 +8,9 @@ void HamiltPW::ch_mock()
     return;
 }
 
-void HamiltPW::hk_mock()
+void HamiltPW::hk_mock(const int ik)
 {
+    this->hpw->init_k(ik);
     return;
 }