Refactor: add allocation of HSolver in esolver_ks_lcao and fix bug of DiagoDavid parameter pass

dyzheng · dyzheng · commit 71e31d7885af · 2022-05-20T16:29:00.000+08:00
diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp
@@ -26,6 +26,10 @@
 #endif
 #include "../src_pw/H_Ewald_pw.h"
 
+#include "module_hamilt/hamilt_lcao.h"
+#include "module_elecstate/elecstate_lcao.h"
+#include "module_hsolver/hsolver_lcao.h"
+
 namespace ModuleESolver
 {
 
@@ -188,6 +192,75 @@ namespace ModuleESolver
             srho.begin(is, GlobalC::CHR, GlobalC::pw, GlobalC::Pgrid, GlobalC::symm);
         }
 
+        //init Psi, HSolver, ElecState, Hamilt
+        if(this->phsol != nullptr)
+        {
+            if(this->phsol->classname != "HSolverLCAO")
+            {
+                delete this->phsol;
+                this->phsol = nullptr;
+            }
+        }
+        else
+        {
+            this->phsol = new hsolver::HSolverLCAO();
+            this->phsol->method = GlobalV::KS_SOLVER;
+        }
+        if(this->pelec != nullptr)
+        {
+            if(this->pelec->classname != "ElecStateLCAO")
+            {
+                delete this->pelec;
+                this->pelec = nullptr;
+            }
+        }
+        else
+        {
+            this->pelec = new elecstate::ElecStateLCAO(
+                    (Charge*)(&(GlobalC::CHR)),
+                    &(GlobalC::kv), 
+                    GlobalC::kv.nks,
+                    GlobalV::NBANDS,
+                    &(this->LOC),
+                    &(this->UHM),
+                    &(this->LOWF));
+        }
+        if(this->phami != nullptr)
+        {
+            if(this->phami->classname != "HamiltLCAO")
+            {
+                delete this->phami;
+                this->phami = nullptr;
+            }
+        }
+        else
+        {
+            // three cases for hamilt class
+            if(GlobalV::GAMMA_ONLY_LOCAL)
+            {
+                this->phami = new hamilt::HamiltLCAO<double, double>(
+                    &(this->UHM.GG),
+                    &(this->UHM.genH),
+                    &(this->LM) );
+            }
+            // non-collinear spin case would not use the second template now, 
+            // would add this feature in the future
+            /*else if(GlobalV::NSPIN==4)
+            {
+                this->phami = new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(
+                    &(this->UHM.GK),
+                    &(this->UHM.genH),
+                    &(this->LM) );
+            }*/ 
+            else
+            {
+                this->phami = new hamilt::HamiltLCAO<std::complex<double>, double>(
+                    &(this->UHM.GK),
+                    &(this->UHM.genH),
+                    &(this->LM) );
+            }
+        }
+
         ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf");
         return;
     }
diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
@@ -199,37 +199,43 @@ namespace ModuleESolver
         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)
+        if(this->phsol != nullptr)
         {
-            this->phsol = new hsolver::HSolverPW(pbas);
+            if(this->phsol->classname != "HSolverPW")
+            {
+                delete this->phsol;
+                this->phsol = nullptr;
+            }
         }
-        else if(this->phsol->classname != "HSolverPW")
+        else
         {
-            delete[] this->phsol;
             this->phsol = new hsolver::HSolverPW(pbas);
+            this->phsol->method = GlobalV::KS_SOLVER;
         }
-        this->phsol->method = GlobalV::KS_SOLVER;
-        if(this->pelec == nullptr)
+        if(this->pelec != nullptr)
         {
-            this->pelec = new elecstate::ElecStatePW( pbas, (Charge*)(&(GlobalC::CHR)), GlobalV::NBANDS);
+            if(this->pelec->classname != "ElecStatePW")
+            {
+                delete this->pelec;
+                this->pelec = nullptr;
+            }
         }
-        else if(this->pelec->classname != "ElecStatePW")
+        else
         {
-            delete[] this->pelec;
             this->pelec = new elecstate::ElecStatePW( pbas, (Charge*)(&(GlobalC::CHR)), GlobalV::NBANDS);
         }
-        Hamilt_PW* hpw = &(GlobalC::hm.hpw);
-        if(this->phami == nullptr)
+        if(this->phami != nullptr)
         {
-            this->phami = new hamilt::HamiltPW(hpw);
+            if(this->phami->classname != "HamiltPW")
+            {
+                delete this->phami;
+                this->phami = nullptr;
+            }
         }
