diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index ace9097bb0..c9b3b69601 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -12,7 +12,6 @@ #include "module_ri/exx_abfs-construct_orbs.h" #include "module_ri/exx_abfs-io.h" #include "module_ri/conv_coulomb_pot_k.h" -#include "module_ri/conv_coulomb_pot_k-template.h" #include "module_base/tool_title.h" #include "module_base/timer.h" #include "module_ri/serialization_cereal.h" @@ -71,14 +70,19 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in) case Conv_Coulomb_Pot_K::Ccp_Type::Ccp: return {}; case Conv_Coulomb_Pot_K::Ccp_Type::Hf: - return {}; + { + // 4/3 * pi * Rcut^3 = V_{supercell} = V_{unitcell} * Nk + const int nspin0 = (GlobalV::NSPIN==2) ? 2 : 1; + const double hf_Rcut = std::pow(0.75 * this->p_kv->nkstot_full/nspin0 * GlobalC::ucell.omega / (ModuleBase::PI), 1.0/3.0); + return {{"hf_Rcut", hf_Rcut}}; + } case Conv_Coulomb_Pot_K::Ccp_Type::Hse: return {{"hse_omega", this->info.hse_omega}}; default: throw std::domain_error(std::string(__FILE__)+" line "+std::to_string(__LINE__)); break; } }; - this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, this->info.ccp_type, get_ccp_parameter(), this->info.ccp_rmesh_times, this->p_kv->nkstot_full); + this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, this->info.ccp_type, get_ccp_parameter(), this->info.ccp_rmesh_times); for( size_t T=0; T!=this->abfs.size(); ++T ) diff --git a/source/module_ri/conv_coulomb_pot_k-template.h b/source/module_ri/conv_coulomb_pot_k-template.h deleted file mode 100644 index 9a3d245286..0000000000 --- a/source/module_ri/conv_coulomb_pot_k-template.h +++ /dev/null @@ -1,51 +0,0 @@ -#ifndef CONV_COULOMB_POT_K_TEMPLATE_H -#define CONV_COULOMB_POT_K_TEMPLATE_H - -#include "conv_coulomb_pot_k.h" -#include -#include - -#include "../module_ri/test_code/exx_abfs-construct_orbs-test.h" - - -template< typename T > -T Conv_Coulomb_Pot_K::cal_orbs_ccp( - const T & orbs, - const Ccp_Type &ccp_type, - const std::map ¶meter, - const double rmesh_times, - const int& nks) -{ - T orbs_ccp(orbs.size()); - for( size_t i=0; i!=orbs.size(); ++i ) - orbs_ccp[i] = cal_orbs_ccp(orbs[i], ccp_type, parameter, rmesh_times, nks ); - return orbs_ccp; -} - -extern template -Numerical_Orbital_Lm Conv_Coulomb_Pot_K::cal_orbs_ccp( - const Numerical_Orbital_Lm & orbs, - const Ccp_Type &ccp_type, - const std::map ¶meter, - const double rmesh_times, - const int& nks); - - - -template< typename T > -double Conv_Coulomb_Pot_K::get_rmesh_proportion( - const T & orbs, - const double psi_threshold) -{ - double rmesh_proportion=0; - for( const auto &orb : orbs ) - rmesh_proportion = std::max(rmesh_proportion, get_rmesh_proportion(orb,psi_threshold)); - return rmesh_proportion; -} - -extern template -double Conv_Coulomb_Pot_K::get_rmesh_proportion( - const Numerical_Orbital_Lm & orbs, - const double psi_threshold); - -#endif \ No newline at end of file diff --git a/source/module_ri/conv_coulomb_pot_k.cpp b/source/module_ri/conv_coulomb_pot_k.cpp index 9f573509ee..62dd582a44 100644 --- a/source/module_ri/conv_coulomb_pot_k.cpp +++ b/source/module_ri/conv_coulomb_pot_k.cpp @@ -2,104 +2,109 @@ #include "../module_base/constants.h" #include "../module_basis/module_ao/ORB_atomic_lm.h" #include "../module_hamilt_pw/hamilt_pwdft/global.h" -std::vector Conv_Coulomb_Pot_K::cal_psi_ccp( const std::vector & psif ) + +namespace Conv_Coulomb_Pot_K { - std::vector psik2_ccp(psif.size()); - for( size_t ik=0; ik Conv_Coulomb_Pot_K::cal_psi_hf(const int& nks, const std::vector &psif, - const std::vector &k_radial, - const double omega = 0) -{ - const int nspin0 = (GlobalV::NSPIN==2) ? 2 : 1; - const double Rc = std::pow(0.75 * nks/nspin0 * GlobalC::ucell.omega / (ModuleBase::PI), 1.0/3.0); - std::vector psik2_ccp(psif.size()); - for (size_t ik = 0; ik < psif.size(); ++ik) - psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * (1 - std::cos(k_radial[ik] * Rc)); - return psik2_ccp; -} + std::vector cal_psi_ccp( + const std::vector & psif) + { + std::vector psik2_ccp(psif.size()); + for( size_t ik=0; ik cal_psi_hf( + const std::vector &psif, + const std::vector &k_radial, + const double hf_Rcut) + { + std::vector psik2_ccp(psif.size()); + for (size_t ik = 0; ik < psif.size(); ++ik) + psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * (1 - std::cos(k_radial[ik] * hf_Rcut)); + return psik2_ccp; + } -std::vector Conv_Coulomb_Pot_K::cal_psi_hse( - const std::vector & psif, - const std::vector & k_radial, - const double omega) -{ - std::vector psik2_ccp(psif.size()); - for( size_t ik=0; ik cal_psi_hse( + const std::vector & psif, + const std::vector & k_radial, + const double hse_omega) + { + std::vector psik2_ccp(psif.size()); + for( size_t ik=0; ik -Numerical_Orbital_Lm Conv_Coulomb_Pot_K::cal_orbs_ccp( - const Numerical_Orbital_Lm &orbs, - const Ccp_Type &ccp_type, - const std::map ¶meter, - const double rmesh_times, - const int& nks) -{ - std::vector psik2_ccp; - switch(ccp_type) + + template<> + Numerical_Orbital_Lm cal_orbs_ccp( + const Numerical_Orbital_Lm &orbs, + const Ccp_Type &ccp_type, + const std::map ¶meter, + const double rmesh_times) { - case Ccp_Type::Ccp: - psik2_ccp = cal_psi_ccp( orbs.get_psif() ); break; - case Ccp_Type::Hf: - psik2_ccp = cal_psi_hf(nks, orbs.get_psif(), orbs.get_k_radial()); break; - case Ccp_Type::Hse: - psik2_ccp = cal_psi_hse( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega") ); break; - default: - throw( ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__) ); break; - } + std::vector psik2_ccp; + switch(ccp_type) + { + case Ccp_Type::Ccp: + psik2_ccp = cal_psi_ccp( orbs.get_psif() ); break; + case Ccp_Type::Hf: + psik2_ccp = cal_psi_hf( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hf_Rcut")); break; + case Ccp_Type::Hse: + psik2_ccp = cal_psi_hse( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega") ); break; + default: + throw( ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__) ); break; + } - const double dr = orbs.get_rab().back(); - const int Nr = (static_cast(orbs.getNr()*rmesh_times)) | 1; - std::vector rab(Nr); - for( size_t ir=0; ir r_radial(Nr); - for( size_t ir=0; ir(orbs.getNr()*rmesh_times)) | 1; + std::vector rab(Nr); + for( size_t ir=0; ir r_radial(Nr); + for( size_t ir=0; ir -double Conv_Coulomb_Pot_K::get_rmesh_proportion( - const Numerical_Orbital_Lm &orbs, - const double psi_threshold) -{ - for(int ir=orbs.getNr()-1; ir>=0; --ir) + Numerical_Orbital_Lm orbs_ccp; + orbs_ccp.set_orbital_info( + orbs.getLabel(), + orbs.getType(), + orbs.getL(), + orbs.getChi(), + Nr, + ModuleBase::GlobalFunc::VECTOR_TO_PTR(rab), + ModuleBase::GlobalFunc::VECTOR_TO_PTR(r_radial), + Numerical_Orbital_Lm::Psi_Type::Psik2, + ModuleBase::GlobalFunc::VECTOR_TO_PTR(psik2_ccp), + orbs.getNk(), + orbs.getDk(), + orbs.getDruniform(), + false, + true, GlobalV::CAL_FORCE); + return orbs_ccp; + } + + template<> + double get_rmesh_proportion( + const Numerical_Orbital_Lm &orbs, + const double psi_threshold) { - if(std::abs(orbs.getPsi(ir))>=psi_threshold) - return static_cast(ir)/orbs.getNr(); + for(int ir=orbs.getNr()-1; ir>=0; --ir) + { + if(std::abs(orbs.getPsi(ir))>=psi_threshold) + return static_cast(ir)/orbs.getNr(); + } + return 0.0; } - return 0.0; + } diff --git a/source/module_ri/conv_coulomb_pot_k.h b/source/module_ri/conv_coulomb_pot_k.h index 9adec9d915..d464a53f91 100644 --- a/source/module_ri/conv_coulomb_pot_k.h +++ b/source/module_ri/conv_coulomb_pot_k.h @@ -5,40 +5,37 @@ #include #include -class Conv_Coulomb_Pot_K +namespace Conv_Coulomb_Pot_K { -public: + enum class Ccp_Type{ // parameter: + Ccp, // + Hf, // "hf_Rcut" + Hse}; // "hse_omega" - enum class Ccp_Type{ // parameter: - Ccp, // - Hf, // - Hse}; // "hse_omega" - - template static T cal_orbs_ccp( + template T cal_orbs_ccp( const T &orbs, const Ccp_Type &ccp_type, const std::map ¶meter, - const double rmesh_times, - const int& nks); - -private: - - template< typename T > static double get_rmesh_proportion( + const double rmesh_times); + + //private: + template< typename T > double get_rmesh_proportion( const T &orbs, const double psi_threshold); - -private: - static std::vector cal_psi_ccp( const std::vector & psif ); - - static std::vector cal_psi_hf(const int& nks, const std::vector &psif, - const std::vector &k_radial, - const double omega); - - static std::vector cal_psi_hse( + //private: + std::vector cal_psi_ccp( + const std::vector & psif); + std::vector cal_psi_hf( + const std::vector &psif, + const std::vector &k_radial, + const double hf_Rcut); + std::vector cal_psi_hse( const std::vector & psif, const std::vector & k_radial, - const double omega); -}; + const double hse_omega); +} + +#include "conv_coulomb_pot_k.hpp" #endif \ No newline at end of file diff --git a/source/module_ri/conv_coulomb_pot_k.hpp b/source/module_ri/conv_coulomb_pot_k.hpp new file mode 100644 index 0000000000..5ca3abe5c8 --- /dev/null +++ b/source/module_ri/conv_coulomb_pot_k.hpp @@ -0,0 +1,37 @@ +#ifndef CONV_COULOMB_POT_K_HPP +#define CONV_COULOMB_POT_K_HPP + +#include "conv_coulomb_pot_k.h" +#include +#include + +namespace Conv_Coulomb_Pot_K +{ + + template< typename T > + std::vector cal_orbs_ccp( + const std::vector & orbs, + const Ccp_Type &ccp_type, + const std::map ¶meter, + const double rmesh_times) + { + std::vector orbs_ccp(orbs.size()); + for( size_t i=0; i!=orbs.size(); ++i ) + orbs_ccp[i] = cal_orbs_ccp(orbs[i], ccp_type, parameter, rmesh_times); + return orbs_ccp; + } + + template< typename T > + double get_rmesh_proportion( + const std::vector & orbs, + const double psi_threshold) + { + double rmesh_proportion=0; + for( const auto &orb : orbs ) + rmesh_proportion = std::max(rmesh_proportion, get_rmesh_proportion(orb,psi_threshold)); + return rmesh_proportion; + } + +} + +#endif \ No newline at end of file