Merge pull request #1002 from deepmodeling/HSolver

Refactor: add HSolver module in LCAO

Author: caic99

source/Makefile.Objects

@@ -136,11 +136,7 @@ ORB_table_alpha.o\
 ORB_gen_tables.o\
 local_orbital_wfc.o\
 local_orbital_charge.o\
-ELEC_cbands_k.o\
-ELEC_cbands_gamma.o\
 ELEC_evolve.o\
-ELEC_scf.o\
-ELEC_nscf.o\
 LOOP_cell.o\
 LOOP_ions.o\
 run_md_lcao.o\

source/input_conv.cpp

@@ -29,6 +29,9 @@
 #include "module_base/timer.h"
 #include "module_surchem/efield.h"
 
+#include "module_elecstate/elecstate_lcao.h"
+#include "module_hsolver/hsolver_lcao.h"
+
 void Input_Conv::Convert(void)
 {
     ModuleBase::TITLE("Input_Conv", "Convert");
@@ -436,9 +439,9 @@ void Input_Conv::Convert(void)
     GlobalC::en.out_proj_band = INPUT.out_proj_band;
 #ifdef __LCAO
     Local_Orbital_Charge::out_dm = INPUT.out_dm;
-    Pdiag_Double::out_mat_hs = INPUT.out_mat_hs;
-    Pdiag_Double::out_mat_hsR = INPUT.out_mat_hs2; // LiuXh add 2019-07-16
-    Pdiag_Double::out_wfc_lcao = INPUT.out_wfc_lcao;
+    hsolver::HSolverLCAO::out_mat_hs = INPUT.out_mat_hs;
+    hsolver::HSolverLCAO::out_mat_hsR = INPUT.out_mat_hs2; // LiuXh add 2019-07-16
+    elecstate::ElecStateLCAO::out_wfc_lcao = INPUT.out_wfc_lcao;
 #endif
 
     GlobalC::en.dos_emin_ev = INPUT.dos_emin_ev;

source/module_elecstate/cal_dm.h

@@ -0,0 +1,115 @@
+#ifndef CAL_DM_H
+#define CAL_DM_H
+
+#include "math_tools.h"
+#include "module_base/matrix.h"
+#include "module_base/complexmatrix.h"
+#include "src_lcao/local_orbital_charge.h"
+
+namespace elecstate
+{
+
+// for gamma_only(double case) and multi-k(complex<double> case)
+inline void cal_dm(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& wg, const psi::Psi<double>& wfc, std::vector<ModuleBase::matrix>& dm)
+{
+    ModuleBase::TITLE("elecstate", "cal_dm");
+
+    //dm.resize(wfc.get_nk(), ParaV->ncol, ParaV->nrow);
+    const int nbands_local = wfc.get_nbands();
+    const int nbasis_local = wfc.get_nbasis();
+
+    // dm = wfc.T * wg * wfc.conj()
+    // dm[is](iw1,iw2) = \sum_{ib} wfc[is](ib,iw1).T * wg(is,ib) * wfc[is](ib,iw2).conj()
+    for (int ik = 0; ik < wfc.get_nk(); ++ik)
+    {
+        wfc.fix_k(ik);
+        //dm.fix_k(ik);
+        dm[ik].create(ParaV->ncol, ParaV->nrow);
+        // wg_wfc(ib,iw) = wg[ib] * wfc(ib,iw);
+        psi::Psi<double> wg_wfc(wfc, 1);
+
+        int ib_global = 0;
+        for (int ib_local = 0; ib_local < nbands_local; ++ib_local)
+        {
+            while (ib_local != ParaV->trace_loc_col[ib_global])
+            {
+                ++ib_global;
+                if (ib_global >= wg.nc)
+                {
+                    break;
+                    ModuleBase::WARNING_QUIT("ElecStateLCAO::cal_dm", "please check trace_loc_col!");
+                }
+            }
+            if (ib_global >= wg.nc) continue;
+            const double wg_local = wg(ik, ib_global);
+            double* wg_wfc_pointer = &(wg_wfc(0, ib_local, 0));
+            BlasConnector::scal(nbasis_local, wg_local, wg_wfc_pointer, 1);
+        }
+
+        // C++: dm(iw1,iw2) = wfc(ib,iw1).T * wg_wfc(ib,iw2)
+#ifdef __MPI
+        psiMulPsiMpi(wg_wfc, wfc, dm[ik], ParaV->desc_wfc, ParaV->desc);
+#else
+        psiMulPsi(wg_wfc, wfc, dm[ik]);
+#endif
+    }
+
+    return;
+}
+
+inline void cal_dm(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& wg, const psi::Psi<std::complex<double>>& wfc, std::vector<ModuleBase::ComplexMatrix>& dm)
+{
+    ModuleBase::TITLE("elecstate", "cal_dm");
+
+    //dm.resize(wfc.get_nk(), ParaV->ncol, ParaV->nrow);
+    const int nbands_local = wfc.get_nbands();
+    const int nbasis_local = wfc.get_nbasis();
+
+    // dm = wfc.T * wg * wfc.conj()
+    // dm[is](iw1,iw2) = \sum_{ib} wfc[is](ib,iw1).T * wg(is,ib) * wfc[is](ib,iw2).conj()
+    for (int ik = 0; ik < wfc.get_nk(); ++ik)
+    {
+        wfc.fix_k(ik);
+        //dm.fix_k(ik);
+        dm[ik].create(ParaV->ncol, ParaV->nrow);
+        // wg_wfc(ib,iw) = wg[ib] * wfc(ib,iw);
+        psi::Psi<std::complex<double>> wg_wfc(1, wfc.get_nbands(), wfc.get_nbasis(), nullptr);
+        const std::complex<double>* pwfc = wfc.get_pointer();
+        std::complex<double>* pwg_wfc = wg_wfc.get_pointer();
+        for(int i = 0;i<wg_wfc.size();++i)
+        {
+            pwg_wfc[i] = conj(pwfc[i]);
+        }
+
+        int ib_global = 0;
+        for (int ib_local = 0; ib_local < nbands_local; ++ib_local)
+        {
+            while (ib_local != ParaV->trace_loc_col[ib_global])
+            {
+                ++ib_global;
+                if (ib_global >= wg.nc)
+                {
+                    break;
+                    ModuleBase::WARNING_QUIT("ElecStateLCAO::cal_dm", "please check trace_loc_col!");
+                }
+            }
+            if (ib_global >= wg.nc) continue;
+            const double wg_local = wg(ik, ib_global);
+            std::complex<double>* wg_wfc_pointer = &(wg_wfc(0, ib_local, 0));
+            BlasConnector::scal(nbasis_local, wg_local, wg_wfc_pointer, 1);
+        }
+
+        // C++: dm(iw1,iw2) = wfc(ib,iw1).T * wg_wfc(ib,iw2)
+#ifdef __MPI
+        psiMulPsiMpi(wg_wfc, wfc, dm[ik], ParaV->desc_wfc, ParaV->desc);
+#else
+        psiMulPsi(wg_wfc, wfc, dm[ik]);
+#endif
+    }
+
+    return;
+}
+
+}//namespace elecstate
+
+#endif

source/module_elecstate/elecstate.cpp

@@ -20,6 +20,18 @@ void ElecState::calculate_weights()
 {
     ModuleBase::TITLE("ElecState", "calculate_weights");
 
+    if (GlobalV::ocp == 1)
+    {
+        for (int ik = 0; ik < GlobalC::kv.nks; ik++)
+        {
+            for (int ib = 0; ib < GlobalV::NBANDS; ib++)
+            {
+                this->wg(ik, ib) = GlobalV::ocp_kb[ik * GlobalV::NBANDS + ib];
+            }
+        }
+        return;
+    }
+
     // for test
     //	std::cout << " gaussian_broadening = " << use_gaussian_broadening << std::endl;
     //	std::cout << " tetrahedron_method = " << use_tetrahedron_method << std::endl;
@@ -191,7 +203,7 @@ void ElecState::calculate_weights()
 
 void ElecState::calEBand()
 {
-    ModuleBase::TITLE("ElecStatePW", "calEBand");
+    ModuleBase::TITLE("ElecState", "calEBand");
     //calculate ebands using wg and ekb
     this->eband = 0.0;
     for (int ik = 0; ik < this->ekb.nr; ++ik)

source/module_elecstate/elecstate_lcao.cpp

@@ -1,58 +1,13 @@
 #include "elecstate_lcao.h"
 
-#include "math_tools.h"
+#include "cal_dm.h"
 #include "module_base/timer.h"
 #include "module_gint/grid_technique.h"
 
 namespace elecstate
 {
 int ElecStateLCAO::out_wfc_lcao = 0;
 
-// for gamma_only(double case) and multi-k(complex<double> case)
-template <typename T> void ElecStateLCAO::cal_dm(const ModuleBase::matrix& wg, const psi::Psi<T>& wfc, psi::Psi<T>& dm)
-{
-    ModuleBase::TITLE("ElecStateLCAO", "cal_dm");
-
-    dm.resize(wfc.get_nk(), this->loc->ParaV->ncol, this->loc->ParaV->nrow);
-    const int nbands_local = wfc.get_nbands();
-    const int nbasis_local = wfc.get_nbasis();
-
-    // dm = wfc.T * wg * wfc.conj()
-    // dm[is](iw1,iw2) = \sum_{ib} wfc[is](ib,iw1).T * wg(is,ib) * wfc[is](ib,iw2).conj()
-    for (int ik = 0; ik < wfc.get_nk(); ++ik)
-    {
-        wfc.fix_k(ik);
-        dm.fix_k(ik);
-        // wg_wfc(ib,iw) = wg[ib] * wfc(ib,iw);
-        psi::Psi<T> wg_wfc(wfc, 1);
-
-        int ib_global = 0;
-        for (int ib_local = 0; ib_local < nbands_local; ++ib_local)
-        {
-            while (ib_local != this->loc->ParaV->trace_loc_col[ib_global])
-            {
-                ++ib_global;
-                if (ib_global >= wg.nc)
-                {
-                    ModuleBase::WARNING_QUIT("ElecStateLCAO::cal_dm", "please check trace_loc_col!");
-                }
-            }
-            const double wg_local = wg(ik, ib_global);
-            T* wg_wfc_pointer = &(wg_wfc(0, ib_local, 0));
-            BlasConnector::scal(nbasis_local, wg_local, wg_wfc_pointer, 1);
-        }
-
-        // C++: dm(iw1,iw2) = wfc(ib,iw1).T * wg_wfc(ib,iw2)
-#ifdef __MPI
-        psiMulPsiMpi(wg_wfc, wfc, dm, this->loc->ParaV->desc_wfc, this->loc->ParaV->desc);
-#else
-        psiMulPsi(wg_wfc, wfc, dm);
-#endif
-    }
-
-    return;
-}
-
 // multi-k case
 void ElecStateLCAO::psiToRho(const psi::Psi<std::complex<double>>& psi)
 {
@@ -65,12 +20,12 @@ void ElecStateLCAO::psiToRho(const psi::Psi<std::complex<double>>& psi)
     ModuleBase::GlobalFunc::NOTE("Calculate the density matrix.");
 
     // this part for calculating dm_k in 2d-block format, not used for charge now
-    psi::Psi<std::complex<double>> dm_k_2d(psi.get_nk(), psi.get_nbasis(), psi.get_nbasis());
+//    psi::Psi<std::complex<double>> dm_k_2d();
 
     if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx"
         || GlobalV::KS_SOLVER == "lapack") // Peize Lin test 2019-05-15
     {
-        this->cal_dm(this->wg, psi, dm_k_2d);
+        cal_dm(this->loc->ParaV, this->wg, psi, this->loc->dm_k);
     }
 
     // this part for steps:
@@ -82,11 +37,29 @@ void ElecStateLCAO::psiToRho(const psi::Psi<std::complex<double>>& psi)
         for (int ik = 0; ik < psi.get_nk(); ik++)
         {
             psi.fix_k(ik);
-            this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi.get_pointer(), this->lowf->wfc_k_grid[ik], ik);
+            this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi.get_pointer(), this->lowf->wfc_k_grid[ik], ik, this->ekb, this->wg);
+            //added by zhengdy-soc, rearrange the wfc_k_grid from [up,down,up,down...] to [up,up...down,down...],
+            if(GlobalV::NSPIN==4)
+            {
+                int row = GlobalC::GridT.lgd;
+                std::vector<std::complex<double>> tmp(row);
+                for(int ib=0; ib<GlobalV::NBANDS; ib++)
+                {
+                    for(int iw=0; iw<row / GlobalV::NPOL; iw++)
+                    {
+                        tmp[iw] = this->lowf->wfc_k_grid[ik][ib][iw * GlobalV::NPOL];
+                        tmp[iw + row / GlobalV::NPOL] = this->lowf->wfc_k_grid[ik][ib][iw * GlobalV::NPOL + 1];
+                    }
+                    for(int iw=0; iw<row; iw++)
+                    {
+                        this->lowf->wfc_k_grid[ik][ib][iw] = tmp[iw];
+                    }
+                }
+            }
         }
     }
 
-    this->loc->cal_dk_k(GlobalC::GridT);
+    this->loc->cal_dk_k(GlobalC::GridT, this->wg);
     for (int is = 0; is < GlobalV::NSPIN; is++)
     {
         ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], this->charge->nrxx); // mohan 2009-11-10