-        else if(this->phami->classname != "HamiltPW")
+        else
         {
-            delete[] this->phami;
-            this->phami = new hamilt::HamiltPW(hpw);
+            this->phami = new hamilt::HamiltPW(&(GlobalC::hm.hpw));
         }
-        //initial psi
-        //GlobalC::wf.evc_transform_psi();
     } 
 
     void ESolver_KS_PW::eachiterinit(const int istep, const int iter)
diff --git a/source/module_hamilt/hamilt_lcao.cpp b/source/module_hamilt/hamilt_lcao.cpp
@@ -11,6 +11,8 @@ namespace hamilt
 template class HamiltLCAO<double, double>;
 // case for nspin<4, multi-k-points
 template class HamiltLCAO<std::complex<double>, double>;
+// case for nspin == 4, non-collinear spin case
+//template class HamiltLCAO<std::complex<double>, std::complex<double>>;
 
 template <typename T, typename T1> void HamiltLCAO<T, T1>::constructFixedReal()
 {
@@ -83,8 +85,9 @@ template <typename T, typename T1> void HamiltLCAO<T, T1>::getMatrix(MatrixBlock
 }
 
 // case for nspin==4
-/*void HamiltLCAO<std::complex<double>, std::complex<double>>::matrix(MatrixBlock<std::complex<double>> hk_in,
-MatrixBlock<std::complex<double>> sk_in)const
+/*template <>
+void HamiltLCAO<std::complex<double>, std::complex<double>>::matrix(MatrixBlock<std::complex<double>> &hk_in,
+                                                      MatrixBlock<std::complex<double>> &sk_in)
 {
     this->getMatrix(hk_in, sk_in);
 }*/
@@ -102,41 +105,16 @@ template <> void HamiltLCAO<double, double>::matrix(MatrixBlock<double> &hk_in,
     this->getMatrix(hk_in, sk_in);
 }
 
-// nspin==4 case not supported yet
-/*void HamiltLCAO<std::complex<double>, std::complex<double>>::constructHamilt(const int iter, const MatrixBlock<double>
-rho)
-{
-    this->constructFixedReal();
-    this->constructUpdateReal();
-}*/
 template <typename T, typename T1>
 void HamiltLCAO<T, T1>::constructHamilt(const int iter, const MatrixBlock<double> rho)
 {
     this->constructFixedReal();
     this->constructUpdateReal();
 }
-/*template<>
-void HamiltLCAO<double, double>::constructHamilt(const int iter, const MatrixBlock<double> rho)
-{
-    this->constructFixedReal();
-    this->constructUpdateReal();
-}*/
-
-/*void HamiltLCAO<std::complex<double>, std::complex<double>>::updateHk(const int ik)
-{
-    this->hk_fixed_mock(ik);
-    this->hk_update_mock(ik);
-};*/
 template <typename T, typename T1> void HamiltLCAO<T, T1>::updateHk(const int ik)
 {
     this->hk_fixed_mock(ik);
     this->hk_update_mock(ik);
 }
-/*template<>
-void HamiltLCAO<double, double>::updateHk(const int ik)
-{
-    this->hk_fixed_mock(ik);
-    this->hk_update_mock(ik);
-}*/
 
 } // namespace hamilt
