From 87816025a5c5016da60571d522af55a17795b70d Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Thu, 12 Sep 2024 02:51:56 +0800 Subject: [PATCH 1/3] remove GlobalC::ORB in module_ri --- source/module_esolver/esolver_ks_lcao.cpp | 9 ++++--- source/module_esolver/lcao_before_scf.cpp | 4 +-- source/module_lr/esolver_lrtd_lcao.cpp | 4 +-- source/module_ri/ABFs_Construct-PCA.cpp | 3 ++- source/module_ri/ABFs_Construct-PCA.h | 5 ++-- source/module_ri/Exx_LRI.h | 5 ++-- source/module_ri/Exx_LRI.hpp | 14 +++++----- source/module_ri/Exx_LRI_interface.h | 4 +-- source/module_ri/Exx_LRI_interface.hpp | 6 ++--- source/module_ri/LRI_CV.h | 2 ++ source/module_ri/LRI_CV.hpp | 10 ++++--- source/module_ri/LRI_CV_Tools.h | 4 +-- source/module_ri/LRI_CV_Tools.hpp | 6 ++--- source/module_ri/Matrix_Orbs11.cpp | 14 +++++----- source/module_ri/Matrix_Orbs11.h | 1 + source/module_ri/Matrix_Orbs21.cpp | 14 +++++----- source/module_ri/Matrix_Orbs21.h | 1 + source/module_ri/Matrix_Orbs22.cpp | 14 +++++----- source/module_ri/Matrix_Orbs22.h | 1 + source/module_ri/RPA_LRI.h | 9 ++++--- source/module_ri/RPA_LRI.hpp | 15 ++++++----- source/module_ri/exx_abfs-construct_orbs.cpp | 8 +++--- source/module_ri/exx_abfs-construct_orbs.h | 2 ++ source/module_ri/exx_abfs-jle.cpp | 28 ++++++++++---------- source/module_ri/exx_abfs-jle.h | 4 +-- source/module_ri/exx_opt_orb-print.cpp | 11 ++++---- source/module_ri/exx_opt_orb.cpp | 27 +++++++++++-------- source/module_ri/exx_opt_orb.h | 6 +++-- 28 files changed, 131 insertions(+), 100 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 1c79f13405..d12bf618be 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -204,8 +204,8 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell { XC_Functional::set_xc_first_loop(ucell); // initialize 2-center radial tables for EXX-LRI - if (GlobalC::exx_info.info_ri.real_number) { this->exx_lri_double->init(MPI_COMM_WORLD, this->kv); } - else { this->exx_lri_complex->init(MPI_COMM_WORLD, this->kv); } + if (GlobalC::exx_info.info_ri.real_number) { this->exx_lri_double->init(MPI_COMM_WORLD, this->kv, orb_); } + else { this->exx_lri_complex->init(MPI_COMM_WORLD, this->kv, orb_); } } } #endif @@ -1175,8 +1175,9 @@ void ESolver_KS_LCAO::after_scf(const int istep) RPA_LRI rpa_lri_double(GlobalC::exx_info.info_ri); rpa_lri_double.cal_postSCF_exx(*dynamic_cast*>(this->pelec)->get_DM(), MPI_COMM_WORLD, - this->kv); - rpa_lri_double.init(MPI_COMM_WORLD, this->kv); + this->kv, + orb_); + rpa_lri_double.init(MPI_COMM_WORLD, this->kv, orb_.cutoffs()); rpa_lri_double.out_for_RPA(this->pv, *(this->psi), this->pelec); } #endif diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 9df28c64ed..e23d848d78 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -210,11 +210,11 @@ void ESolver_KS_LCAO::before_scf(const int istep) #ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf if (GlobalC::exx_info.info_ri.real_number) { - this->exd->exx_beforescf(this->kv, *this->p_chgmix, GlobalC::ucell, this->pv); + this->exd->exx_beforescf(this->kv, *this->p_chgmix, GlobalC::ucell, this->pv, orb_); } else { - this->exc->exx_beforescf(this->kv, *this->p_chgmix, GlobalC::ucell, this->pv); + this->exc->exx_beforescf(this->kv, *this->p_chgmix, GlobalC::ucell, this->pv, orb_); } #endif // __EXX diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 42333f52f6..8df436d3bf 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -200,7 +200,7 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol } else // construct C, V from scratch { this->exx_lri = std::make_shared>(exx_info.info_ri); - this->exx_lri->init(MPI_COMM_WORLD, this->kv); + this->exx_lri->init(MPI_COMM_WORLD, this->kv, ks_sol.orb_); this->exx_lri->cal_exx_ions(); } } @@ -365,7 +365,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu if ((xc_kernel == "hf" || xc_kernel == "hse") && this->input.lr_solver != "spectrum") { this->exx_lri = std::make_shared>(exx_info.info_ri); - this->exx_lri->init(MPI_COMM_WORLD, this->kv); + this->exx_lri->init(MPI_COMM_WORLD, this->kv, orb); this->exx_lri->cal_exx_ions(); } // else diff --git a/source/module_ri/ABFs_Construct-PCA.cpp b/source/module_ri/ABFs_Construct-PCA.cpp index 4c1471560c..cf4c7b8fdf 100644 --- a/source/module_ri/ABFs_Construct-PCA.cpp +++ b/source/module_ri/ABFs_Construct-PCA.cpp @@ -67,6 +67,7 @@ namespace PCA } std::vector, RI::Tensor>>> cal_PCA( + const LCAO_Orbitals& orb, const std::vector>> &lcaos, const std::vector>> &abfs, const double kmesh_times ) @@ -89,7 +90,7 @@ namespace PCA GlobalC::exx_info.info_ri.abfs_Lmax = std::max( GlobalC::exx_info.info_ri.abfs_Lmax, static_cast(abfs[T].size())-1 ); Matrix_Orbs21 m_abfslcaos_lcaos; - m_abfslcaos_lcaos.init( 1, kmesh_times, 1 ); + m_abfslcaos_lcaos.init( 1, orb, kmesh_times, 1 ); m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos ); std::map>> delta_R; diff --git a/source/module_ri/ABFs_Construct-PCA.h b/source/module_ri/ABFs_Construct-PCA.h index b172b578e3..bbd76cb61e 100644 --- a/source/module_ri/ABFs_Construct-PCA.h +++ b/source/module_ri/ABFs_Construct-PCA.h @@ -1,7 +1,7 @@ #ifndef ABFS_CONSTRUCT_PCA_H #define ABFS_CONSTRUCT_PCA_H -#include "../module_basis/module_ao/ORB_atomic_lm.h" +#include "../module_basis/module_ao/ORB_read.h" #include #include @@ -15,10 +15,11 @@ namespace ABFs_Construct namespace PCA { extern std::vector,RI::Tensor>>> cal_PCA( + const LCAO_Orbitals &orb, const std::vector>> &lcaos, const std::vector>> &abfs, // abfs must be orthonormal const double kmesh_times ); } } -#endif // ABFS_CONSTRUCT_PCA_H \ No newline at end of file +#endif // ABFS_CONSTRUCT_PCA_H diff --git a/source/module_ri/Exx_LRI.h b/source/module_ri/Exx_LRI.h index 0085569a64..30effc2447 100644 --- a/source/module_ri/Exx_LRI.h +++ b/source/module_ri/Exx_LRI.h @@ -54,7 +54,7 @@ class Exx_LRI Exx_LRI operator=(Exx_LRI&&); - void init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in); + void init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, const LCAO_Orbitals& orb); void cal_exx_force(); void cal_exx_stress(); std::vector> get_abfs_nchis() const; @@ -69,6 +69,7 @@ class Exx_LRI const Exx_Info::Exx_Info_RI &info; MPI_Comm mpi_comm; const K_Vectors *p_kv = nullptr; + std::vector orb_cutoff_; std::vector>> lcaos; std::vector>> abfs; @@ -96,4 +97,4 @@ class Exx_LRI #include "Exx_LRI.hpp" -#endif \ No newline at end of file +#endif diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index cbe0f91027..9e159fc844 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -26,7 +26,7 @@ #include template -void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in) +void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, const LCAO_Orbitals& orb) { ModuleBase::TITLE("Exx_LRI","init"); ModuleBase::timer::tick("Exx_LRI", "init"); @@ -49,19 +49,20 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in) this->mpi_comm = mpi_comm_in; this->p_kv = &kv_in; + this->orb_cutoff_ = orb.cutoffs(); - this->lcaos = Exx_Abfs::Construct_Orbs::change_orbs( GlobalC::ORB, this->info.kmesh_times ); + this->lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->info.kmesh_times ); // #ifdef __MPI // Exx_Abfs::Util::bcast( this->info.files_abfs, 0, this->mpi_comm ); // #endif const std::vector>> - abfs_same_atom = Exx_Abfs::Construct_Orbs::abfs_same_atom( this->lcaos, this->info.kmesh_times, this->info.pca_threshold ); + abfs_same_atom = Exx_Abfs::Construct_Orbs::abfs_same_atom( orb, this->lcaos, this->info.kmesh_times, this->info.pca_threshold ); if(this->info.files_abfs.empty()) { this->abfs = abfs_same_atom; } else { - this->abfs = Exx_Abfs::IO::construct_abfs( abfs_same_atom, GlobalC::ORB, this->info.files_abfs, this->info.kmesh_times ); + this->abfs = Exx_Abfs::IO::construct_abfs( abfs_same_atom, orb, this->info.files_abfs, this->info.kmesh_times ); } Exx_Abfs::Construct_Orbs::print_orbs_size(this->abfs, GlobalV::ofs_running); @@ -92,6 +93,7 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in) } this->cv.set_orbitals( + orb, this->lcaos, this->abfs, this->abfs_ccp, this->info.kmesh_times, this->info.ccp_rmesh_times ); @@ -126,7 +128,7 @@ void Exx_LRI::cal_exx_ions() this->exx_lri.set_parallel(this->mpi_comm, atoms_pos, latvec, period); // std::max(3) for gamma_only, list_A2 should contain cell {-1,0,1}. In the future distribute will be neighbour. - const std::array period_Vs = LRI_CV_Tools::cal_latvec_range(1+this->info.ccp_rmesh_times); + const std::array period_Vs = LRI_CV_Tools::cal_latvec_range(1+this->info.ccp_rmesh_times, orb_cutoff_); const std::pair, std::vector>>>> list_As_Vs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false); @@ -152,7 +154,7 @@ void Exx_LRI::cal_exx_ions() } } - const std::array period_Cs = LRI_CV_Tools::cal_latvec_range(2); + const std::array period_Cs = LRI_CV_Tools::cal_latvec_range(2, orb_cutoff_); const std::pair, std::vector>>>> list_As_Cs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false); diff --git a/source/module_ri/Exx_LRI_interface.h b/source/module_ri/Exx_LRI_interface.h index 33f46cf814..d1b64c56df 100644 --- a/source/module_ri/Exx_LRI_interface.h +++ b/source/module_ri/Exx_LRI_interface.h @@ -49,7 +49,7 @@ class Exx_LRI_Interface // Processes in ESolver_KS_LCAO /// @brief in beforescf: set xc type, opt_orb, do DM mixing - void exx_beforescf(const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const Parallel_2D& pv); + void exx_beforescf(const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const Parallel_2D& pv, const LCAO_Orbitals& orb); /// @brief in eachiterinit: do DM mixing and calculate Hexx when entering 2nd SCF void exx_eachiterinit(const elecstate::DensityMatrix& dm/**< double should be Tdata if complex-PBE-DM is supported*/, @@ -79,4 +79,4 @@ class Exx_LRI_Interface #include "Exx_LRI_interface.hpp" -#endif \ No newline at end of file +#endif diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index e4d8c1d8d8..fcf641ca40 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -37,7 +37,7 @@ void Exx_LRI_Interface::read_Hexxs_cereal(const std::string& file_name } template -void Exx_LRI_Interface::exx_beforescf(const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const Parallel_2D& pv) +void Exx_LRI_Interface::exx_beforescf(const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const Parallel_2D& pv, const LCAO_Orbitals& orb) { #ifdef __MPI if (GlobalC::exx_info.info_global.cal_exx) @@ -72,7 +72,7 @@ void Exx_LRI_Interface::exx_beforescf(const K_Vectors& kv, const Charg { //program should be stopped after this judgement Exx_Opt_Orb exx_opt_orb; - exx_opt_orb.generate_matrix(kv); + exx_opt_orb.generate_matrix(kv, orb); ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf"); return; } @@ -264,4 +264,4 @@ bool Exx_LRI_Interface::exx_after_converge( restart_reset(); return true; } -#endif \ No newline at end of file +#endif diff --git a/source/module_ri/LRI_CV.h b/source/module_ri/LRI_CV.h index 26250d9d56..252259c119 100644 --- a/source/module_ri/LRI_CV.h +++ b/source/module_ri/LRI_CV.h @@ -34,6 +34,7 @@ class LRI_CV ~LRI_CV(); void set_orbitals( + const LCAO_Orbitals& orb, const std::vector>> &lcaos_in, const std::vector>> &abfs_in, const std::vector>> &abfs_ccp_in, @@ -59,6 +60,7 @@ class LRI_CV size_t get_index_abfs_size(const size_t &iat){return this->index_abfs[iat].count_size; } private: + std::vector orb_cutoff_; std::vector>> lcaos; std::vector>> abfs; std::vector>> abfs_ccp; diff --git a/source/module_ri/LRI_CV.hpp b/source/module_ri/LRI_CV.hpp index 5f312b0b2f..8640088130 100644 --- a/source/module_ri/LRI_CV.hpp +++ b/source/module_ri/LRI_CV.hpp @@ -37,6 +37,7 @@ LRI_CV::~LRI_CV() template void LRI_CV::set_orbitals( + const LCAO_Orbitals& orb, const std::vector>> &lcaos_in, const std::vector>> &abfs_in, const std::vector>> &abfs_ccp_in, @@ -46,6 +47,7 @@ void LRI_CV::set_orbitals( ModuleBase::TITLE("LRI_CV", "set_orbitals"); ModuleBase::timer::tick("LRI_CV", "set_orbitals"); + this->orb_cutoff_ = orb.cutoffs(); this->lcaos = lcaos_in; this->abfs = abfs_in; this->abfs_ccp = abfs_ccp_in; @@ -59,11 +61,11 @@ void LRI_CV::set_orbitals( range_abfs = Exx_Abfs::Abfs_Index::construct_range( abfs ); this->index_abfs = ModuleBase::Element_Basis_Index::construct_index( range_abfs ); - this->m_abfs_abfs.init( 2, kmesh_times, (1+this->ccp_rmesh_times)/2.0 ); + this->m_abfs_abfs.init( 2, orb, kmesh_times, (1+this->ccp_rmesh_times)/2.0 ); this->m_abfs_abfs.init_radial( this->abfs_ccp, this->abfs ); this->m_abfs_abfs.init_radial_table(); - this->m_abfslcaos_lcaos.init( 1, kmesh_times, 1 ); + this->m_abfslcaos_lcaos.init( 1, orb, kmesh_times, 1 ); this->m_abfslcaos_lcaos.init_radial( this->abfs_ccp, this->lcaos, this->lcaos ); this->m_abfslcaos_lcaos.init_radial_table(); @@ -101,8 +103,8 @@ auto LRI_CV::cal_datas( const ModuleBase::Vector3 tau0 = GlobalC::ucell.atoms[it0].tau[ia0]; const ModuleBase::Vector3 tau1 = GlobalC::ucell.atoms[it1].tau[ia1]; const double Rcut = std::min( - GlobalC::ORB.Phi[it0].getRcut() * rmesh_times + GlobalC::ORB.Phi[it1].getRcut(), - GlobalC::ORB.Phi[it1].getRcut() * rmesh_times + GlobalC::ORB.Phi[it0].getRcut()); + orb_cutoff_[it0] * rmesh_times + orb_cutoff_[it1], + orb_cutoff_[it1] * rmesh_times + orb_cutoff_[it0]); const Abfs::Vector3_Order R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*GlobalC::ucell.latvec); if( R_delta.norm()*GlobalC::ucell.lat0 < Rcut ) { diff --git a/source/module_ri/LRI_CV_Tools.h b/source/module_ri/LRI_CV_Tools.h index 57cfa3092b..f7ccf1f8cf 100644 --- a/source/module_ri/LRI_CV_Tools.h +++ b/source/module_ri/LRI_CV_Tools.h @@ -78,7 +78,7 @@ namespace LRI_CV_Tools std::map>> && ds_in); template - extern std::array cal_latvec_range(const double &rcut_times); + extern std::array cal_latvec_range(const double &rcut_times, const std::vector& orb_cutoff); template extern std::map,RI::Tensor>>> @@ -97,4 +97,4 @@ namespace LRI_CV_Tools #include "LRI_CV_Tools.hpp" -#endif \ No newline at end of file +#endif diff --git a/source/module_ri/LRI_CV_Tools.hpp b/source/module_ri/LRI_CV_Tools.hpp index 29d260979b..3e720ba9f4 100644 --- a/source/module_ri/LRI_CV_Tools.hpp +++ b/source/module_ri/LRI_CV_Tools.hpp @@ -245,11 +245,11 @@ LRI_CV_Tools::change_order(std::map>> template std::array -LRI_CV_Tools::cal_latvec_range(const double &rcut_times) +LRI_CV_Tools::cal_latvec_range(const double &rcut_times, const std::vector& orb_cutoff) { double Rcut_max = 0; for(int T=0; T proj = ModuleBase::Mathzone::latvec_projection( std::array,3>{GlobalC::ucell.a1, GlobalC::ucell.a2, GlobalC::ucell.a3}); const ModuleBase::Vector3 latvec_times = Rcut_max * rcut_times / (proj * GlobalC::ucell.lat0); @@ -356,4 +356,4 @@ LRI_CV_Tools::cal_dMRs( } return dMRs; } -#endif \ No newline at end of file +#endif diff --git a/source/module_ri/Matrix_Orbs11.cpp b/source/module_ri/Matrix_Orbs11.cpp index d42dc53378..17fdb35238 100644 --- a/source/module_ri/Matrix_Orbs11.cpp +++ b/source/module_ri/Matrix_Orbs11.cpp @@ -9,24 +9,24 @@ #include "module_base/tool_title.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -void Matrix_Orbs11::init(const int mode, const double kmesh_times, const double rmesh_times) +void Matrix_Orbs11::init(const int mode, const LCAO_Orbitals& orb, const double kmesh_times, const double rmesh_times) { ModuleBase::TITLE("Matrix_Orbs11", "init"); ModuleBase::timer::tick("Matrix_Orbs11", "init"); int Lmax_used, Lmax; - const int ntype = GlobalC::ORB.get_ntype(); + const int ntype = orb.get_ntype(); int lmax_orb = -1, lmax_beta = -1; for (int it = 0; it < ntype; it++) { - lmax_orb = std::max(lmax_orb, GlobalC::ORB.Phi[it].getLmax()); + lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); lmax_beta = std::max(lmax_beta, GlobalC::ucell.infoNL.Beta[it].getLmax()); } - const double dr = GlobalC::ORB.get_dR(); - const double dk = GlobalC::ORB.get_dk(); - const int kmesh = GlobalC::ORB.get_kmesh() * kmesh_times + 1; - int Rmesh = static_cast(GlobalC::ORB.get_Rmax() * rmesh_times / dr) + 4; + const double dr = orb.get_dR(); + const double dk = orb.get_dk(); + const int kmesh = orb.get_kmesh() * kmesh_times + 1; + int Rmesh = static_cast(orb.get_Rmax() * rmesh_times / dr) + 4; Rmesh += 1 - Rmesh % 2; Center2_Orb::init_Table_Spherical_Bessel(2, diff --git a/source/module_ri/Matrix_Orbs11.h b/source/module_ri/Matrix_Orbs11.h index 4ea51f3abd..8861198029 100644 --- a/source/module_ri/Matrix_Orbs11.h +++ b/source/module_ri/Matrix_Orbs11.h @@ -25,6 +25,7 @@ class Matrix_Orbs11 // 1: // 2: void init(const int mode, + const LCAO_Orbitals& orb, const double kmesh_times, // extend Kcut, keep dK const double rmesh_times); // extend Rcut, keep dR diff --git a/source/module_ri/Matrix_Orbs21.cpp b/source/module_ri/Matrix_Orbs21.cpp index a6dd81d45c..a772598ce8 100644 --- a/source/module_ri/Matrix_Orbs21.cpp +++ b/source/module_ri/Matrix_Orbs21.cpp @@ -9,23 +9,23 @@ #include "module_base/tool_title.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -void Matrix_Orbs21::init(const int mode, const double kmesh_times, const double rmesh_times) +void Matrix_Orbs21::init(const int mode, const LCAO_Orbitals& orb, const double kmesh_times, const double rmesh_times) { ModuleBase::TITLE("Matrix_Orbs21", "init"); ModuleBase::timer::tick("Matrix_Orbs21", "init"); int Lmax_used, Lmax; - const int ntype = GlobalC::ORB.get_ntype(); + const int ntype = orb.get_ntype(); int lmax_orb = -1, lmax_beta = -1; for (int it = 0; it < ntype; it++) { - lmax_orb = std::max(lmax_orb, GlobalC::ORB.Phi[it].getLmax()); + lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); lmax_beta = std::max(lmax_beta, GlobalC::ucell.infoNL.Beta[it].getLmax()); } - const double dr = GlobalC::ORB.get_dR(); - const double dk = GlobalC::ORB.get_dk(); - const int kmesh = GlobalC::ORB.get_kmesh() * kmesh_times + 1; - int Rmesh = static_cast(GlobalC::ORB.get_Rmax() * rmesh_times / dr) + 4; + const double dr = orb.get_dR(); + const double dk = orb.get_dk(); + const int kmesh = orb.get_kmesh() * kmesh_times + 1; + int Rmesh = static_cast(orb.get_Rmax() * rmesh_times / dr) + 4; Rmesh += 1 - Rmesh % 2; Center2_Orb::init_Table_Spherical_Bessel(3, diff --git a/source/module_ri/Matrix_Orbs21.h b/source/module_ri/Matrix_Orbs21.h index 780a745826..680b11b0da 100644 --- a/source/module_ri/Matrix_Orbs21.h +++ b/source/module_ri/Matrix_Orbs21.h @@ -23,6 +23,7 @@ class Matrix_Orbs21 // mode: // 1: void init(const int mode, + const LCAO_Orbitals& orb, const double kmesh_times, // extend Kcut, keep dK const double rmesh_times); // extend Rcut, keep dR diff --git a/source/module_ri/Matrix_Orbs22.cpp b/source/module_ri/Matrix_Orbs22.cpp index c882c8cb11..027f252614 100644 --- a/source/module_ri/Matrix_Orbs22.cpp +++ b/source/module_ri/Matrix_Orbs22.cpp @@ -9,23 +9,23 @@ #include "module_base/tool_title.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -void Matrix_Orbs22::init(const int mode, const double kmesh_times, const double rmesh_times) +void Matrix_Orbs22::init(const int mode, const LCAO_Orbitals& orb, const double kmesh_times, const double rmesh_times) { ModuleBase::TITLE("Matrix_Orbs22", "init"); ModuleBase::timer::tick("Matrix_Orbs22", "init"); int Lmax_used, Lmax; - const int ntype = GlobalC::ORB.get_ntype(); + const int ntype = orb.get_ntype(); int lmax_orb = -1, lmax_beta = -1; for (int it = 0; it < ntype; it++) { - lmax_orb = std::max(lmax_orb, GlobalC::ORB.Phi[it].getLmax()); + lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); lmax_beta = std::max(lmax_beta, GlobalC::ucell.infoNL.Beta[it].getLmax()); } - const double dr = GlobalC::ORB.get_dR(); - const double dk = GlobalC::ORB.get_dk(); - const int kmesh = GlobalC::ORB.get_kmesh() * kmesh_times + 1; - int Rmesh = static_cast(GlobalC::ORB.get_Rmax() * rmesh_times / dr) + 4; + const double dr = orb.get_dR(); + const double dk = orb.get_dk(); + const int kmesh = orb.get_kmesh() * kmesh_times + 1; + int Rmesh = static_cast(orb.get_Rmax() * rmesh_times / dr) + 4; Rmesh += 1 - Rmesh % 2; Center2_Orb::init_Table_Spherical_Bessel(4, diff --git a/source/module_ri/Matrix_Orbs22.h b/source/module_ri/Matrix_Orbs22.h index 382b9bbe69..0c7f08be3f 100644 --- a/source/module_ri/Matrix_Orbs22.h +++ b/source/module_ri/Matrix_Orbs22.h @@ -23,6 +23,7 @@ class Matrix_Orbs22 // mode: // 1: void init(const int mode, + const LCAO_Orbitals& orb, const double kmesh_times, // extend Kcut, keep dK const double rmesh_times); // extend Rcut, keep dR diff --git a/source/module_ri/RPA_LRI.h b/source/module_ri/RPA_LRI.h index 70222cf8c2..b169b0884b 100644 --- a/source/module_ri/RPA_LRI.h +++ b/source/module_ri/RPA_LRI.h @@ -37,11 +37,12 @@ template class RPA_LRI { } ~RPA_LRI(){}; - void init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in); + void init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, const std::vector& orb_cutoff); void cal_rpa_cv(); void cal_postSCF_exx(const elecstate::DensityMatrix& dm, const MPI_Comm& mpi_comm_in, - const K_Vectors& kv); + const K_Vectors& kv, + const LCAO_Orbitals& orb); void out_for_RPA(const Parallel_Orbitals& parav, const psi::Psi& psi, const elecstate::ElecState* pelec); @@ -62,6 +63,8 @@ template class RPA_LRI const Exx_Info::Exx_Info_RI &info; const K_Vectors *p_kv=nullptr; MPI_Comm mpi_comm; + std::vector orb_cutoff_; + std::vector>> lcaos; std::vector>> abfs; std::vector>> abfs_ccp; @@ -77,4 +80,4 @@ template class RPA_LRI Exx_LRI exx_lri_rpa(GlobalC::exx_info.info_ri); #include "RPA_LRI.hpp" -#endif \ No newline at end of file +#endif diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp index c20fa860f5..08b287ae2d 100644 --- a/source/module_ri/RPA_LRI.hpp +++ b/source/module_ri/RPA_LRI.hpp @@ -8,17 +8,19 @@ #include #include #include +#include #include "module_ri/module_exx_symmetry/symmetry_rotation.h" #include "RPA_LRI.h" #include "module_parameter/parameter.h" template -void RPA_LRI::init(const MPI_Comm& mpi_comm_in, const K_Vectors& kv_in) +void RPA_LRI::init(const MPI_Comm& mpi_comm_in, const K_Vectors& kv_in, const std::vector& orb_cutoff) { ModuleBase::TITLE("RPA_LRI", "init"); ModuleBase::timer::tick("RPA_LRI", "init"); this->mpi_comm = mpi_comm_in; + this->orb_cutoff_ = orb_cutoff; this->lcaos = exx_lri_rpa.lcaos; this->abfs = exx_lri_rpa.abfs; this->abfs_ccp = exx_lri_rpa.abfs_ccp; @@ -38,7 +40,7 @@ void RPA_LRI::cal_rpa_cv() } const std::array period = {p_kv->nmp[0], p_kv->nmp[1], p_kv->nmp[2]}; - const std::array period_Vs = LRI_CV_Tools::cal_latvec_range(1 + this->info.ccp_rmesh_times); + const std::array period_Vs = LRI_CV_Tools::cal_latvec_range(1 + this->info.ccp_rmesh_times, orb_cutoff_); const std::pair, std::vector>>>> list_As_Vs = RI::Distribute_Equally::distribute_atoms(this->mpi_comm, atoms, period_Vs, 2, false); @@ -49,7 +51,7 @@ void RPA_LRI::cal_rpa_cv() }); this->Vs_period = RI::RI_Tools::cal_period(Vs, period); - const std::array period_Cs = LRI_CV_Tools::cal_latvec_range(2); + const std::array period_Cs = LRI_CV_Tools::cal_latvec_range(2, orb_cutoff_); const std::pair, std::vector>>>> list_As_Cs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false); @@ -71,7 +73,8 @@ void RPA_LRI::cal_rpa_cv() template void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix& dm, const MPI_Comm& mpi_comm_in, - const K_Vectors& kv) + const K_Vectors& kv, + const LCAO_Orbitals& orb) { Mix_DMk_2D mix_DMk_2D; bool exx_spacegroup_symmetry = (GlobalV::NSPIN < 4 && ModuleSymmetry::Symmetry::symm_flag == 1); @@ -102,7 +105,7 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix GlobalC::exx_info.info_global.hybrid_alpha = 1; GlobalC::exx_info.info_ri.ccp_rmesh_times = PARAM.inp.rpa_ccp_rmesh_times; - exx_lri_rpa.init(mpi_comm_in, kv); + exx_lri_rpa.init(mpi_comm_in, kv, orb); exx_lri_rpa.cal_exx_ions(); if (exx_spacegroup_symmetry) @@ -456,4 +459,4 @@ void RPA_LRI::out_coulomb_k() // // rpa_lri.set_Cs(Cs); // } -#endif \ No newline at end of file +#endif diff --git a/source/module_ri/exx_abfs-construct_orbs.cpp b/source/module_ri/exx_abfs-construct_orbs.cpp index be5c9eb085..55eed08664 100644 --- a/source/module_ri/exx_abfs-construct_orbs.cpp +++ b/source/module_ri/exx_abfs-construct_orbs.cpp @@ -77,6 +77,7 @@ std::vector>> Exx_Abfs::Construct_ // P = f * Y std::vector>> Exx_Abfs::Construct_Orbs::abfs_same_atom( + const LCAO_Orbitals& orb, const std::vector>> &orbs, const double kmesh_times_mot, const double times_threshold ) @@ -99,7 +100,7 @@ std::vector>> Exx_Abfs::Construct_ #endif const std::vector>>> - abfs_same_atom_pca_psi = pca( abfs_same_atom, orbs, kmesh_times_mot, times_threshold ); + abfs_same_atom_pca_psi = pca( orb, abfs_same_atom, orbs, kmesh_times_mot, times_threshold ); #if TEST_EXX_LCAO==1 print_orbs(abfs_same_atom_pca_psi,"abfs_same_atom_pca_psi.dat"); @@ -255,6 +256,7 @@ std::vector>>> Exx_Abfs::Construct_O } std::vector>>> Exx_Abfs::Construct_Orbs::pca( + const LCAO_Orbitals& orb, const std::vector>> &abfs, const std::vector>> &orbs, const double kmesh_times_mot, @@ -264,7 +266,7 @@ std::vector>>> Exx_Abfs::Construct_O return std::vector>>>(abfs.size()); const std::vector,RI::Tensor>>> - eig = ABFs_Construct::PCA::cal_PCA( orbs, abfs, kmesh_times_mot ); + eig = ABFs_Construct::PCA::cal_PCA( orb, orbs, abfs, kmesh_times_mot ); const std::vector>>> psis = get_psi( abfs ); std::vector>>> psis_new( psis.size() ); @@ -484,4 +486,4 @@ void Exx_Abfs::Construct_Orbs::print_orbs_size( } os<>> abfs_same_atom( + const LCAO_Orbitals& orb, const std::vector>> &lcaos, const double kmesh_times_mot, const double times_threshold=0); @@ -40,6 +41,7 @@ class Exx_Abfs::Construct_Orbs const double norm_threshold = std::numeric_limits::min() ); static std::vector>>> pca( + const LCAO_Orbitals& orb, const std::vector>> &abfs, const std::vector>> &orbs, const double kmesh_times_mot, diff --git a/source/module_ri/exx_abfs-jle.cpp b/source/module_ri/exx_abfs-jle.cpp index 4e693a1eee..9a43c6632d 100644 --- a/source/module_ri/exx_abfs-jle.cpp +++ b/source/module_ri/exx_abfs-jle.cpp @@ -12,7 +12,7 @@ int Exx_Abfs::Jle::Lmax = 2; double Exx_Abfs::Jle::Ecut_exx = 60; double Exx_Abfs::Jle::tolerence = 1.0e-12; -void Exx_Abfs::Jle::init_jle( const double kmesh_times ) +void Exx_Abfs::Jle::init_jle( const double kmesh_times, const LCAO_Orbitals& orb ) { jle.resize( GlobalC::ucell.ntype ); @@ -22,35 +22,35 @@ void Exx_Abfs::Jle::init_jle( const double kmesh_times ) for (int L=0; L <= Lmax ; ++L) { const size_t ecut_number - = static_cast( sqrt( Ecut_exx ) * GlobalC::ORB.Phi[T].getRcut() / ModuleBase::PI ); // Rydberg Unit. + = static_cast( sqrt( Ecut_exx ) * orb.Phi[T].getRcut() / ModuleBase::PI ); // Rydberg Unit. jle[T][L].resize( ecut_number ); std::vector en(ecut_number, 0.0); - ModuleBase::Sphbes::Spherical_Bessel_Roots(ecut_number, L, tolerence, ModuleBase::GlobalFunc::VECTOR_TO_PTR(en), GlobalC::ORB.Phi[T].getRcut()); + ModuleBase::Sphbes::Spherical_Bessel_Roots(ecut_number, L, tolerence, ModuleBase::GlobalFunc::VECTOR_TO_PTR(en), orb.Phi[T].getRcut()); for(size_t E=0; E!=ecut_number; ++E) { - std::vector jle_r( GlobalC::ORB.Phi[T].PhiLN(0,0).getNr() ); + std::vector jle_r( orb.Phi[T].PhiLN(0,0).getNr() ); ModuleBase::Sphbes::Spherical_Bessel( - GlobalC::ORB.Phi[T].PhiLN(0,0).getNr(), - GlobalC::ORB.Phi[T].PhiLN(0,0).getRadial(), + orb.Phi[T].PhiLN(0,0).getNr(), + orb.Phi[T].PhiLN(0,0).getRadial(), en[E], L, ModuleBase::GlobalFunc::VECTOR_TO_PTR(jle_r)); jle[T][L][E].set_orbital_info( - GlobalC::ORB.Phi[T].PhiLN(0,0).getLabel(), - GlobalC::ORB.Phi[T].PhiLN(0,0).getType(), + orb.Phi[T].PhiLN(0,0).getLabel(), + orb.Phi[T].PhiLN(0,0).getType(), L, E, // N? - GlobalC::ORB.Phi[T].PhiLN(0,0).getNr(), - GlobalC::ORB.Phi[T].PhiLN(0,0).getRab(), - GlobalC::ORB.Phi[T].PhiLN(0,0).getRadial(), + orb.Phi[T].PhiLN(0,0).getNr(), + orb.Phi[T].PhiLN(0,0).getRab(), + orb.Phi[T].PhiLN(0,0).getRadial(), Numerical_Orbital_Lm::Psi_Type::Psi, ModuleBase::GlobalFunc::VECTOR_TO_PTR(jle_r), - static_cast(GlobalC::ORB.Phi[T].PhiLN(0,0).getNk() * kmesh_times) | 1, - GlobalC::ORB.Phi[T].PhiLN(0,0).getDk(), - GlobalC::ORB.Phi[T].PhiLN(0,0).getDruniform(), + static_cast(orb.Phi[T].PhiLN(0,0).getNk() * kmesh_times) | 1, + orb.Phi[T].PhiLN(0,0).getDk(), + orb.Phi[T].PhiLN(0,0).getDruniform(), false, true, PARAM.inp.cal_force); } diff --git a/source/module_ri/exx_abfs-jle.h b/source/module_ri/exx_abfs-jle.h index 43d3a9d80f..967df31de6 100644 --- a/source/module_ri/exx_abfs-jle.h +++ b/source/module_ri/exx_abfs-jle.h @@ -4,7 +4,7 @@ #include "exx_abfs.h" #include -#include "../module_basis/module_ao/ORB_atomic_lm.h" +#include "../module_basis/module_ao/ORB_read.h" class Exx_Abfs::Jle { @@ -15,7 +15,7 @@ class Exx_Abfs::Jle std::vector< Numerical_Orbital_Lm>>> jle; - void init_jle( const double kmesh_times ); + void init_jle( const double kmesh_times, const LCAO_Orbitals& orb ); static bool generate_matrix; static int Lmax; diff --git a/source/module_ri/exx_opt_orb-print.cpp b/source/module_ri/exx_opt_orb-print.cpp index 0a39a03f13..52505505d9 100644 --- a/source/module_ri/exx_opt_orb-print.cpp +++ b/source/module_ri/exx_opt_orb-print.cpp @@ -9,6 +9,7 @@ void Exx_Opt_Orb::print_matrix( const std::vector>> &matrix_S, const RI::Tensor &matrix_V, const size_t TA, const size_t IA, const size_t TB, const size_t IB, + const std::vector& orb_cutoff, const ModuleBase::Element_Basis_Index::Range &range_jles, const ModuleBase::Element_Basis_Index::IndexLNM &index_jles) const { @@ -64,9 +65,9 @@ void Exx_Opt_Orb::print_matrix( ofs << Exx_Abfs::Jle::Ecut_exx << " ecutwfc_jlq" << std::endl;//mohan modify 2009-09-08 if(TA==TB) - ofs << GlobalC::ORB.Phi[TA].getRcut() << " rcut_Jlq" << std::endl; + ofs << orb_cutoff[TA] << " rcut_Jlq" << std::endl; else - ofs << GlobalC::ORB.Phi[TA].getRcut() << " " << GlobalC::ORB.Phi[TB].getRcut() << " rcut_Jlq" << std::endl; + ofs << orb_cutoff[TA] << " " << orb_cutoff[TB] << " rcut_Jlq" << std::endl; // mohan add 'smooth' and 'smearing_sigma' 2009-08-28 ofs << 0 << " smooth" << std::endl; @@ -90,8 +91,8 @@ void Exx_Opt_Orb::print_matrix( const size_t nwfc = (TA==TB && IA==IB) ? cal_sum_M(TA) : cal_sum_M(TA)+cal_sum_M(TB); ofs << nwfc << " nwfc" << std::endl; - const size_t ecut_numberA = static_cast( sqrt( Exx_Abfs::Jle::Ecut_exx ) * GlobalC::ORB.Phi[TA].getRcut() / ModuleBase::PI ); // Rydberg Unit - const size_t ecut_numberB = static_cast( sqrt( Exx_Abfs::Jle::Ecut_exx ) * GlobalC::ORB.Phi[TB].getRcut() / ModuleBase::PI ); // Rydberg Unit + const size_t ecut_numberA = static_cast( sqrt( Exx_Abfs::Jle::Ecut_exx ) * orb_cutoff[TA] / ModuleBase::PI ); // Rydberg Unit + const size_t ecut_numberB = static_cast( sqrt( Exx_Abfs::Jle::Ecut_exx ) * orb_cutoff[TB] / ModuleBase::PI ); // Rydberg Unit if(TA==TB) ofs << ecut_numberA << " ne" << std::endl; else @@ -202,4 +203,4 @@ void Exx_Opt_Orb::print_matrix( print_S(ofs); print_V(ofs); ofs.close(); -} \ No newline at end of file +} diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index 3eeb975cbd..7f2d27f17d 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -13,8 +13,9 @@ #include "../module_ri/test_code/element_basis_index-test.h" #include "../module_ri/test_code/test_function.h" +#include -void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const +void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) const { // std::ofstream ofs_mpi(GlobalC::exx_lcao.test_dir.process+"time_"+ModuleBase::GlobalFunc::TO_STRING(GlobalV::MY_RANK),std::ofstream::app); @@ -22,15 +23,15 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const // ofs_mpi<<"memory:\t"<>> - lcaos = Exx_Abfs::Construct_Orbs::change_orbs( GlobalC::ORB, this->kmesh_times ); + lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->kmesh_times ); const std::vector>> - abfs = Exx_Abfs::Construct_Orbs::abfs_same_atom( lcaos, this->kmesh_times, GlobalC::exx_info.info_ri.pca_threshold ); + abfs = Exx_Abfs::Construct_Orbs::abfs_same_atom( orb, lcaos, this->kmesh_times, GlobalC::exx_info.info_ri.pca_threshold ); // ofs_mpi<<"memory:\t"<kmesh_times ); + jle.init_jle( this->kmesh_times, orb ); // ofs_mpi<<"memory:\t"< std::map>>>> { Matrix_Orbs22 m_lcaoslcaos_lcaoslcaos; - m_lcaoslcaos_lcaoslcaos.init( 1, this->kmesh_times, 1 ); + m_lcaoslcaos_lcaoslcaos.init( 1, orb, this->kmesh_times, 1 ); m_lcaoslcaos_lcaoslcaos.init_radial( lcaos, lcaos, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 m_lcaoslcaos_lcaoslcaos.init_radial_table(radial_R); @@ -87,7 +88,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const const auto ms_lcaoslcaos_jys = [&]() -> std::map>>>>> { Matrix_Orbs21 m_jyslcaos_lcaos; - m_jyslcaos_lcaos.init( 1, this->kmesh_times, 1 ); + m_jyslcaos_lcaos.init( 1, orb, this->kmesh_times, 1 ); m_jyslcaos_lcaos.init_radial( jle.jle, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 m_jyslcaos_lcaos.init_radial_table(radial_R); @@ -103,7 +104,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const const auto ms_jys_jys = [&]() -> std::map>>>> { Matrix_Orbs11 m_jys_jys; - m_jys_jys.init( 2, this->kmesh_times, 1 ); + m_jys_jys.init( 2, orb, this->kmesh_times, 1 ); m_jys_jys.init_radial( jle.jle, jle.jle ); #if TEST_EXX_RADIAL>=1 m_jys_jys.init_radial_table(radial_R); @@ -119,7 +120,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const const auto ms_abfs_abfs = [&]() -> std::map>>>> { Matrix_Orbs11 m_abfs_abfs; - m_abfs_abfs.init( 2, this->kmesh_times, 1 ); + m_abfs_abfs.init( 2, orb, this->kmesh_times, 1 ); m_abfs_abfs.init_radial( abfs, abfs ); #if TEST_EXX_RADIAL>=1 m_abfs_abfs.init_radial_table(radial_R); @@ -135,7 +136,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const const auto ms_lcaoslcaos_abfs = [&]() -> std::map>>>>> { Matrix_Orbs21 m_abfslcaos_lcaos; - m_abfslcaos_lcaos.init( 1, this->kmesh_times, 1 ); + m_abfslcaos_lcaos.init( 1, orb, this->kmesh_times, 1 ); m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 m_abfslcaos_lcaos.init_radial_table(radial_R); @@ -151,7 +152,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const const auto ms_jys_abfs = [&]() -> std::map>>>> { Matrix_Orbs11 m_jys_abfs; - m_jys_abfs.init( 2, this->kmesh_times, 1 ); + m_jys_abfs.init( 2, orb, this->kmesh_times, 1 ); m_jys_abfs.init_radial( jle.jle, abfs ); #if TEST_EXX_RADIAL>=1 m_jys_abfs.init_radial_table(radial_R); @@ -212,6 +213,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const m_jys_jys_proj, m_lcaoslcaos_lcaoslcaos_proj, TA, IA, TB, IB, + orb.cutoffs(), range_jys, index_jys ); } else @@ -222,6 +224,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const {{ms_jys_jys.at(T).at(I).at(T).at(I)}}, ms_lcaoslcaos_lcaoslcaos.at(T).at(I).at(T).at(I), TA, IA, TB, IB, + orb.cutoffs(), range_jys, index_jys ); } } @@ -278,6 +281,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const m_jys_jys_proj, m_lcaoslcaos_lcaoslcaos_proj, TA, IA, TB, IB, + orb.cutoffs(), range_jys, index_jys ); } else @@ -289,6 +293,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv) const {ms_jys_jys.at(TB).at(IB).at(TA).at(IA), ms_jys_jys.at(TB).at(IB).at(TB).at(IB)}}, ms_lcaoslcaos_lcaoslcaos.at(TA).at(IA).at(TB).at(IB), TA, IA, TB, IB, + orb.cutoffs(), range_jys, index_jys ); } } @@ -364,4 +369,4 @@ std::map>> Exx_Opt_Orb::get_radial_R() c radial_R[TB][TA].insert( delta_R ); } return radial_R; -} \ No newline at end of file +} diff --git a/source/module_ri/exx_opt_orb.h b/source/module_ri/exx_opt_orb.h index efcbb6c248..ba6910f2ba 100644 --- a/source/module_ri/exx_opt_orb.h +++ b/source/module_ri/exx_opt_orb.h @@ -4,6 +4,7 @@ #include "../module_base/matrix.h" #include "../module_base/element_basis_index.h" #include "module_cell/klist.h" +#include "module_basis/module_ao/ORB_read.h" #include #include #include @@ -12,7 +13,7 @@ class Exx_Opt_Orb { public: - void generate_matrix(const K_Vectors &kv) const; + void generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) const; private: std::vector>> cal_I( const std::map>>>> &ms, @@ -29,10 +30,11 @@ class Exx_Opt_Orb const std::vector>> &matrix_S, const RI::Tensor &matrix_V, const size_t TA, const size_t IA, const size_t TB, const size_t IB, + const std::vector& orb_cutoff, const ModuleBase::Element_Basis_Index::Range &range_jles, const ModuleBase::Element_Basis_Index::IndexLNM &index_jles) const; std::map>> get_radial_R() const; int kmesh_times = 4; }; -#endif \ No newline at end of file +#endif From 9c4a7b63169b775a1d7be48ddb7ae0150c8a9c20 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Wed, 11 Sep 2024 19:35:56 +0000 Subject: [PATCH 2/3] [pre-commit.ci lite] apply automatic fixes --- source/module_ri/ABFs_Construct-PCA.cpp | 26 +++++--- source/module_ri/Matrix_Orbs11.cpp | 84 ++++++++++++++++--------- source/module_ri/exx_opt_orb.cpp | 12 ++-- 3 files changed, 79 insertions(+), 43 deletions(-) diff --git a/source/module_ri/ABFs_Construct-PCA.cpp b/source/module_ri/ABFs_Construct-PCA.cpp index cf4c7b8fdf..31b1731bce 100644 --- a/source/module_ri/ABFs_Construct-PCA.cpp +++ b/source/module_ri/ABFs_Construct-PCA.cpp @@ -42,10 +42,13 @@ namespace PCA ModuleBase::TITLE("ABFs_Construct::PCA::get_sub_matrix"); assert(m.shape.size() == 3); RI::Tensor m_sub({ m.shape[0], m.shape[1], range[T][L].N }); - for (std::size_t ir=0; ir!=m.shape[0]; ++ir) - for (std::size_t jr=0; jr!=m.shape[1]; ++jr) - for (std::size_t N=0; N!=range[T][L].N; ++N) + for (std::size_t ir=0; ir!=m.shape[0]; ++ir) { + for (std::size_t jr=0; jr!=m.shape[1]; ++jr) { + for (std::size_t N=0; N!=range[T][L].N; ++N) { m_sub(ir, jr, N) = m(ir, jr, index[T][L][N][0]); +} +} +} m_sub = m_sub.reshape({ m.shape[0] * m.shape[1], range[T][L].N }); return m_sub; } @@ -57,11 +60,13 @@ namespace PCA for( std::size_t ic=0; ic!=m.shape[1]; ++ic ) { double sum=0; - for( std::size_t ir=0; ir!=m.shape[0]; ++ir ) + for( std::size_t ir=0; ir!=m.shape[0]; ++ir ) { sum += m(ir,ic); +} const double mean = sum/m.shape[0]; - for( std::size_t ir=0; ir!=m.shape[0]; ++ir ) + for( std::size_t ir=0; ir!=m.shape[0]; ++ir ) { m_new(ir,ic) = m(ir,ic) - mean; +} } return m_new; } @@ -86,16 +91,18 @@ namespace PCA const int Lmax_bak = GlobalC::exx_info.info_ri.abfs_Lmax; GlobalC::exx_info.info_ri.abfs_Lmax = std::numeric_limits::min(); - for( std::size_t T=0; T!=abfs.size(); ++T ) + for( std::size_t T=0; T!=abfs.size(); ++T ) { GlobalC::exx_info.info_ri.abfs_Lmax = std::max( GlobalC::exx_info.info_ri.abfs_Lmax, static_cast(abfs[T].size())-1 ); +} Matrix_Orbs21 m_abfslcaos_lcaos; m_abfslcaos_lcaos.init( 1, orb, kmesh_times, 1 ); m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos ); std::map>> delta_R; - for( std::size_t it=0; it!=abfs.size(); ++it ) + for( std::size_t it=0; it!=abfs.size(); ++it ) { delta_R[it][it] = {0.0}; +} m_abfslcaos_lcaos.init_radial_table(delta_R); GlobalC::exx_info.info_ri.abfs_Lmax = Lmax_bak; @@ -133,10 +140,11 @@ namespace PCA { for (int ic = 0; ic != m.shape[1]; ++ic) { - if (std::abs(m(ir, ic)) > threshold) + if (std::abs(m(ir, ic)) > threshold) { os << m(ir, ic) << "\t"; - else + } else { os << 0 << "\t"; +} } os << std::endl; } diff --git a/source/module_ri/Matrix_Orbs11.cpp b/source/module_ri/Matrix_Orbs11.cpp index 17fdb35238..76aab8a6e4 100644 --- a/source/module_ri/Matrix_Orbs11.cpp +++ b/source/module_ri/Matrix_Orbs11.cpp @@ -56,15 +56,21 @@ void Matrix_Orbs11::init_radial(const std::vectorMGT))); + Center2_Orb::Orb11(orb_A[TA][LA][NA], orb_B[TB][LB][NB], psb_, this->MGT))); +} +} +} +} +} +} ModuleBase::timer::tick("Matrix_Orbs11", "init_radial"); } @@ -72,18 +78,24 @@ void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& { ModuleBase::TITLE("Matrix_Orbs11", "init_radial"); ModuleBase::timer::tick("Matrix_Orbs11", "init_radial"); - for (size_t TA = 0; TA != orb_A.get_ntype(); ++TA) - for (size_t TB = 0; TB != orb_B.get_ntype(); ++TB) - for (int LA = 0; LA <= orb_A.Phi[TA].getLmax(); ++LA) - for (size_t NA = 0; NA != orb_A.Phi[TA].getNchi(LA); ++NA) - for (int LB = 0; LB <= orb_B.Phi[TB].getLmax(); ++LB) - for (size_t NB = 0; NB != orb_B.Phi[TB].getNchi(LB); ++NB) + for (size_t TA = 0; TA != orb_A.get_ntype(); ++TA) { + for (size_t TB = 0; TB != orb_B.get_ntype(); ++TB) { + for (int LA = 0; LA <= orb_A.Phi[TA].getLmax(); ++LA) { + for (size_t NA = 0; NA != orb_A.Phi[TA].getNchi(LA); ++NA) { + for (int LB = 0; LB <= orb_B.Phi[TB].getLmax(); ++LB) { + for (size_t NB = 0; NB != orb_B.Phi[TB].getNchi(LB); ++NB) { center2_orb11_s[TA][TB][LA][NA][LB].insert( std::make_pair(NB, Center2_Orb::Orb11(orb_A.Phi[TA].PhiLN(LA, NA), orb_B.Phi[TB].PhiLN(LB, NB), psb_, - this->MGT))); + this->MGT))); +} +} +} +} +} +} ModuleBase::timer::tick("Matrix_Orbs11", "init_radial"); } @@ -91,13 +103,19 @@ void Matrix_Orbs11::init_radial_table() { ModuleBase::TITLE("Matrix_Orbs11", "init_radial_table"); ModuleBase::timer::tick("Matrix_Orbs11", "init_radial_table"); - for (auto& coA: center2_orb11_s) - for (auto& coB: coA.second) - for (auto& coC: coB.second) - for (auto& coD: coC.second) - for (auto& coE: coD.second) - for (auto& coF: coE.second) - coF.second.init_radial_table(); + for (auto& coA: center2_orb11_s) { + for (auto& coB: coA.second) { + for (auto& coC: coB.second) { + for (auto& coD: coC.second) { + for (auto& coE: coD.second) { + for (auto& coF: coE.second) { + coF.second.init_radial_table(); +} +} +} +} +} +} ModuleBase::timer::tick("Matrix_Orbs11", "init_radial_table"); } @@ -105,7 +123,7 @@ void Matrix_Orbs11::init_radial_table(const std::map(position); - for (size_t i = 0; i != 4; ++i) - radials.insert(iq + i); + for (size_t i = 0; i != 4; ++i) { + radials.insert(iq + i); +} } - for (auto& coC: *center2_orb11_sAB) - for (auto& coD: coC.second) - for (auto& coE: coD.second) - for (auto& coF: coE.second) - coF.second.init_radial_table(radials); + for (auto& coC: *center2_orb11_sAB) { + for (auto& coD: coC.second) { + for (auto& coE: coD.second) { + for (auto& coF: coE.second) { + coF.second.init_radial_table(radials); +} +} +} +} } - } + } +} ModuleBase::timer::tick("Matrix_Orbs11", "init_radial_table"); } diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index 7f2d27f17d..d714de1a9d 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -36,8 +36,9 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) // ofs_mpi<<"memory:\t"<(abfs[T].size())-1 ); +} const ModuleBase::Element_Basis_Index::Range range_lcaos = Exx_Abfs::Abfs_Index::construct_range( lcaos ); const ModuleBase::Element_Basis_Index::IndexLNM index_lcaos = ModuleBase::Element_Basis_Index::construct_index( range_lcaos ); @@ -357,9 +358,9 @@ std::map>> Exx_Opt_Orb::get_radial_R() c { ModuleBase::TITLE("Exx_Opt_Orb::get_radial_R"); std::map>> radial_R; - for( size_t TA=0; TA!=GlobalC::ucell.ntype; ++TA ) - for( size_t IA=0; IA!=GlobalC::ucell.atoms[TA].na; ++IA ) - for( size_t TB=0; TB!=GlobalC::ucell.ntype; ++TB ) + for( size_t TA=0; TA!=GlobalC::ucell.ntype; ++TA ) { + for( size_t IA=0; IA!=GlobalC::ucell.atoms[TA].na; ++IA ) { + for( size_t TB=0; TB!=GlobalC::ucell.ntype; ++TB ) { for( size_t IB=0; IB!=GlobalC::ucell.atoms[TB].na; ++IB ) { const ModuleBase::Vector3 &tauA = GlobalC::ucell.atoms[TA].tau[IA]; @@ -368,5 +369,8 @@ std::map>> Exx_Opt_Orb::get_radial_R() c radial_R[TA][TB].insert( delta_R ); radial_R[TB][TA].insert( delta_R ); } +} +} +} return radial_R; } From a86763d9a85105222e9d229d22628edeb09d67fa Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Thu, 12 Sep 2024 17:56:25 +0800 Subject: [PATCH 3/3] remove GlobalC::ORB --- source/module_basis/module_ao/ORB_read.cpp | 9 -- source/module_basis/module_ao/ORB_read.h | 10 --- source/module_esolver/esolver_ks_lcao.cpp | 5 +- source/module_esolver/esolver_ks_lcao.h | 2 +- .../module_esolver/esolver_ks_lcao_tddft.cpp | 1 + source/module_esolver/lcao_gets.cpp | 7 +- .../hamilt_lcaodft/FORCE_STRESS.cpp | 1 + .../hamilt_lcaodft/hamilt_lcao.cpp | 16 +++- .../hamilt_lcaodft/hamilt_lcao.h | 3 +- .../operator_lcao/deepks_lcao.cpp | 11 +-- .../operator_lcao/dftu_force_stress.hpp | 1 - .../operator_lcao/dftu_lcao.cpp | 7 +- .../hamilt_lcaodft/operator_lcao/dftu_lcao.h | 3 + .../operator_lcao/ekinetic_new.cpp | 6 +- .../operator_lcao/ekinetic_new.h | 3 + .../operator_lcao/nonlocal_new.cpp | 6 +- .../operator_lcao/nonlocal_new.h | 4 + .../operator_lcao/overlap_new.cpp | 7 +- .../operator_lcao/overlap_new.h | 4 + .../operator_lcao/td_ekinetic_lcao.cpp | 7 +- .../operator_lcao/td_ekinetic_lcao.h | 3 + .../operator_lcao/td_nonlocal_lcao.cpp | 9 +- .../operator_lcao/td_nonlocal_lcao.h | 2 + .../operator_lcao/test/test_T_NL_cd.cpp | 2 + .../operator_lcao/test/test_dftu.cpp | 7 +- .../operator_lcao/test/test_ekineticnew.cpp | 4 +- .../operator_lcao/test/test_nonlocalnew.cpp | 5 +- .../operator_lcao/test/test_overlapnew.cpp | 4 +- .../operator_lcao/test/test_overlapnew_cd.cpp | 2 +- .../operator_lcao/test/tmp_mocks.cpp | 6 +- .../operator_lcao/veff_lcao.cpp | 3 +- .../hamilt_lcaodft/operator_lcao/veff_lcao.h | 11 ++- .../module_tddft/td_current.cpp | 13 ++- .../module_tddft/td_current.h | 3 + source/module_hsolver/diago_dav_subspace.cpp | 84 +++++++++++++++++-- source/module_io/td_current_io.cpp | 3 +- source/module_io/td_current_io.h | 1 + source/module_io/write_eband_terms.hpp | 13 +-- source/module_io/write_vxc.hpp | 4 +- source/module_lr/esolver_lrtd_lcao.cpp | 4 +- source/module_lr/esolver_lrtd_lcao.h | 1 + source/module_lr/hamilt_casida.h | 5 +- .../operator_casida/operator_lr_hxc.h | 9 +- 43 files changed, 205 insertions(+), 106 deletions(-) diff --git a/source/module_basis/module_ao/ORB_read.cpp b/source/module_basis/module_ao/ORB_read.cpp index 1d6f9195a1..faac0fb42e 100644 --- a/source/module_basis/module_ao/ORB_read.cpp +++ b/source/module_basis/module_ao/ORB_read.cpp @@ -15,10 +15,6 @@ //============================== /// PLEASE avoid using 'ORB' as global variable // mohan note 2021-03-23 -namespace GlobalC -{ -LCAO_Orbitals ORB; -} LCAO_Orbitals::LCAO_Orbitals() { @@ -41,11 +37,6 @@ LCAO_Orbitals::~LCAO_Orbitals() delete[] Alpha; } -const LCAO_Orbitals& LCAO_Orbitals::get_const_instance() -{ - return GlobalC::ORB; -} - std::vector LCAO_Orbitals::cutoffs() const { std::vector cutoffs(ntype); for (int it = 0; it < ntype; ++it) { diff --git a/source/module_basis/module_ao/ORB_read.h b/source/module_basis/module_ao/ORB_read.h index e714086fce..576759ac7b 100644 --- a/source/module_basis/module_ao/ORB_read.h +++ b/source/module_basis/module_ao/ORB_read.h @@ -22,9 +22,6 @@ class LCAO_Orbitals LCAO_Orbitals(); ~LCAO_Orbitals(); - // static function to get global instance - static const LCAO_Orbitals& get_const_instance(); - void init( std::ofstream& ofs_in, const int& ntype, @@ -134,11 +131,4 @@ class LCAO_Orbitals friend class TwoCenterBundle; // for the sake of TwoCenterBundle::to_LCAO_Orbitals }; -/// PLEASE avoid using 'ORB' as global variable -/// -///mohan note 2021 - 03 - 23 -namespace GlobalC -{ -extern LCAO_Orbitals ORB; -} #endif diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index d12bf618be..b343b63793 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -63,7 +63,7 @@ namespace ModuleESolver //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -ESolver_KS_LCAO::ESolver_KS_LCAO(): orb_(GlobalC::ORB) +ESolver_KS_LCAO::ESolver_KS_LCAO() { this->classname = "ESolver_KS_LCAO"; this->basisname = "LCAO"; @@ -463,6 +463,7 @@ void ESolver_KS_LCAO::after_all_runners() this->GG, this->GK, this->kv, + orb_.cutoffs(), this->pelec->wg, GlobalC::GridD #ifdef __EXX @@ -490,6 +491,7 @@ void ESolver_KS_LCAO::after_all_runners() this->kv, this->pelec->wg, GlobalC::GridD, + orb_.cutoffs(), this->two_center_bundle_ #ifdef __EXX , this->exx_lri_double ? &this->exx_lri_double->Hexxs : nullptr @@ -1266,6 +1268,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) this->kv.kvec_d, &hR, &GlobalC::ucell, + orb_.cutoffs(), &GlobalC::GridD, two_center_bundle_.kinetic_orb.get()); diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 735504dc4c..4d0cc10514 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -79,7 +79,7 @@ class ESolver_KS_LCAO : public ESolver_KS { TwoCenterBundle two_center_bundle_; // temporary introduced during removing GlobalC::ORB - LCAO_Orbitals& orb_; + LCAO_Orbitals orb_; // Temporarily store the stress to unify the interface with PW, // because it's hard to seperate force and stress calculation in LCAO. diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 39275e4659..21215b81e5 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -425,6 +425,7 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) kv, two_center_bundle_.overlap_orb.get(), tmp_DM->get_paraV_pointer(), + orb_, this->RA); } ESolver_KS_LCAO, double>::after_scf(istep); diff --git a/source/module_esolver/lcao_gets.cpp b/source/module_esolver/lcao_gets.cpp index b1ddf1ff3a..56b070f0b0 100644 --- a/source/module_esolver/lcao_gets.cpp +++ b/source/module_esolver/lcao_gets.cpp @@ -63,7 +63,8 @@ void ESolver_KS_LCAO, double>::get_S(void) this->p_hamilt = new hamilt::HamiltLCAO, double>( &this->pv, this->kv, - *(two_center_bundle_.overlap_orb)); + *(two_center_bundle_.overlap_orb), + orb_.cutoffs()); dynamic_cast, double>*>( this->p_hamilt->ops) ->contributeHR(); @@ -103,7 +104,9 @@ void ESolver_KS_LCAO, std::complex>::get_S(void) std::complex>( &this->pv, this->kv, - *(two_center_bundle_.overlap_orb)); + *(two_center_bundle_.overlap_orb), + orb_.cutoffs() + ); dynamic_cast< hamilt::OperatorLCAO, std::complex>*>( this->p_hamilt->ops) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 623a8af8e3..fd63fbfd42 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -257,6 +257,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, GlobalC::ucell, &GlobalC::GridD, two_center_bundle.overlap_orb_onsite.get(), + orb.cutoffs(), &GlobalC::dftu); tmp_dftu.cal_force_stress(isforce, isstress, force_dftu, stress_dftu); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index 8f94ba6135..692197a3ca 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -6,6 +6,7 @@ #include "module_base/timer.h" #include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#include #ifdef __DEEPKS #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" @@ -41,7 +42,7 @@ namespace hamilt { template -HamiltLCAO::HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& kv_in, const TwoCenterIntegrator& intor_overlap_orb) +HamiltLCAO::HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& kv_in, const TwoCenterIntegrator& intor_overlap_orb, const std::vector& orb_cutoff) { this->classname = "HamiltLCAO"; @@ -57,6 +58,7 @@ HamiltLCAO::HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& this->hR, this->sR, &GlobalC::ucell, + orb_cutoff, &GlobalC::GridD, &intor_overlap_orb); } @@ -130,6 +132,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, this->sR, &GlobalC::ucell, + orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb.get()); @@ -140,6 +143,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, &GlobalC::ucell, + orb.cutoffs(), &GlobalC::GridD, two_center_bundle.kinetic_orb.get()); this->getOperator()->add(ekinetic); @@ -153,6 +157,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, &GlobalC::ucell, + orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb_beta.get()); this->getOperator()->add(nonlocal); @@ -174,6 +179,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, pot_in, this->hR, // no explicit call yet &GlobalC::ucell, + orb.cutoffs(), &GlobalC::GridD ); this->getOperator()->add(veff); @@ -215,6 +221,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, GlobalC::ucell, &GlobalC::GridD, two_center_bundle.overlap_orb_onsite.get(), + orb.cutoffs(), &GlobalC::dftu); } this->getOperator()->add(dftu); @@ -240,6 +247,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, pot_in, this->hR, &GlobalC::ucell, + orb.cutoffs(), &GlobalC::GridD); // reset spin index and real space Hamiltonian matrix int start_spin = -1; @@ -256,6 +264,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, this->sR, &GlobalC::ucell, + orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb.get()); if (this->getOperator() == nullptr) @@ -275,6 +284,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, &GlobalC::ucell, + orb.cutoffs(), &GlobalC::GridD, two_center_bundle.kinetic_orb.get()); this->getOperator()->add(ekinetic); @@ -288,6 +298,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, &GlobalC::ucell, + orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb_beta.get()); // TDDFT velocity gague will calculate full non-local potential including the original one and the @@ -326,6 +337,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, kv, &GlobalC::ucell, + orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb.get()); this->getOperator()->add(td_ekinetic); @@ -334,6 +346,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, &GlobalC::ucell, + orb, &GlobalC::GridD); this->getOperator()->add(td_nonlocal); } @@ -355,6 +368,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, GlobalC::ucell, &GlobalC::GridD, two_center_bundle.overlap_orb_onsite.get(), + orb.cutoffs(), &GlobalC::dftu); } this->getOperator()->add(dftu); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h index 63a3b16d24..740fff3008 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h @@ -9,6 +9,7 @@ #include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" #include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp" +#include #ifdef __EXX #include "module_ri/Exx_LRI.h" #endif @@ -43,7 +44,7 @@ class HamiltLCAO : public Hamilt /** * @brief Constructor of vacuum Operators, only HR and SR will be initialed as empty HContainer */ - HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& kv_in, const TwoCenterIntegrator& intor_overlap_orb); + HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& kv_in, const TwoCenterIntegrator& intor_overlap_orb, const std::vector& orb_cutoff); ~HamiltLCAO() { diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp index 559ec5ac33..9c0e803b23 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp @@ -89,12 +89,11 @@ void hamilt::DeePKS>::initialize_HR(Grid_Driver* Gr const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad1]; const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; // choose the real adjacent atoms - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // Note: the distance of atoms should less than the cutoff radius, // When equal, the theoretical value of matrix element is zero, // but the calculated value is not zero due to the numerical error, which would lead to result changes. if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0 - < orb.Phi[T1].getRcut() + orb.Alpha[0].getRcut()) + < ptr_orb_->Phi[T1].getRcut() + ptr_orb_->Alpha[0].getRcut()) { is_adj[ad1] = true; } @@ -310,8 +309,6 @@ void hamilt::DeePKS>::calculate_HR() const Parallel_Orbitals* paraV = this->H_V_delta->get_paraV(); const int npol = this->ucell->get_npol(); - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); - // 1. calculate for each pair of atoms for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) { @@ -327,9 +324,9 @@ void hamilt::DeePKS>::calculate_HR() if (!GlobalV::deepks_equiv) { int ib = 0; - for (int L0 = 0; L0 <= orb.Alpha[0].getLmax(); ++L0) + for (int L0 = 0; L0 <= ptr_orb_->Alpha[0].getLmax(); ++L0) { - for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) + for (int N0 = 0; N0 < ptr_orb_->Alpha[0].getNchi(L0); ++N0) { const int inl = GlobalC::ld.get_inl(T0, I0, L0, N0); const double* pgedm = GlobalC::ld.get_gedms(inl); @@ -354,7 +351,7 @@ void hamilt::DeePKS>::calculate_HR() int nproj = 0; for (int il = 0; il < GlobalC::ld.get_lmaxd() + 1; il++) { - nproj += (2 * il + 1) * orb.Alpha[0].getNchi(il); + nproj += (2 * il + 1) * ptr_orb_->Alpha[0].getNchi(il); } for (int iproj = 0; iproj < nproj; iproj++) { diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_force_stress.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_force_stress.hpp index 77a51ffd16..e83c811a03 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_force_stress.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_force_stress.hpp @@ -64,7 +64,6 @@ void DFTU>::cal_force_stress(const bool cal_force, const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad]; const Atom* atom1 = &ucell->atoms[T1]; - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); auto all_indexes = paraV->get_indexes_row(iat1); auto col_indexes = paraV->get_indexes_col(iat1); // insert col_indexes into all_indexes to get universal set with no repeat elements diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.cpp index 9dd08f8440..c1ec37e3b1 100755 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.cpp @@ -18,8 +18,9 @@ hamilt::DFTU>::DFTU(HS_Matrix_K* hsk_in, const UnitCell& ucell_in, Grid_Driver* GridD_in, const TwoCenterIntegrator* intor, + const std::vector& orb_cutoff, ModuleDFTU::DFTU* dftu_in) - : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), intor_(intor) + : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), intor_(intor), orb_cutoff_(orb_cutoff) { this->cal_type = calculation_type::lcao_dftu; this->ucell = &ucell_in; @@ -71,12 +72,11 @@ void hamilt::DFTU>::initialize_HR(Grid_Driver* Grid const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad1]; const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; // choose the real adjacent atoms - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // Note: the distance of atoms should less than the cutoff radius, // When equal, the theoretical value of matrix element is zero, // but the calculated value is not zero due to the numerical error, which would lead to result changes. if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0 - < orb.Phi[T1].getRcut() + PARAM.inp.onsite_radius) + < orb_cutoff_[T1] + PARAM.inp.onsite_radius) { is_adj[ad1] = true; } @@ -120,7 +120,6 @@ void hamilt::DFTU>::cal_nlm_all(const Parallel_Orbi const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad]; const Atom* atom1 = &ucell->atoms[T1]; - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); auto all_indexes = paraV->get_indexes_row(iat1); auto col_indexes = paraV->get_indexes_col(iat1); // insert col_indexes into all_indexes to get universal set with no repeat elements diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h index f17471ac91..f60a172d62 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h @@ -32,6 +32,7 @@ class DFTU> : public OperatorLCAO const UnitCell& ucell_in, Grid_Driver* gridD_in, const TwoCenterIntegrator* intor, + const std::vector& orb_cutoff, ModuleDFTU::DFTU* dftu_in); ~DFTU>(); @@ -56,6 +57,8 @@ class DFTU> : public OperatorLCAO const TwoCenterIntegrator* intor_ = nullptr; + std::vector orb_cutoff_; + /// @brief the number of spin components, 1 for no-spin, 2 for collinear spin case and 4 for non-collinear spin case int nspin = 0; /// @brief the current spin index for nspin==2 to calculate spin-up and spin-down separately diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp index 0702a14d77..a68311ec4c 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp @@ -13,9 +13,10 @@ hamilt::EkineticNew>::EkineticNew( const std::vector>& kvec_d_in, hamilt::HContainer* hR_in, const UnitCell* ucell_in, + const std::vector& orb_cutoff, Grid_Driver* GridD_in, const TwoCenterIntegrator* intor) - : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), intor_(intor) + : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), orb_cutoff_(orb_cutoff), intor_(intor) { this->cal_type = calculation_type::lcao_fixed; this->ucell = ucell_in; @@ -66,12 +67,11 @@ void hamilt::EkineticNew>::initialize_HR(Grid_Drive } const ModuleBase::Vector3& R_index2 = adjs.box[ad1]; // choose the real adjacent atoms - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // Note: the distance of atoms should less than the cutoff radius, // When equal, the theoretical value of matrix element is zero, // but the calculated value is not zero due to the numerical error, which would lead to result changes. if (this->ucell->cal_dtau(iat1, iat2, R_index2).norm() * this->ucell->lat0 - < orb.Phi[T1].getRcut() + orb.Phi[T2].getRcut()) + < orb_cutoff_[T1] + orb_cutoff_[T2]) { is_adj[ad1] = true; } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.h index 2cb05b53eb..960985559f 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.h @@ -6,6 +6,7 @@ #include "module_cell/unitcell.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include namespace hamilt { @@ -43,6 +44,7 @@ class EkineticNew> : public OperatorLCAO const std::vector>& kvec_d_in, HContainer* hR_in, const UnitCell* ucell_in, + const std::vector& orb_cutoff, Grid_Driver* GridD_in, const TwoCenterIntegrator* intor); @@ -61,6 +63,7 @@ class EkineticNew> : public OperatorLCAO private: const UnitCell* ucell = nullptr; + std::vector orb_cutoff_; hamilt::HContainer* HR_fixed = nullptr; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.cpp index 7714cea230..17aea64b3d 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.cpp @@ -15,9 +15,10 @@ hamilt::NonlocalNew>::NonlocalNew( const std::vector>& kvec_d_in, hamilt::HContainer* hR_in, const UnitCell* ucell_in, + const std::vector& orb_cutoff, Grid_Driver* GridD_in, const TwoCenterIntegrator* intor) - : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), intor_(intor) + : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), orb_cutoff_(orb_cutoff), intor_(intor) { this->cal_type = calculation_type::lcao_fixed; this->ucell = ucell_in; @@ -66,12 +67,11 @@ void hamilt::NonlocalNew>::initialize_HR(Grid_Drive const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad1]; const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; // choose the real adjacent atoms - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // Note: the distance of atoms should less than the cutoff radius, // When equal, the theoretical value of matrix element is zero, // but the calculated value is not zero due to the numerical error, which would lead to result changes. if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0 - < orb.Phi[T1].getRcut() + this->ucell->infoNL.Beta[T0].get_rcut_max()) + < orb_cutoff_[T1] + this->ucell->infoNL.Beta[T0].get_rcut_max()) { is_adj[ad1] = true; } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h index 92fdc19bde..8206a28895 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h @@ -8,6 +8,7 @@ #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" #include +#include namespace hamilt { @@ -42,6 +43,7 @@ class NonlocalNew> : public OperatorLCAO const std::vector>& kvec_d_in, hamilt::HContainer* hR_in, const UnitCell* ucell_in, + const std::vector& orb_cutoff, Grid_Driver* GridD_in, const TwoCenterIntegrator* intor); ~NonlocalNew>(); @@ -57,6 +59,8 @@ class NonlocalNew> : public OperatorLCAO private: const UnitCell* ucell = nullptr; + std::vector orb_cutoff_; + hamilt::HContainer* HR_fixed = nullptr; // the following variable is introduced temporarily during LCAO refactoring diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp index 02737a7c2e..c29edfc0e5 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp @@ -5,6 +5,7 @@ #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h" +#include template hamilt::OverlapNew>::OverlapNew(HS_Matrix_K* hsk_in, @@ -12,9 +13,10 @@ hamilt::OverlapNew>::OverlapNew(HS_Matrix_K* hs hamilt::HContainer* hR_in, hamilt::HContainer* SR_in, const UnitCell* ucell_in, + const std::vector& orb_cutoff, Grid_Driver* GridD_in, const TwoCenterIntegrator* intor) - : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), intor_(intor) + : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), orb_cutoff_(orb_cutoff), intor_(intor) { this->cal_type = calculation_type::lcao_overlap; this->ucell = ucell_in; @@ -53,12 +55,11 @@ void hamilt::OverlapNew>::initialize_SR(Grid_Driver } const ModuleBase::Vector3& R_index = adjs.box[ad]; // choose the real adjacent atoms - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // Note: the distance of atoms should less than the cutoff radius, // When equal, the theoretical value of matrix element is zero, // but the calculated value is not zero due to the numerical error, which would lead to result changes. if (this->ucell->cal_dtau(iat1, iat2, R_index).norm() * this->ucell->lat0 - >= orb.Phi[T1].getRcut() + orb.Phi[T2].getRcut()) + >= orb_cutoff_[T1] + orb_cutoff_[T2]) { continue; } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.h index 8f6016cd21..93699cecea 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.h @@ -6,6 +6,7 @@ #include "module_cell/unitcell.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include namespace hamilt { @@ -41,6 +42,7 @@ class OverlapNew> : public OperatorLCAO hamilt::HContainer* hR_in, hamilt::HContainer* SR_in, const UnitCell* ucell_in, + const std::vector& orb_cutoff, Grid_Driver* GridD_in, const TwoCenterIntegrator* intor); @@ -53,6 +55,8 @@ class OverlapNew> : public OperatorLCAO private: const UnitCell* ucell = nullptr; + std::vector orb_cutoff_; + hamilt::HContainer* SR = nullptr; const TwoCenterIntegrator* intor_ = nullptr; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp index 9cac7a186b..a3dc0728f8 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp @@ -18,9 +18,10 @@ TDEkinetic>::TDEkinetic(HS_Matrix_K* hsk_in, hamilt::HContainer* hR_in, const K_Vectors* kv_in, const UnitCell* ucell_in, + const std::vector& orb_cutoff, Grid_Driver* GridD_in, const TwoCenterIntegrator* intor) - : kv(kv_in), OperatorLCAO(hsk_in, kv_in->kvec_d, hR_in), intor_(intor) + : OperatorLCAO(hsk_in, kv_in->kvec_d, hR_in), orb_cutoff_(orb_cutoff), kv(kv_in), intor_(intor) { this->ucell = ucell_in; this->cal_type = calculation_type::lcao_tddft_velocity; @@ -131,7 +132,6 @@ void TDEkinetic>::cal_HR_IJR(const int& iat1, std::complex* hr_mat_p, std::complex** current_mat_p) { - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // --------------------------------------------- // get info of orbitals of atom1 and atom2 from ucell // --------------------------------------------- @@ -276,12 +276,11 @@ void TDEkinetic>::initialize_HR(Grid_Driver* GridD) } const ModuleBase::Vector3& R_index2 = adjs.box[ad1]; // choose the real adjacent atoms - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // Note: the distance of atoms should less than the cutoff radius, // When equal, the theoretical value of matrix element is zero, // but the calculated value is not zero due to the numerical error, which would lead to result changes. if (this->ucell->cal_dtau(iat1, iat2, R_index2).norm() * this->ucell->lat0 - < orb.Phi[T1].getRcut() + orb.Phi[T2].getRcut()) + < orb_cutoff_[T1] + orb_cutoff_[T2]) { is_adj[ad1] = true; } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h index c2c01c6ced..c58082d07b 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h @@ -7,6 +7,7 @@ #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" #include "module_hamilt_lcao/module_tddft/td_velocity.h" #include "operator_lcao.h" +#include namespace hamilt @@ -38,6 +39,7 @@ class TDEkinetic> : public OperatorLCAO hamilt::HContainer* hR_in, const K_Vectors* kv_in, const UnitCell* ucell_in, + const std::vector& orb_cutoff, Grid_Driver* GridD_in, const TwoCenterIntegrator* intor); ~TDEkinetic(); @@ -82,6 +84,7 @@ class TDEkinetic> : public OperatorLCAO private: const UnitCell* ucell = nullptr; + std::vector orb_cutoff_; HContainer* SR = nullptr; /// @brief Store real space hamiltonian. TD term should include imaginary part, thus it has to be complex type. Only shared between TD operators. diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp index 50e703fa27..545a8fd95d 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp @@ -17,8 +17,9 @@ hamilt::TDNonlocal>::TDNonlocal(HS_Matrix_K* hs const std::vector>& kvec_d_in, hamilt::HContainer* hR_in, const UnitCell* ucell_in, + const LCAO_Orbitals& orb, Grid_Driver* GridD_in) - : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in) + : hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), orb_(orb) { this->cal_type = calculation_type::lcao_tddft_velocity; this->ucell = ucell_in; @@ -75,12 +76,11 @@ void hamilt::TDNonlocal>::initialize_HR(Grid_Driver const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad1]; const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; // choose the real adjacent atoms - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // Note: the distance of atoms should less than the cutoff radius, // When equal, the theoretical value of matrix element is zero, // but the calculated value is not zero due to the numerical error, which would lead to result changes. if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0 - < orb.Phi[T1].getRcut() + this->ucell->infoNL.Beta[T0].get_rcut_max()) + < orb_.Phi[T1].getRcut() + this->ucell->infoNL.Beta[T0].get_rcut_max()) { is_adj[ad1] = true; } @@ -155,7 +155,6 @@ void hamilt::TDNonlocal>::calculate_HR() const int iat1 = ucell->itia2iat(T1, I1); const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad]; const Atom* atom1 = &ucell->atoms[T1]; - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); auto all_indexes = paraV->get_indexes_row(iat1); auto col_indexes = paraV->get_indexes_col(iat1); all_indexes.insert(all_indexes.end(), col_indexes.begin(), col_indexes.end()); @@ -171,7 +170,7 @@ void hamilt::TDNonlocal>::calculate_HR() // snap_psibeta_half_tddft() are used to calculate // and as well if current are needed - module_tddft::snap_psibeta_half_tddft(orb, + module_tddft::snap_psibeta_half_tddft(orb_, this->ucell->infoNL, nlm, tau1 * this->ucell->lat0, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.h index 1fbcd48029..bef490a73a 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.h @@ -38,6 +38,7 @@ class TDNonlocal> : public OperatorLCAO const std::vector>& kvec_d_in, hamilt::HContainer* hR_in, const UnitCell* ucell_in, + const LCAO_Orbitals& orb, Grid_Driver* GridD_in); ~TDNonlocal>(); @@ -52,6 +53,7 @@ class TDNonlocal> : public OperatorLCAO private: const UnitCell* ucell = nullptr; + const LCAO_Orbitals& orb_; HContainer* HR = nullptr; /// @brief Store real space hamiltonian. TD term should include imaginary part, thus it has to be complex type. Only diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_T_NL_cd.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_T_NL_cd.cpp index 99d9218180..2a6ae0275a 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_T_NL_cd.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_T_NL_cd.cpp @@ -141,6 +141,7 @@ TEST_F(TNLTest, testTVNLcd2cd) kvec_d_in, HR, &ucell, + {1.0}, &gd, &intor_); hamilt::Operator>* op1 @@ -148,6 +149,7 @@ TEST_F(TNLTest, testTVNLcd2cd) kvec_d_in, HR, &ucell, + {1.0}, &gd, &intor_); // merge two Operators to a chain diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp index c1f49215d8..7372654a29 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp @@ -153,8 +153,6 @@ TEST_F(DFTUTest, constructHRd2d) hamilt::HS_Matrix_K hsk(paraV, true); hsk.set_zero_hk(); Grid_Driver gd(0, 0); - // check some input values - EXPECT_EQ(LCAO_Orbitals::get_const_instance().Phi[0].getRcut(), 1.0); // reset HR and DMR const double factor = 1.0 / test_nw / test_nw / test_size / test_size; for (int i = 0; i < DMR->get_nnr(); i++) @@ -164,7 +162,7 @@ TEST_F(DFTUTest, constructHRd2d) } std::chrono::high_resolution_clock::time_point start_time = std::chrono::high_resolution_clock::now(); hamilt::DFTU> - op(&hsk, kvec_d_in, HR, ucell, &gd, &intor_, &GlobalC::dftu); + op(&hsk, kvec_d_in, HR, ucell, &gd, &intor_, {1.0}, &GlobalC::dftu); std::chrono::high_resolution_clock::time_point end_time = std::chrono::high_resolution_clock::now(); std::chrono::duration elapsed_time = std::chrono::duration_cast>(end_time - start_time); @@ -221,7 +219,6 @@ TEST_F(DFTUTest, constructHRd2cd) hamilt::HS_Matrix_K> hsk(paraV, true); hsk.set_zero_hk(); Grid_Driver gd(0, 0); - EXPECT_EQ(LCAO_Orbitals::get_const_instance().Phi[0].getRcut(), 1.0); // reset HR and DMR const double factor = 0.5 / test_nw / test_nw / test_size / test_size; for (int i = 0; i < DMR->get_nnr(); i++) @@ -230,7 +227,7 @@ TEST_F(DFTUTest, constructHRd2cd) HR->get_wrapper()[i] = 0.0; } hamilt::DFTU, double>> - op(&hsk, kvec_d_in, HR, ucell, &gd, &intor_, &GlobalC::dftu); + op(&hsk, kvec_d_in, HR, ucell, &gd, &intor_, {1.0}, &GlobalC::dftu); op.contributeHR(); // check the occupations of dftu for spin-up for (int iat = 0; iat < test_size; iat++) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_ekineticnew.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_ekineticnew.cpp index 04a9ae288d..3176f3c8d3 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_ekineticnew.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_ekineticnew.cpp @@ -108,7 +108,7 @@ TEST_F(EkineticNewTest, constructHRd2d) hsk.set_zero_hk(); Grid_Driver gd(0, 0); hamilt::EkineticNew> - op(&hsk, kvec_d_in, HR, &ucell, &gd, &intor_); + op(&hsk, kvec_d_in, HR, &ucell, {1.0}, &gd, &intor_); op.contributeHR(); // check the value of HR for (int iap = 0; iap < HR->size_atom_pairs(); ++iap) @@ -158,7 +158,7 @@ TEST_F(EkineticNewTest, constructHRd2cd) hsk.set_zero_hk(); Grid_Driver gd(0, 0); hamilt::EkineticNew, double>> - op(&hsk, kvec_d_in, HR, &ucell, &gd, &intor_); + op(&hsk, kvec_d_in, HR, &ucell, {1.0}, &gd, &intor_); op.contributeHR(); // check the value of HR for (int iap = 0; iap < HR->size_atom_pairs(); ++iap) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_nonlocalnew.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_nonlocalnew.cpp index 2a4c9b5ade..001bbc1b5f 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_nonlocalnew.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_nonlocalnew.cpp @@ -134,10 +134,9 @@ TEST_F(NonlocalNewTest, constructHRd2d) Grid_Driver gd(0, 0); // check some input values EXPECT_EQ(ucell.infoNL.Beta[0].get_rcut_max(), 1.0); - EXPECT_EQ(LCAO_Orbitals::get_const_instance().Phi[0].getRcut(), 1.0); std::chrono::high_resolution_clock::time_point start_time = std::chrono::high_resolution_clock::now(); hamilt::NonlocalNew> - op(&hsk, kvec_d_in, HR, &ucell, &gd, &intor_); + op(&hsk, kvec_d_in, HR, &ucell, {1.0}, &gd, &intor_); std::chrono::high_resolution_clock::time_point end_time = std::chrono::high_resolution_clock::now(); std::chrono::duration elapsed_time = std::chrono::duration_cast>(end_time - start_time); @@ -207,7 +206,7 @@ TEST_F(NonlocalNewTest, constructHRd2cd) hsk.set_zero_hk(); Grid_Driver gd(0, 0); hamilt::NonlocalNew, double>> - op(&hsk, kvec_d_in, HR, &ucell, &gd, &intor_); + op(&hsk, kvec_d_in, HR, &ucell, {1.0}, &gd, &intor_); op.contributeHR(); // check the value of HR for (int iap = 0; iap < HR->size_atom_pairs(); ++iap) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew.cpp index 7d1b256e69..89fdcd5134 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew.cpp @@ -108,7 +108,7 @@ TEST_F(OverlapNewTest, constructHRd2d) hsk.set_zero_sk(); Grid_Driver gd(0, 0); hamilt::OverlapNew> - op(&hsk, kvec_d_in, nullptr, SR, &ucell, &gd, &intor_); + op(&hsk, kvec_d_in, nullptr, SR, &ucell, {1.0}, &gd, &intor_); op.contributeHR(); // check the value of SR for (int iap = 0; iap < SR->size_atom_pairs(); ++iap) @@ -142,7 +142,7 @@ TEST_F(OverlapNewTest, constructHRd2cd) hsk.set_zero_sk(); Grid_Driver gd(0, 0); hamilt::OverlapNew, double>> - op(&hsk, kvec_d_in, nullptr, SR, &ucell, &gd, &intor_); + op(&hsk, kvec_d_in, nullptr, SR, &ucell, {1.0}, &gd, &intor_); op.contributeHR(); // check the value of SR for (int iap = 0; iap < SR->size_atom_pairs(); ++iap) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew_cd.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew_cd.cpp index ffbe3bbb8a..15ad0893ad 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew_cd.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew_cd.cpp @@ -108,7 +108,7 @@ TEST_F(OverlapNewTest, constructHRcd2cd) hsk.set_zero_sk(); Grid_Driver gd(0, 0); hamilt::OverlapNew, std::complex>> - op(&hsk, kvec_d_in, nullptr, SR, &ucell, &gd, &intor_); + op(&hsk, kvec_d_in, nullptr, SR, &ucell, {1.0}, &gd, &intor_); op.contributeHR(); // check the value of SR for (int iap = 0; iap < SR->size_atom_pairs(); ++iap) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/tmp_mocks.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/tmp_mocks.cpp index dfaced0221..38f09cf72c 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/tmp_mocks.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/tmp_mocks.cpp @@ -147,10 +147,6 @@ void TwoCenterIntegrator::snap( } #include "module_basis/module_ao/ORB_read.h" -const LCAO_Orbitals& LCAO_Orbitals::get_const_instance() { - static LCAO_Orbitals instance; - return instance; -} LCAO_Orbitals::LCAO_Orbitals() { this->Phi = new Numerical_Orbital[1]; } LCAO_Orbitals::~LCAO_Orbitals() { delete[] Phi; } @@ -205,4 +201,4 @@ void Numerical_Orbital::set_orbital_info(const int&, const std::string&, const int&, const int*, - const int&) {} \ No newline at end of file + const int&) {} diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp index 0a2600771c..23755e497b 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp @@ -38,12 +38,11 @@ void Veff>::initialize_HR(const UnitCell* ucell_in, } const ModuleBase::Vector3& R_index2 = adjs.box[ad1]; // choose the real adjacent atoms - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // Note: the distance of atoms should less than the cutoff radius, // When equal, the theoretical value of matrix element is zero, // but the calculated value is not zero due to the numerical error, which would lead to result changes. if (ucell_in->cal_dtau(iat1, iat2, R_index2).norm() * ucell_in->lat0 - < orb.Phi[T1].getRcut() + orb.Phi[T2].getRcut()) + < orb_cutoff_[T1] + orb_cutoff_[T2]) { hamilt::AtomPair tmp(iat1, iat2, R_index2, paraV); this->hR->insert_pair(tmp); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h index 9c97e1020f..ee92f87777 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h @@ -7,6 +7,7 @@ #include "operator_lcao.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_cell/unitcell.h" +#include namespace hamilt { @@ -40,8 +41,9 @@ class Veff> : public OperatorLCAO elecstate::Potential* pot_in, hamilt::HContainer* hR_in, const UnitCell* ucell_in, + const std::vector& orb_cutoff, Grid_Driver* GridD_in) - : GK(GK_in), pot(pot_in), ucell(ucell_in), + : GK(GK_in), orb_cutoff_(orb_cutoff), pot(pot_in), ucell(ucell_in), gd(GridD_in), OperatorLCAO(hsk_in, kvec_d_in, hR_in) { this->cal_type = calculation_type::lcao_gint; @@ -59,8 +61,9 @@ class Veff> : public OperatorLCAO elecstate::Potential* pot_in, hamilt::HContainer* hR_in, const UnitCell* ucell_in, + const std::vector& orb_cutoff, Grid_Driver* GridD_in) - : GG(GG_in), pot(pot_in), OperatorLCAO(hsk_in, kvec_d_in, hR_in) + : GG(GG_in), orb_cutoff_(orb_cutoff), pot(pot_in), OperatorLCAO(hsk_in, kvec_d_in, hR_in) { this->cal_type = calculation_type::lcao_gint; this->initialize_HR(ucell_in, GridD_in); @@ -87,6 +90,8 @@ class Veff> : public OperatorLCAO // used for gamma only algorithms. Gint_Gamma* GG = nullptr; + std::vector orb_cutoff_; + // Charge calculating method in LCAO base and contained grid base calculation: DM_R, DM, pvpR_reduced elecstate::Potential* pot = nullptr; @@ -104,4 +109,4 @@ class Veff> : public OperatorLCAO }; } // namespace hamilt -#endif \ No newline at end of file +#endif diff --git a/source/module_hamilt_lcao/module_tddft/td_current.cpp b/source/module_hamilt_lcao/module_tddft/td_current.cpp index 6b8483fe5c..08491769f1 100644 --- a/source/module_hamilt_lcao/module_tddft/td_current.cpp +++ b/source/module_hamilt_lcao/module_tddft/td_current.cpp @@ -10,8 +10,9 @@ TD_current::TD_current(const UnitCell* ucell_in, Grid_Driver* GridD_in, const Parallel_Orbitals* paraV, + const LCAO_Orbitals& orb, const TwoCenterIntegrator* intor) - : ucell(ucell_in), Grid(GridD_in), paraV(paraV) , intor_(intor) + : ucell(ucell_in), paraV(paraV) , orb_(orb), Grid(GridD_in), intor_(intor) { // for length gague, the A(t) = 0 for all the time. this->cart_At = ModuleBase::Vector3(0,0,0); @@ -54,12 +55,11 @@ void TD_current::initialize_vcomm_r(Grid_Driver* GridD, const Parallel_Orbitals* const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad1]; const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; // choose the real adjacent atoms - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // Note: the distance of atoms should less than the cutoff radius, // When equal, the theoretical value of matrix element is zero, // but the calculated value is not zero due to the numerical error, which would lead to result changes. if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0 - < orb.Phi[T1].getRcut() + this->ucell->infoNL.Beta[T0].get_rcut_max()) + < orb_.Phi[T1].getRcut() + this->ucell->infoNL.Beta[T0].get_rcut_max()) { is_adj[ad1] = true; } @@ -131,12 +131,11 @@ void TD_current::initialize_grad_term(Grid_Driver* GridD, const Parallel_Orbital } const ModuleBase::Vector3& R_index2 = adjs.box[ad1]; // choose the real adjacent atoms - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // Note: the distance of atoms should less than the cutoff radius, // When equal, the theoretical value of matrix element is zero, // but the calculated value is not zero due to the numerical error, which would lead to result changes. if (this->ucell->cal_dtau(iat1, iat2, R_index2).norm() * this->ucell->lat0 - < orb.Phi[T1].getRcut() + orb.Phi[T2].getRcut()) + < orb_.Phi[T1].getRcut() + orb_.Phi[T2].getRcut()) { is_adj[ad1] = true; } @@ -205,7 +204,6 @@ void TD_current::calculate_vcomm_r() const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad]; const Atom* atom1 = &ucell->atoms[T1]; - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); auto all_indexes = paraV->get_indexes_row(iat1); #ifdef _OPENMP if(atom_row_list.find(iat1) == atom_row_list.end()) @@ -229,7 +227,7 @@ void TD_current::calculate_vcomm_r() // snap_psibeta_half_tddft() are used to calculate // and as well if current are needed - module_tddft::snap_psibeta_half_tddft(orb, + module_tddft::snap_psibeta_half_tddft(orb_, this->ucell->infoNL, nlm, tau1 * this->ucell->lat0, @@ -427,7 +425,6 @@ void TD_current::cal_grad_IJR(const int& iat1, const ModuleBase::Vector3& dtau, std::complex** current_mat_p) { - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // --------------------------------------------- // get info of orbitals of atom1 and atom2 from ucell // --------------------------------------------- diff --git a/source/module_hamilt_lcao/module_tddft/td_current.h b/source/module_hamilt_lcao/module_tddft/td_current.h index f4a1fe2804..9c3b2d362e 100644 --- a/source/module_hamilt_lcao/module_tddft/td_current.h +++ b/source/module_hamilt_lcao/module_tddft/td_current.h @@ -17,6 +17,7 @@ class TD_current TD_current(const UnitCell* ucell_in, Grid_Driver* GridD_in, const Parallel_Orbitals* paraV, + const LCAO_Orbitals& orb, const TwoCenterIntegrator* intor); ~TD_current(); @@ -33,6 +34,8 @@ class TD_current const Parallel_Orbitals* paraV = nullptr; + const LCAO_Orbitals& orb_; + Grid_Driver* Grid = nullptr; /// @brief Store real space hamiltonian. TD term should include imaginary part, thus it has to be complex type. Only shared between TD operators. std::vector>*> current_term = {nullptr, nullptr, nullptr}; diff --git a/source/module_hsolver/diago_dav_subspace.cpp b/source/module_hsolver/diago_dav_subspace.cpp index 3c5311adb2..f6d99d8447 100644 --- a/source/module_hsolver/diago_dav_subspace.cpp +++ b/source/module_hsolver/diago_dav_subspace.cpp @@ -20,8 +20,8 @@ Diago_DavSubspace::Diago_DavSubspace(const std::vector& precond const int& diag_nmax_in, const bool& need_subspace_in, const diag_comm_info& diag_comm_in) - : precondition(precondition_in), n_band(nband_in), dim(nbasis_in), nbase_x(nband_in * david_ndim_in), - diag_thr(diag_thr_in), iter_nmax(diag_nmax_in), is_subspace(need_subspace_in), diag_comm(diag_comm_in) + : n_band(nband_in), dim(nbasis_in), nbase_x(nband_in * david_ndim_in), + diag_thr(diag_thr_in), iter_nmax(diag_nmax_in), is_subspace(need_subspace_in), diag_comm(diag_comm_in), precondition(precondition_in) { this->device = base_device::get_device_type(this->ctx); @@ -97,7 +97,7 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, // convflag[m] = true if the m th band is convergent std::vector convflag(this->n_band, false); - // unconv[m] store the number of the m th unconvergent band + // unconv[m] store the index of the m th unconvergent band std::vector unconv(this->n_band); // the dimension of the reduced psi set @@ -108,6 +108,9 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, ModuleBase::timer::tick("Diago_DavSubspace", "first"); + // copy psi_in to psi_in_iter for the first n_band bands + // NOTE: not the first nbase_x (nband * david_ndim) bands + // since psi_in_iter is zero-initialized, the rest bands are zero for (int m = 0; m < this->n_band; m++) { unconv[m] = m; @@ -119,10 +122,15 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, this->dim); } + // compute h*psi_in_iter + // NOTE: bands after the first n_band should yield zero hpsi_func(this->hphi, this->psi_in_iter, this->nbase_x, this->dim, 0, this->nbase_x - 1); + // at this stage, notconv = n_band and nbase = 0 + // note that nbase of cal_elem is an inout parameter: nbase := nbase + notconv this->cal_elem(this->dim, nbase, this->notconv, this->psi_in_iter, this->hphi, this->hcc, this->scc); + // generalized eigenvalue problem (hcc, scc) for the first n_band bands this->diag_zhegvx(nbase, this->n_band, this->hcc, @@ -146,6 +154,7 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, { dav_iter++; + // nbase is the effective subspace size this->cal_grad(hpsi_func, this->dim, nbase, @@ -221,6 +230,26 @@ int Diago_DavSubspace::diag_once(const HPsiFunc& hpsi_func, { // overall convergence or last iteration: exit the iteration + // update this->psi_in_iter according to psi_in + for (size_t i = 0; i < this->n_band; i++) + { + syncmem_complex_op()(this->ctx, + this->ctx, + this->psi_in_iter + i * this->dim, + psi_in + i * psi_in_dmax, + this->dim); + } + + this->refresh(this->dim, + this->n_band, + nbase, + eigenvalue_in_hsolver, + this->psi_in_iter, + this->hphi, + this->hcc, + this->scc, + this->vcc); + ModuleBase::timer::tick("Diago_DavSubspace", "last"); break; } @@ -283,6 +312,8 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, } } + // psi_iter[:, nbase:nbase+notconv] + // = psi_iter[:, :nbase] * vcc[:nbase, :notconv] gemm_op()(this->ctx, 'N', 'N', @@ -298,6 +329,9 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, psi_iter + nbase * this->dim, this->dim); + // for psi_iter[:, nbase:nbase+notconv], + // each column is multiplied by the corresponding (minus) eigenvalue + // NOTE: eigenvalue_iter[m] correspond to psi_iter[:, nbase+m] (to be verified) for (int m = 0; m < notconv; m++) { @@ -310,6 +344,12 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, e_temp_cpu.data()); } + // psi_iter[:, nbase:nbase+notconv] += hphi[:, :nbase] * vcc[:nbase, :notconv] + // + // in terms of input, psi_iter should be + // psi_iter[:, nbase:nbase+notconv] := + // hpsi[:, :nbase] * vcc[:nbase, :notconv] - e[]*psi_iter[:, nbase:nbase+notconv]*vcc[:nbase, :notconv] + // TODO two gemm operations can be combined into one gemm_op()(this->ctx, 'N', 'N', @@ -325,13 +365,31 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, psi_iter + (nbase) * this->dim, this->dim); + // FIXME + // (QE dav solver -> isolver = 0) + // In QE, for davidson, h_diag is computed by g2kin + v_of_0 + // this is passed to g_psi, where the residual vector is element-wise + // multiplied by (h_diag - e * s_siag) + // + // There is a macro "TEST_NEW_PRECONDITIONING" in QE, which + // takes x = h_diag - e * s_siag and define + // denm = 0.5 * ( 1 + x + sqrt(1 + (x-1)^2) ) + // which then precondition the residual vector by psi := psi ./ denm + // note that this preconditioning recovers x when x is large, and + // approaches 0.5*(1+sqrt(2)) when x is close to 0. + // + // + // In the current version of ABACUS, "precondition" is simply g2kin, + // and the tested preconditioning of QE is always active. + // // "precondition!!!" std::vector pre(this->dim, 0.0); for (int m = 0; m < notconv; m++) { for (size_t i = 0; i < this->dim; i++) { - double x = this->precondition[i] - (*eigenvalue_iter)[m]; + //double x = this->precondition[i] - (*eigenvalue_iter)[m]; + double x = std::abs(this->precondition[i] - (*eigenvalue_iter)[m]); pre[i] = 0.5 * (1.0 + x + sqrt(1 + (x - 1.0) * (x - 1.0))); } vector_div_vector_op()(this->ctx, @@ -342,6 +400,8 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, } // "normalize!!!" in order to improve numerical stability of subspace diagonalization + + // column-wise normalize psi_iter[:, nbase:nbase+notconv] std::vector psi_norm(notconv, 0.0); for (size_t i = 0; i < notconv; i++) { @@ -349,7 +409,7 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, this->dim, psi_iter + (nbase + i) * this->dim, psi_iter + (nbase + i) * this->dim, - false); + false); // FIXME why reduce is false? assert(psi_norm[i] > 0.0); psi_norm[i] = sqrt(psi_norm[i]); @@ -360,6 +420,7 @@ void Diago_DavSubspace::cal_grad(const HPsiFunc& hpsi_func, psi_norm[i]); } + // update hpsi[:, nbase:nbase+notconv] hpsi_func(&hphi[nbase * this->dim], psi_iter, this->nbase_x, this->dim, nbase, nbase + notconv - 1); ModuleBase::timer::tick("Diago_DavSubspace", "cal_grad"); @@ -377,6 +438,9 @@ void Diago_DavSubspace::cal_elem(const int& dim, { ModuleBase::timer::tick("Diago_DavSubspace", "cal_elem"); + // hcc[:, nbase:] = psi_iter.H * hpsi_iter[:, nbase:] + // while psi_iter has nbase_x columns, only the first nbase + notconv columns are + // involved in the matrix multiplication gemm_op()(this->ctx, 'C', 'N', @@ -392,6 +456,7 @@ void Diago_DavSubspace::cal_elem(const int& dim, &hcc[nbase * this->nbase_x], this->nbase_x); + // scc[:, :nbase] := psi_iter.H * psi_iter[:, nbase:] gemm_op()(this->ctx, 'C', 'N', @@ -471,20 +536,28 @@ void Diago_DavSubspace::cal_elem(const int& dim, const size_t last_nbase = nbase; // init: last_nbase = 0 nbase = nbase + notconv; + + // NOTE: nbase is an in & out parameter! nbase is updated to nbase + notconv + for (size_t i = 0; i < nbase; i++) { if (i >= last_nbase) { + // ensure the diagonal elements of hcc and scc are real + // NOTE: set_real_tocomplex convert a real to itself, and a complex to a new complex + // whose real part is the input's real part and imaginary part is zero hcc[i * this->nbase_x + i] = set_real_tocomplex(hcc[i * this->nbase_x + i]); scc[i * this->nbase_x + i] = set_real_tocomplex(scc[i * this->nbase_x + i]); } for (size_t j = std::max(i + 1, last_nbase); j < nbase; j++) { + // ensure the hermicity of hcc and scc hcc[i * this->nbase_x + j] = get_conj(hcc[j * this->nbase_x + i]); scc[i * this->nbase_x + j] = get_conj(scc[j * this->nbase_x + i]); } } + // make hcc[nbase:, nbase:] and scc[nbase:, nbase:] to be zero for (size_t i = nbase; i < this->nbase_x; i++) { for (size_t j = nbase; j < this->nbase_x; j++) @@ -743,6 +816,7 @@ int Diago_DavSubspace::diag(const HPsiFunc& hpsi_func, do { + printf("enter diag... is_subspace = %d, ntry = %d\n", this->is_subspace, ntry); if (this->is_subspace || ntry > 0) { this->diagH_subspace(psi_in, eigenvalue_in_hsolver, hpsi_func, this->n_band, this->dim, psi_in_dmax); diff --git a/source/module_io/td_current_io.cpp b/source/module_io/td_current_io.cpp index 9bd08297e8..ddd93216c7 100644 --- a/source/module_io/td_current_io.cpp +++ b/source/module_io/td_current_io.cpp @@ -122,6 +122,7 @@ void ModuleIO::write_current(const int istep, const K_Vectors& kv, const TwoCenterIntegrator* intor, const Parallel_Orbitals* pv, + const LCAO_Orbitals& orb, Record_adj& ra) { @@ -131,7 +132,7 @@ void ModuleIO::write_current(const int istep, std::vector>*> current_term = {nullptr, nullptr, nullptr}; if (!TD_Velocity::tddft_velocity) { - cal_current = new TD_current(&GlobalC::ucell, &GlobalC::GridD, pv, intor); + cal_current = new TD_current(&GlobalC::ucell, &GlobalC::GridD, pv, orb, intor); cal_current->calculate_vcomm_r(); cal_current->calculate_grad_term(); for (int dir = 0; dir < 3; dir++) diff --git a/source/module_io/td_current_io.h b/source/module_io/td_current_io.h index 17f34928d6..74acc03134 100644 --- a/source/module_io/td_current_io.h +++ b/source/module_io/td_current_io.h @@ -16,6 +16,7 @@ void write_current(const int istep, const K_Vectors& kv, const TwoCenterIntegrator* intor, const Parallel_Orbitals* pv, + const LCAO_Orbitals& orb, Record_adj& ra); /// @brief calculate sum_n[𝜌_(𝑛𝑘,𝜇𝜈)] for current calculation diff --git a/source/module_io/write_eband_terms.hpp b/source/module_io/write_eband_terms.hpp index 9cde57469b..97ca13f2f9 100644 --- a/source/module_io/write_eband_terms.hpp +++ b/source/module_io/write_eband_terms.hpp @@ -23,6 +23,7 @@ namespace ModuleIO const K_Vectors& kv, const ModuleBase::matrix& wg, Grid_Driver& gd, + const std::vector& orb_cutoff, const TwoCenterBundle& two_center_bundle #ifdef __EXX , std::vector>>>* Hexxd = nullptr @@ -60,7 +61,7 @@ namespace ModuleIO hamilt::HContainer kinetic_R_ao(pv); if_gamma_fix(kinetic_R_ao); hamilt::EkineticNew> kinetic_op(&kinetic_k_ao, kv.kvec_d, - &kinetic_R_ao, &ucell, &gd, two_center_bundle.kinetic_orb.get()); + &kinetic_R_ao, &ucell, orb_cutoff, &gd, two_center_bundle.kinetic_orb.get()); kinetic_op.contributeHR(); std::vector> e_orb_kinetic; for (int ik = 0;ik < kv.get_nks();++ik) @@ -86,7 +87,7 @@ namespace ModuleIO hamilt::HContainer v_pp_local_R_ao(pv); if_gamma_fix(v_pp_local_R_ao); std::vector> e_orb_pp_local; - hamilt::Veff> v_pp_local_op(gint, &v_pp_local_k_ao, kv.kvec_d, &pot_local, &v_pp_local_R_ao, &ucell, &gd); + hamilt::Veff> v_pp_local_op(gint, &v_pp_local_k_ao, kv.kvec_d, &pot_local, &v_pp_local_R_ao, &ucell, orb_cutoff, &gd); v_pp_local_op.contributeHR(); for (int ik = 0;ik < kv.get_nks();++ik) { @@ -109,7 +110,7 @@ namespace ModuleIO if_gamma_fix(v_pp_nonlocal_R_ao); std::vector> e_orb_pp_nonlocal; hamilt::NonlocalNew> v_pp_nonlocal_op(&v_pp_nonlocal_k_ao, kv.kvec_d, - &v_pp_nonlocal_R_ao, &ucell, &gd, two_center_bundle.overlap_orb_beta.get()); + &v_pp_nonlocal_R_ao, &ucell, orb_cutoff, &gd, two_center_bundle.overlap_orb_beta.get()); v_pp_nonlocal_op.contributeHR(); for (int ik = 0;ik < kv.get_nks();++ik) { @@ -141,7 +142,7 @@ namespace ModuleIO for (int is = 0; is < nspin0; ++is) { v_hartree_op[is] = new hamilt::Veff>(gint, - &v_hartree_k_ao, kv.kvec_d, &pot_hartree, &v_hartree_R_ao[is], &ucell, &gd); + &v_hartree_k_ao, kv.kvec_d, &pot_hartree, &v_hartree_R_ao[is], &ucell, orb_cutoff, &gd); v_hartree_op[is]->contributeHR(); } std::vector> e_orb_hartree; @@ -163,7 +164,7 @@ namespace ModuleIO // 5. xc (including exx) if (!PARAM.inp.out_mat_xc) // avoid duplicate output { - write_Vxc(nspin, nbasis, drank, pv, psi, ucell, sf, rho_basis, rhod_basis, vloc, chg, gint_gamma, gint_k, kv, wg, gd + write_Vxc(nspin, nbasis, drank, pv, psi, ucell, sf, rho_basis, rhod_basis, vloc, chg, gint_gamma, gint_k, kv, orb_cutoff, wg, gd #ifdef __EXX , Hexxd, Hexxc #endif @@ -171,4 +172,4 @@ namespace ModuleIO } } } -#endif \ No newline at end of file +#endif diff --git a/source/module_io/write_vxc.hpp b/source/module_io/write_vxc.hpp index 5314946941..61527de4b3 100644 --- a/source/module_io/write_vxc.hpp +++ b/source/module_io/write_vxc.hpp @@ -195,6 +195,7 @@ void write_Vxc(const int nspin, Gint_Gamma& gint_gamma, // mohan add 2024-04-01 Gint_k& gint_k, // mohan add 2024-04-01 const K_Vectors& kv, + const std::vector& orb_cutoff, const ModuleBase::matrix& wg, Grid_Driver& gd #ifdef __EXX @@ -242,6 +243,7 @@ void write_Vxc(const int nspin, potxc, &vxcs_R_ao[is], &ucell, + orb_cutoff, &gd); vxcs_op_ao[is]->contributeHR(); @@ -336,4 +338,4 @@ void write_Vxc(const int nspin, } } } // namespace ModuleIO -#endif \ No newline at end of file +#endif diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 8df436d3bf..e9f161ea53 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -206,6 +206,7 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol } #endif this->pelec = new elecstate::ElecStateLCAO(); + orb_cutoff_ = ks_sol.orb_.cutoffs(); } template @@ -242,6 +243,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu LCAO_Orbitals orb; two_center_bundle_.to_LCAO_Orbitals(orb, inp.lcao_ecut, inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax); + orb_cutoff_ = orb.cutoffs(); this->set_dimension(); // setup 2d-block distribution for AO-matrix and KS wfc @@ -388,7 +390,7 @@ void LR::ESolver_LR::runner(int istep, UnitCell& cell) for (int is = 0;is < nspin;++is) { if (nspin == 2) { std::cout << "Calculating " << spin_type[is] << " excitations" << std::endl; } - hamilt::Hamilt* phamilt = new HamiltCasidaLR(xc_kernel, nspin, this->nbasis, this->nocc, this->nvirt, this->ucell, GlobalC::GridD, this->psi_ks, this->eig_ks, + hamilt::Hamilt* phamilt = new HamiltCasidaLR(xc_kernel, nspin, this->nbasis, this->nocc, this->nvirt, this->ucell, orb_cutoff_, GlobalC::GridD, this->psi_ks, this->eig_ks, #ifdef __EXX this->exx_lri, this->exx_info.info_global.hybrid_alpha, #endif diff --git a/source/module_lr/esolver_lrtd_lcao.h b/source/module_lr/esolver_lrtd_lcao.h index 9b25651073..fa7e7a35a4 100644 --- a/source/module_lr/esolver_lrtd_lcao.h +++ b/source/module_lr/esolver_lrtd_lcao.h @@ -50,6 +50,7 @@ namespace LR protected: const Input_para& input; const UnitCell& ucell; + std::vector orb_cutoff_; // not to use ElecState because 2-particle state is quite different from 1-particle state. // implement a independent one (ExcitedState) to pack physical properties if needed. diff --git a/source/module_lr/hamilt_casida.h b/source/module_lr/hamilt_casida.h index 6766975d5b..9eaeb918fb 100644 --- a/source/module_lr/hamilt_casida.h +++ b/source/module_lr/hamilt_casida.h @@ -20,6 +20,7 @@ namespace LR const int& nocc, const int& nvirt, const UnitCell& ucell_in, + const std::vector& orb_cutoff, Grid_Driver& gd_in, const psi::Psi* psi_ks_in, const ModuleBase::matrix& eig_ks, @@ -42,7 +43,7 @@ namespace LR this->ops = new OperatorLRDiag(eig_ks, pX_in, nk, nocc, nvirt); //add Hxc operator OperatorLRHxc* lr_hxc = new OperatorLRHxc(nspin, naos, nocc, nvirt, psi_ks_in, - this->DM_trans, gint_in, pot_in, ucell_in, gd_in, kv_in, pX_in, pc_in, pmat_in); + this->DM_trans, gint_in, pot_in, ucell_in, orb_cutoff, gd_in, kv_in, pX_in, pc_in, pmat_in); this->ops->add(lr_hxc); #ifdef __EXX if (xc_kernel == "hf" || xc_kernel == "hse") @@ -74,4 +75,4 @@ namespace LR /// Hxc+Exx: size=nbands, store the result of each bands for common use std::vector>> DM_trans; }; -} \ No newline at end of file +} diff --git a/source/module_lr/operator_casida/operator_lr_hxc.h b/source/module_lr/operator_casida/operator_lr_hxc.h index 9f679ba2eb..1c5ac17591 100644 --- a/source/module_lr/operator_casida/operator_lr_hxc.h +++ b/source/module_lr/operator_casida/operator_lr_hxc.h @@ -23,6 +23,7 @@ namespace LR typename TGint::type* gint_in, std::weak_ptr pot_in, const UnitCell& ucell_in, + const std::vector& orb_cutoff, Grid_Driver& gd_in, const K_Vectors& kv_in, Parallel_2D* pX_in, @@ -30,7 +31,7 @@ namespace LR Parallel_Orbitals* pmat_in) : nspin(nspin), naos(naos), nocc(nocc), nvirt(nvirt), psi_ks(psi_ks_in), DM_trans(DM_trans_in), gint(gint_in), pot(pot_in), - ucell(ucell_in), gd(gd_in), kv(kv_in), + ucell(ucell_in), orb_cutoff_(orb_cutoff), gd(gd_in), kv(kv_in), pX(pX_in), pc(pc_in), pmat(pmat_in) { ModuleBase::TITLE("OperatorLRHxc", "OperatorLRHxc"); @@ -65,8 +66,7 @@ namespace LR int iat2 = this->ucell.itia2iat(T2, I2); if (pmat->get_row_size(iat1) <= 0 || pmat->get_col_size(iat2) <= 0) { continue; } const ModuleBase::Vector3& R_index = adjs.box[ad]; - const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); - if (ucell.cal_dtau(iat1, iat2, R_index).norm() * this->ucell.lat0 >= orb.Phi[T1].getRcut() + orb.Phi[T2].getRcut()) { continue; } + if (ucell.cal_dtau(iat1, iat2, R_index).norm() * this->ucell.lat0 >= orb_cutoff_[T1] + orb_cutoff_[T2]) { continue; } hamilt::AtomPair tmp(iat1, iat2, R_index.x, R_index.y, R_index.z, pmat); hR.insert_pair(tmp); } @@ -120,6 +120,7 @@ namespace LR typename TGint::type* gint = nullptr; const UnitCell& ucell; + std::vector orb_cutoff_; Grid_Driver& gd; bool tdm_sym = false; ///< whether transition density matrix is symmetric @@ -127,4 +128,4 @@ namespace LR /// test mutable bool first_print = true; }; -} \ No newline at end of file +}