@@ -119,9 +92,9 @@ void ElecStateLCAO::psiToRho(const psi::Psi<double>& psi)
     {
         ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d");
 
-        psi::Psi<double> dm_gamma_2d(psi.get_nk(), psi.get_nbasis(), psi.get_nbasis());
+        //psi::Psi<double> dm_gamma_2d;
         // caution:wfc and dm
-        this->cal_dm(this->wg, psi, dm_gamma_2d);
+        cal_dm(this->loc->ParaV, this->wg, psi, this->loc->dm_gamma);
 
         ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d");
 
@@ -132,9 +105,10 @@ void ElecStateLCAO::psiToRho(const psi::Psi<double>& psi)
             {
                 psi.fix_k(ik);
                 double** wfc_grid = nullptr; // output but not do "2d-to-grid" conversion
-                this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi.get_pointer(), wfc_grid);
+                this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi.get_pointer(), wfc_grid, this->ekb, this->wg);
             }
-            this->loc->dm2dToGrid(dm_gamma_2d, this->loc->DM[ik]); // transform dm_gamma[is].c to this->loc->DM[is]
+            //this->loc->dm2dToGrid(this->loc->dm_gamma[ik], this->loc->DM[ik]); // transform dm_gamma[is].c to this->loc->DM[is]
+            this->loc->cal_dk_gamma_from_2D_pub();
         }
     }
 