diff --git a/source/module_hamilt/hamilt_lcao.h b/source/module_hamilt/hamilt_lcao.h
@@ -38,9 +38,15 @@ template <typename T> class LocalMatrix
 template <typename T, typename T1> class HamiltLCAO : public Hamilt
 {
   public:
-    HamiltLCAO(Gint_Gamma* GG_in, Gint_k* GK_in, LCAO_gen_fixedH* genH_in, LCAO_Matrix* LM_in)
+    HamiltLCAO(Gint_Gamma* GG_in, LCAO_gen_fixedH* genH_in, LCAO_Matrix* LM_in)
     {
         this->GG = GG_in;
+        this->genH = genH_in;
+        this->LM = LM_in;
+        this->classname = "HamiltLCAO";
+    }
+    HamiltLCAO(Gint_k* GK_in, LCAO_gen_fixedH* genH_in, LCAO_Matrix* LM_in)
+    {
         this->GK = GK_in;
         this->genH = genH_in;
         this->LM = LM_in;
diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp
@@ -45,6 +45,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, psi::Psi<std::complex<double>>& p
     }
     else if (this->method == "dav")
     {
+        DiagoDavid::PW_DIAG_NDIM = GlobalV::PW_DIAG_NDIM;
         if (pdiagh != nullptr)
         {
             if (pdiagh->method != this->method)

Original file line number	Diff line number	Diff line change
`@@ -199,37 +199,43 @@ namespace ModuleESolver`
`199`	`199`	`hsolver::DiagoIterAssist::PW_DIAG_NMAX = GlobalV::PW_DIAG_NMAX;`
`200`	`200`	`hsolver::DiagoIterAssist::PW_DIAG_THR = GlobalV::PW_DIAG_THR;`
`201`	`201`	`const PW_Basis* pbas = &(GlobalC::pw);`
`202`		`- if(this->phsol == nullptr)`
	`202`	`+ if(this->phsol != nullptr)`
`203`	`203`	`{`
`204`		`- this->phsol = new hsolver::HSolverPW(pbas);`
	`204`	`+ if(this->phsol->classname != "HSolverPW")`
	`205`	`+ {`
	`206`	`+ delete this->phsol;`
	`207`	`+ this->phsol = nullptr;`
	`208`	`+ }`
`205`	`209`	`}`
`206`		`- else if(this->phsol->classname != "HSolverPW")`
	`210`	`+ else`
`207`	`211`	`{`
`208`		`- delete[] this->phsol;`
`209`	`212`	`this->phsol = new hsolver::HSolverPW(pbas);`
	`213`	`+ this->phsol->method = GlobalV::KS_SOLVER;`
`210`	`214`	`}`
`211`		`- this->phsol->method = GlobalV::KS_SOLVER;`
`212`		`- if(this->pelec == nullptr)`
	`215`	`+ if(this->pelec != nullptr)`
`213`	`216`	`{`
`214`		`- this->pelec = new elecstate::ElecStatePW( pbas, (Charge*)(&(GlobalC::CHR)), GlobalV::NBANDS);`
	`217`	`+ if(this->pelec->classname != "ElecStatePW")`
	`218`	`+ {`
	`219`	`+ delete this->pelec;`
	`220`	`+ this->pelec = nullptr;`
	`221`	`+ }`
`215`	`222`	`}`
`216`		`- else if(this->pelec->classname != "ElecStatePW")`
	`223`	`+ else`
`217`	`224`	`{`
`218`		`- delete[] this->pelec;`
`219`	`225`	`this->pelec = new elecstate::ElecStatePW( pbas, (Charge*)(&(GlobalC::CHR)), GlobalV::NBANDS);`
`220`	`226`	`}`
`221`		`- Hamilt_PW* hpw = &(GlobalC::hm.hpw);`
`222`		`- if(this->phami == nullptr)`
	`227`	`+ if(this->phami != nullptr)`
`223`	`228`	`{`
`224`		`- this->phami = new hamilt::HamiltPW(hpw);`
	`229`	`+ if(this->phami->classname != "HamiltPW")`
	`230`	`+ {`
	`231`	`+ delete this->phami;`
	`232`	`+ this->phami = nullptr;`
	`233`	`+ }`
`225`	`234`	`}`
`226`		`- else if(this->phami->classname != "HamiltPW")`
	`235`	`+ else`
`227`	`236`	`{`
`228`		`- delete[] this->phami;`
`229`		`- this->phami = new hamilt::HamiltPW(hpw);`
	`237`	`+ this->phami = new hamilt::HamiltPW(&(GlobalC::hm.hpw));`
`230`	`238`	`}`
`231`		`- //initial psi`
`232`		`- //GlobalC::wf.evc_transform_psi();`
`233`	`239`	`}`
`234`	`240`
`235`	`241`	`void ESolver_KS_PW::eachiterinit(const int istep, const int iter)`
Original file line number	Diff line number	Diff line change
`@@ -45,6 +45,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, psi::Psi<std::complex<double>>& p`
`45`	`45`	`}`
`46`	`46`	`else if (this->method == "dav")`
`47`	`47`	`{`
	`48`	`+ DiagoDavid::PW_DIAG_NDIM = GlobalV::PW_DIAG_NDIM;`
`48`	`49`	`if (pdiagh != nullptr)`
`49`	`50`	`{`
`50`	`51`	`if (pdiagh->method != this->method)`