diff --git a/source/module_hsolver/diago_pexsi.cpp b/source/module_hsolver/diago_pexsi.cpp index 4076922626a..f7bd881bed4 100644 --- a/source/module_hsolver/diago_pexsi.cpp +++ b/source/module_hsolver/diago_pexsi.cpp @@ -13,6 +13,9 @@ typedef hamilt::MatrixBlock> matcd; namespace hsolver { +template +std::vector DiagoPexsi::mu_buffer; + template DiagoPexsi::DiagoPexsi(const Parallel_Orbitals* ParaV_in) { @@ -76,7 +79,7 @@ void DiagoPexsi::diag(hamilt::Hamilt* phm_in, psi::Psi& this->totalFreeEnergy = this->ps->get_totalFreeEnergy(); this->totalEnergyH = this->ps->get_totalEnergyH(); this->totalEnergyS = this->ps->get_totalEnergyS(); - this->mu_buffer[ik] = this->ps->get_mu(); + mu_buffer[ik] = this->ps->get_mu(); } template <> diff --git a/source/module_hsolver/diago_pexsi.h b/source/module_hsolver/diago_pexsi.h index 9802273e4b5..657471a9e50 100644 --- a/source/module_hsolver/diago_pexsi.h +++ b/source/module_hsolver/diago_pexsi.h @@ -16,7 +16,7 @@ class DiagoPexsi : public DiagH { private: using Real = typename GetTypeReal::type; - std::vector mu_buffer; + static std::vector mu_buffer; public: DiagoPexsi(const Parallel_Orbitals* ParaV_in); diff --git a/source/module_hsolver/hsolver_lcao.cpp b/source/module_hsolver/hsolver_lcao.cpp index adb193cf336..21f94ecf95e 100644 --- a/source/module_hsolver/hsolver_lcao.cpp +++ b/source/module_hsolver/hsolver_lcao.cpp @@ -52,24 +52,23 @@ void HSolverLCAO::solve(hamilt::Hamilt* pHamilt, // select the method of diagonalization this->method = method_in; -#ifdef __PEXSI +#ifdef __PEXSI // other purification methods should follow this routine if (this->method == "pexsi") { - if (this->pdiagh != nullptr) - { - if (this->pdiagh->method != this->method) - { - delete[] this->pdiagh; - this->pdiagh = nullptr; - } - auto tem = dynamic_cast*>(this->pdiagh); - } - if (this->pdiagh == nullptr) { - DiagoPexsi* tem = new DiagoPexsi(this->ParaV); - this->pdiagh = tem; - // this->pdiagh = dynamic_cast*>(tem); - this->pdiagh->method = this->method; + DiagoPexsi pe(ParaV); + for (int ik = 0; ik < psi.get_nk(); ++ik) { + /// update H(k) for each k point + pHamilt->updateHk(ik); + psi.fix_k(ik); + // solve eigenvector and eigenvalue for H(k) + pe.diag(pHamilt, psi, nullptr); } + auto _pes = dynamic_cast*>(pes); + pes->f_en.eband = pe.totalFreeEnergy; + // maybe eferm could be dealt with in the future + _pes->dmToRho(pe.DM, pe.EDM); + ModuleBase::timer::tick("HSolverLCAO", "solve"); + return; } #endif @@ -100,6 +99,7 @@ void HSolverLCAO::solve(hamilt::Hamilt* pHamilt, pHamilt->updateHk(ik); psi.fix_k(ik); + std::cout << pes->nelec_spin[0] << " " << pes->nelec_spin[1] << std::endl; // solve eigenvector and eigenvalue for H(k) this->hamiltSolvePsiK(pHamilt, psi, &(pes->ekb(ik, 0))); @@ -126,22 +126,7 @@ void HSolverLCAO::solve(hamilt::Hamilt* pHamilt, // calculate charge by psi // called in scf calculation -#ifdef __PEXSI - if (this->method == "pexsi") - { - DiagoPexsi tem = dynamic_cast*>(this->pdiagh); - if (tem == nullptr) - ModuleBase::WARNING_QUIT("HSolverLCAO", "pexsi need debug!"); - elecstate::ElecStateLCAO* _pes - = dynamic_cast*>(pes); - pes->f_en.eband = tem->totalFreeEnergy; - // maybe eferm could be dealt with in the future - _pes->dmToRho(tem->DM, tem->EDM); - } else -#endif - { - pes->psiToRho(psi); - } + pes->psiToRho(psi); ModuleBase::timer::tick("HSolverLCAO", "solve"); } @@ -188,12 +173,6 @@ void HSolverLCAO::hamiltSolvePsiK(hamilt::Hamilt* hm, psi::Psi& la.diag(hm, psi, eigenvalue); #else ModuleBase::WARNING_QUIT("HSolverLCAO::solve", "This method of DiagH is not supported!"); -#endif -#ifdef __PEXSI - else if (this->method == "pexsi") - { - this->pdiagh->diag(hm, psi, eigenvalue); - } #endif } else diff --git a/source/module_hsolver/module_pexsi/simple_pexsi.cpp b/source/module_hsolver/module_pexsi/simple_pexsi.cpp index 075c8031829..d2961c40bd8 100644 --- a/source/module_hsolver/module_pexsi/simple_pexsi.cpp +++ b/source/module_hsolver/module_pexsi/simple_pexsi.cpp @@ -110,11 +110,11 @@ int loadPEXSIOption(MPI_Comm comm, int_para[11] = pexsi::PEXSI_Solver::pexsi_symm; int_para[12] = pexsi::PEXSI_Solver::pexsi_trans; int_para[13] = pexsi::PEXSI_Solver::pexsi_method; - int_para[14] = 2; + int_para[14] = 1; int_para[15] = 0; int_para[16] = pexsi::PEXSI_Solver::pexsi_nproc_pole; - double_para[0] = GlobalV::NSPIN; // pexsi::PEXSI_Solver::pexsi_spin; + double_para[0] = 2;//GlobalV::NSPIN; // pexsi::PEXSI_Solver::pexsi_spin; double_para[1] = pexsi::PEXSI_Solver::pexsi_temp; double_para[2] = pexsi::PEXSI_Solver::pexsi_gap; double_para[3] = pexsi::PEXSI_Solver::pexsi_delta_e; @@ -298,15 +298,15 @@ int simplePEXSI(MPI_Comm comm_PEXSI, int numTotalInertiaIter; // Number of total inertia[out] // LiuXh modify 2021-04-29, add DONE(ofs_running,"xx") for test ModuleBase::timer::tick("Diago_LCAO_Matrix", "PEXSIDFT"); - PPEXSIDFTDriver(plan, // PEXSI plan[in] - options, // PEXSI Options[in] + PPEXSIDFTDriver2(plan, // PEXSI plan[in] + &options, // PEXSI Options[in] numElectronExact, // exact electron number[in] &mu, // chemical potential[out] &nelec, // number of electrons[out] - &muMinInertia, // Lower bound for mu after the last inertia[out] - &muMaxInertia, // Upper bound for mu after the last inertia[out] + // &muMinInertia, // Lower bound for mu after the last inertia[out] + // &muMaxInertia, // Upper bound for mu after the last inertia[out] &numTotalInertiaIter, // Number of total inertia[out] - &numTotalPEXSIIter, // number of total pexsi evaluation procedure[out] + // &numTotalPEXSIIter, // number of total pexsi evaluation procedure[out] &info); // 0: successful; otherwise: unsuccessful // LiuXh modify 2021-04-29, add DONE(ofs_running,"xx") for test ModuleBase::timer::tick("Diago_LCAO_Matrix", "PEXSIDFT");