source/module_elecstate/elecstate_lcao.h

@@ -48,8 +48,6 @@ class ElecStateLCAO : public ElecState
     // calcualte rho for each k
     // void rhoBandK(const psi::Psi<std::complex<double>>& psi);
 
-    template <typename T> void cal_dm(const ModuleBase::matrix& wg, const psi::Psi<T>& wfc, psi::Psi<T>& dm);
-
     Local_Orbital_Charge* loc = nullptr;
     LCAO_Hamilt* uhm = nullptr;
     Local_Orbital_wfc* lowf = nullptr;

source/module_elecstate/elecstate_pw.cpp

@@ -3,12 +3,15 @@
 #include "module_base/constants.h"
 #include "src_parallel/parallel_reduce.h"
 #include "src_pw/global.h"
+#include "module_base/timer.h"
 
 namespace elecstate
 {
 
 void ElecStatePW::psiToRho(const psi::Psi<std::complex<double>>& psi)
 {
+    ModuleBase::TITLE("ElecStatePW", "psiToRho");
+    ModuleBase::timer::tick("ElecStatePW", "psiToRho");
     this->calculate_weights();
 
     this->calEBand();
@@ -28,6 +31,7 @@ void ElecStatePW::psiToRho(const psi::Psi<std::complex<double>>& psi)
         this->updateRhoK(psi);
     }
     this->parallelK();
+    ModuleBase::timer::tick("ElecStatePW", "psiToRho");
 }
 
 void ElecStatePW::updateRhoK(const psi::Psi<std::complex<double>>& psi)