diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp index 5c4af43f40..db9e1270b7 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp @@ -82,13 +82,13 @@ template OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, HContainer*hR_in, const UnitCell& ucell_in, - const K_Vectors& kv_in, - std::vector>>>* Hexxd_in, - std::vector>>>>* Hexxc_in, + const K_Vectors& kv_in, + std::vector>>>* Hexxd_in, + std::vector>>>>* Hexxc_in, Add_Hexx_Type add_hexx_type_in, const int istep, int* two_level_step_in, - const bool restart_in) + const bool restart_in) : OperatorLCAO(hsk_in, kv_in.kvec_d, hR_in), ucell(ucell_in), kv(kv_in), @@ -105,42 +105,75 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, if (PARAM.inp.calculation == "nscf" && GlobalC::exx_info.info_global.cal_exx) { // if nscf, read HexxR first and reallocate hR according to the read-in HexxR - const std::string file_name_exx = PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); - bool all_exist = true; - for (int is=0;is std::vector { - std::ifstream ifs(file_name_exx + "_" + std::to_string(is) + ".csr"); - if (!ifs) { all_exist = false; break; } - } - if (all_exist) + std::vector file_name_list; + for (int irank=0; irank std::vector + { + std::vector file_name_list; + for (int irank=0; irank &file_name_list) -> bool + { + for (const std::string &file_name : file_name_list) + { + std::ifstream ifs(file_name); + if (!ifs.is_open()) + { return false; } + } + return true; + }; + + std::cout<<" Attention: The number of MPI processes must be strictly identical between SCF and NSCF when computing exact-exchange."<add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxd, this->hR); } + ModuleIO::read_Hexxs_csr(file_name_exx_csr, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd); + if (this->add_hexx_type == Add_Hexx_Type::R) + { reallocate_hcontainer(*Hexxd, this->hR); } } else { - ModuleIO::read_Hexxs_csr(file_name_exx, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc); - if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxc, this->hR); } + ModuleIO::read_Hexxs_csr(file_name_exx_csr, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc); + if (this->add_hexx_type == Add_Hexx_Type::R) + { reallocate_hcontainer(*Hexxc, this->hR); } } } - else + else if (check_exist(file_name_list_cereal())) { // Read HexxR in binary format (old version) - const std::string file_name_exx_cereal = PARAM.globalv.global_readin_dir + "HexxR_" + std::to_string(GlobalV::MY_RANK); + const std::string file_name_exx_cereal = PARAM.globalv.global_readin_dir + "HexxR_" + std::to_string(PARAM.globalv.myrank); + std::ifstream ifs(file_name_exx_cereal, std::ios::binary); + if (!ifs) + { ModuleBase::WARNING_QUIT("OperatorEXX", "Can't open EXX file < " + file_name_exx_cereal + " >."); } if (GlobalC::exx_info.info_ri.real_number) { ModuleIO::read_Hexxs_cereal(file_name_exx_cereal, *Hexxd); - if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxd, this->hR); } + if (this->add_hexx_type == Add_Hexx_Type::R) + { reallocate_hcontainer(*Hexxd, this->hR); } } else { ModuleIO::read_Hexxs_cereal(file_name_exx_cereal, *Hexxc); - if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxc, this->hR); } + if (this->add_hexx_type == Add_Hexx_Type::R) + { reallocate_hcontainer(*Hexxc, this->hR); } } } + else + { + ModuleBase::WARNING_QUIT("OperatorEXX", "Can't open EXX file in " + PARAM.globalv.global_readin_dir); + } this->use_cell_nearest = false; } else @@ -207,7 +240,7 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, else if (this->add_hexx_type == Add_Hexx_Type::R) { // read in Hexx(R) - const std::string restart_HR_path = PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); + const std::string restart_HR_path = PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(PARAM.globalv.myrank); bool all_exist = true; for (int is = 0; is < PARAM.inp.nspin; ++is) { @@ -227,7 +260,7 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, else { // Read HexxR in binary format (old version) - const std::string restart_HR_path_cereal = GlobalC::restart.folder + "HexxR_" + std::to_string(GlobalV::MY_RANK); + const std::string restart_HR_path_cereal = GlobalC::restart.folder + "HexxR_" + std::to_string(PARAM.globalv.myrank); if (GlobalC::exx_info.info_ri.real_number) { ModuleIO::read_Hexxs_cereal(restart_HR_path_cereal, *Hexxd); } diff --git a/source/module_io/restart_exx_csr.hpp b/source/module_io/restart_exx_csr.hpp index 6f6f3dd2ae..da6ae6131c 100644 --- a/source/module_io/restart_exx_csr.hpp +++ b/source/module_io/restart_exx_csr.hpp @@ -28,7 +28,12 @@ namespace ModuleIO { const std::vector& R = csr.getRCoordinate(iR); TC dR({ R[0], R[1], R[2] }); - Hexxs[is][iat1][{iat2, dR}] = RI::Tensor({ static_cast(ucell.atoms[ucell.iat2it[iat1]].nw), static_cast(ucell.atoms[ucell.iat2it[iat2]].nw) }); + Hexxs[is][iat1][{iat2, dR}] = RI::Tensor( + { + static_cast(ucell.atoms[ucell.iat2it[iat1]].nw), + static_cast(ucell.atoms[ucell.iat2it[iat2]].nw) + } + ); } } } @@ -44,7 +49,12 @@ namespace ModuleIO const int& npol = ucell.get_npol(); const int& i = ijv.first.first * npol; const int& j = ijv.first.second * npol; - Hexxs.at(is).at(ucell.iwt2iat[i]).at({ ucell.iwt2iat[j], { R[0], R[1], R[2] } })(ucell.iwt2iw[i] / npol, ucell.iwt2iw[j] / npol) = ijv.second; + Hexxs.at(is).at(ucell.iwt2iat[i]).at( + { + ucell.iwt2iat[j], + { R[0], R[1], R[2] } + } + )(ucell.iwt2iw[i] / npol, ucell.iwt2iw[j] / npol) = ijv.second; } } } @@ -57,6 +67,8 @@ namespace ModuleIO ModuleBase::TITLE("Exx_LRI", "read_Hexxs_cereal"); ModuleBase::timer::tick("Exx_LRI", "read_Hexxs_cereal"); std::ifstream ifs(file_name, std::ios::binary); + if(!ifs.is_open()) + { ModuleBase::WARNING_QUIT("read_Hexxs_cereal", file_name+" not found."); } cereal::BinaryInputArchive iar(ifs); iar(Hexxs); ModuleBase::timer::tick("Exx_LRI", "read_Hexxs_cereal"); diff --git a/source/module_ri/ABFs_Construct-PCA.cpp b/source/module_ri/ABFs_Construct-PCA.cpp index 933b8c9804..357348c9ee 100644 --- a/source/module_ri/ABFs_Construct-PCA.cpp +++ b/source/module_ri/ABFs_Construct-PCA.cpp @@ -32,14 +32,14 @@ namespace PCA dsyev_(&jobz, &uplo, &nr, a.ptr(), &nc, w, work.data(), &lwork, &info); } - RI::Tensor get_sub_matrix( + RI::Tensor get_sub_matrix( const RI::Tensor & m, // size: (lcaos, lcaos, abfs) const std::size_t & T, const std::size_t & L, const ModuleBase::Element_Basis_Index::Range & range, const ModuleBase::Element_Basis_Index::IndexLNM & index ) { - ModuleBase::TITLE("ABFs_Construct::PCA::get_sub_matrix"); + 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) { @@ -74,12 +74,12 @@ namespace PCA std::vector, RI::Tensor>>> cal_PCA( const UnitCell &ucell, const LCAO_Orbitals& orb, - const std::vector>> &lcaos, + const std::vector>> &lcaos, const std::vector>> &abfs, const double kmesh_times ) { ModuleBase::TITLE("ABFs_Construct::PCA::cal_PCA"); - + const ModuleBase::Element_Basis_Index::Range range_lcaos = Exx_Abfs::Abfs_Index::construct_range( lcaos ); const ModuleBase::Element_Basis_Index::IndexLNM @@ -107,27 +107,27 @@ namespace PCA m_abfslcaos_lcaos.init_radial_table(delta_R); GlobalC::exx_info.info_ri.abfs_Lmax = Lmax_bak; - + std::vector,RI::Tensor>>> eig(abfs.size()); for( std::size_t T=0; T!=abfs.size(); ++T ) { - const RI::Tensor A = m_abfslcaos_lcaos.cal_overlap_matrix( - T, - T, + const RI::Tensor A = m_abfslcaos_lcaos.cal_overlap_matrix( + T, + T, + ModuleBase::Vector3{0,0,0}, ModuleBase::Vector3{0,0,0}, - ModuleBase::Vector3{0,0,0}, - index_abfs, + index_abfs, index_lcaos, index_lcaos, Matrix_Orbs21::Matrix_Order::A2BA1); - + eig[T].resize(abfs[T].size()); for( std::size_t L=0; L!=abfs[T].size(); ++L ) { const RI::Tensor A_sub = get_sub_matrix( A, T, L, range_abfs, index_abfs ); RI::Tensor mm = A_sub.transpose() * A_sub; std::vector eig_value(mm.shape[0]); - + int info=1; tensor_dsyev('V', 'L', mm, eig_value.data(), info); @@ -158,7 +158,7 @@ namespace PCA eig[T][L] = std::make_pair( eig_value, mm ); } } - + return eig; } diff --git a/source/module_ri/ABFs_Construct-PCA.h b/source/module_ri/ABFs_Construct-PCA.h index af5390ab9f..6ebf6b5b7a 100644 --- a/source/module_ri/ABFs_Construct-PCA.h +++ b/source/module_ri/ABFs_Construct-PCA.h @@ -15,10 +15,10 @@ namespace ABFs_Construct { namespace PCA { - extern std::vector,RI::Tensor>>> cal_PCA( + extern std::vector,RI::Tensor>>> cal_PCA( const UnitCell& ucell, const LCAO_Orbitals &orb, - const std::vector>> &lcaos, + const std::vector>> &lcaos, const std::vector>> &abfs, // abfs must be orthonormal const double kmesh_times ); } diff --git a/source/module_ri/Exx_LRI_interface.h b/source/module_ri/Exx_LRI_interface.h index 2b6fe01b0c..77aaa6c47b 100644 --- a/source/module_ri/Exx_LRI_interface.h +++ b/source/module_ri/Exx_LRI_interface.h @@ -42,9 +42,9 @@ class Exx_LRI_Interface } Exx_LRI_Interface() = delete; - /// read and write Hexxs using cereal - void write_Hexxs_cereal(const std::string& file_name) const; - void read_Hexxs_cereal(const std::string& file_name); + ///// read and write Hexxs using cereal + //void write_Hexxs_cereal(const std::string& file_name) const; + //void read_Hexxs_cereal(const std::string& file_name); std::vector>>>& get_Hexxs() const { return this->exx_ptr->Hexxs; } double &get_Eexx() const { return this->exx_ptr->Eexx; } diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index 0a3bb72419..724835d997 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -18,6 +18,7 @@ #include #include +/* template void Exx_LRI_Interface::write_Hexxs_cereal(const std::string& file_name) const { @@ -32,13 +33,17 @@ void Exx_LRI_Interface::write_Hexxs_cereal(const std::string& file_nam template void Exx_LRI_Interface::read_Hexxs_cereal(const std::string& file_name) { - ModuleBase::TITLE("Exx_LRI", "read_Hexxs_cereal"); - ModuleBase::timer::tick("Exx_LRI", "read_Hexxs_cereal"); - std::ifstream ifs(file_name + "_" + std::to_string(GlobalV::MY_RANK), std::ofstream::binary); + ModuleBase::TITLE("Exx_LRI_Interface", "read_Hexxs_cereal"); + ModuleBase::timer::tick("Exx_LRI_Interface", "read_Hexxs_cereal"); + const std::string file_name_rank = file_name + "_" + std::to_string(GlobalV::MY_RANK); + std::ifstream ifs(file_name_rank, std::ofstream::binary); + if(!ifs.is_open()) + { ModuleBase::WARNING_QUIT("Exx_LRI_Interface", file_name_rank+" not found."); } cereal::BinaryInputArchive iar(ifs); iar(this->exx_ptr->Hexxs); ModuleBase::timer::tick("Exx_LRI", "read_Hexxs_cereal"); } +*/ template void Exx_LRI_Interface::init(const MPI_Comm &mpi_comm, @@ -115,7 +120,8 @@ void Exx_LRI_Interface::exx_before_all_runners(const K_Vectors& kv, co if (this->exx_spacegroup_symmetry) { const std::array& period = RI_Util::get_Born_vonKarmen_period(kv); - this->symrot_.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, + this->symrot_.find_irreducible_sector( + ucell.symm, ucell.atoms, ucell.st, RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat); // this->symrot_.set_Cs_rotation(this->exx_ptr->get_abfs_nchis()); this->symrot_.cal_Ms(kv, ucell, pv); @@ -209,8 +215,19 @@ void Exx_LRI_Interface::exx_eachiterinit(const int istep, { this->mix_DMk_2D.mix(dm_in.get_DMK_vector(), flag_restart); } const std::vector>>> Ds = PARAM.globalv.gamma_only_local - ? RI_2D_Comm::split_m2D_ktoR(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm_in.get_paraV_pointer(), PARAM.inp.nspin) - : RI_2D_Comm::split_m2D_ktoR(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm_in.get_paraV_pointer(), PARAM.inp.nspin, this->exx_spacegroup_symmetry); + ? RI_2D_Comm::split_m2D_ktoR( + ucell, + *this->exx_ptr->p_kv, + this->mix_DMk_2D.get_DMk_gamma_out(), + *dm_in.get_paraV_pointer(), + PARAM.inp.nspin) + : RI_2D_Comm::split_m2D_ktoR( + ucell, + *this->exx_ptr->p_kv, + this->mix_DMk_2D.get_DMk_k_out(), + *dm_in.get_paraV_pointer(), + PARAM.inp.nspin, + this->exx_spacegroup_symmetry); if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { this->cal_exx_elec(Ds, ucell,*dm_in.get_paraV_pointer(), &this->symrot_); } @@ -240,11 +257,10 @@ void Exx_LRI_Interface::exx_hamilt2density(elecstate::ElecState& elec, { if (GlobalV::MY_RANK == 0) { - try { GlobalC::restart.load_disk("Eexx", 0, 1, &this->exx_ptr->Eexx); } + try + { GlobalC::restart.load_disk("Eexx", 0, 1, &this->exx_ptr->Eexx); } catch (const std::exception& e) - { - std::cout << "WARNING: Cannot read Eexx from disk, the energy of the 1st loop will be wrong, sbut it does not influence the subsequent loops." << std::endl; - } + { std::cout << "WARNING: Cannot read Eexx from disk, the energy of the 1st loop will be wrong, sbut it does not influence the subsequent loops." << std::endl; } } Parallel_Common::bcast_double(this->exx_ptr->Eexx); this->exx_ptr->Eexx /= GlobalC::exx_info.info_global.hybrid_alpha; diff --git a/source/module_ri/exx_abfs-abfs_index.cpp b/source/module_ri/exx_abfs-abfs_index.cpp index 2d3966c869..24b45fbc15 100644 --- a/source/module_ri/exx_abfs-abfs_index.cpp +++ b/source/module_ri/exx_abfs-abfs_index.cpp @@ -29,9 +29,9 @@ ModuleBase::Element_Basis_Index::Range range[T].resize( orb[T].size() ); for( size_t L=0; L!=range[T].size(); ++L ) { - range[T][L].N = orb[T][L].size(); - range[T][L].M = 2*L+1; - } + range[T][L].N = orb[T][L].size(); + range[T][L].M = 2*L+1; + } } return range; } diff --git a/source/module_ri/exx_abfs-abfs_index.h b/source/module_ri/exx_abfs-abfs_index.h index 09cc493f2e..c65f87c8b1 100644 --- a/source/module_ri/exx_abfs-abfs_index.h +++ b/source/module_ri/exx_abfs-abfs_index.h @@ -12,7 +12,7 @@ class LCAO_Orbitals; class Exx_Abfs::Abfs_Index { public: - static ModuleBase::Element_Basis_Index::Range construct_range( const LCAO_Orbitals &orb ); + static ModuleBase::Element_Basis_Index::Range construct_range( const LCAO_Orbitals &orb ); static ModuleBase::Element_Basis_Index::Range construct_range( const std::vector>> &orb ); }; diff --git a/source/module_ri/exx_abfs-construct_orbs.h b/source/module_ri/exx_abfs-construct_orbs.h index cd33b19d3d..2908ad628f 100644 --- a/source/module_ri/exx_abfs-construct_orbs.h +++ b/source/module_ri/exx_abfs-construct_orbs.h @@ -12,33 +12,33 @@ class Exx_Abfs::Construct_Orbs { public: - static std::vector>> change_orbs( + static std::vector>> change_orbs( const LCAO_Orbitals &orb_in, const double kmesh_times ); - static std::vector>> change_orbs( + static std::vector>> change_orbs( const std::vector>> &orb_in, const double kmesh_times ); - static std::vector>> abfs_same_atom( + static std::vector>> abfs_same_atom( const UnitCell &ucell, const LCAO_Orbitals& orb, const std::vector>> &lcaos, const double kmesh_times_mot, const double times_threshold=0); - + static void print_orbs_size( const UnitCell& ucell, const std::vector>> &orbs, - std::ostream &os); - + std::ostream &os); + private: - static std::vector>>> psi_mult_psi( + static std::vector>>> psi_mult_psi( const std::vector>> &lcaos ); - - static std::vector>>> psir_mult_psir( + + static std::vector>>> psir_mult_psir( const std::vector>> &lcaos ); - - static std::vector>>> orth( + + static std::vector>>> orth( const std::vector>>> &psis, const std::vector>> &lcaos, const double norm_threshold = std::numeric_limits::min() ); @@ -50,16 +50,16 @@ class Exx_Abfs::Construct_Orbs const std::vector>> &orbs, const double kmesh_times_mot, const double times_threshold ); - - static std::vector>>> div_r( + + static std::vector>>> div_r( const std::vector>>> &psirs, - const std::vector &r_radial ); - + const std::vector &r_radial ); + static std::vector>> orbital( const std::vector>>> &psis, const std::vector>> &orbs_info, - const double kmesh_times); - + const double kmesh_times); + static std::vector>>> get_psi( const std::vector>> &orbs ); }; diff --git a/source/module_ri/exx_abfs-io.cpp b/source/module_ri/exx_abfs-io.cpp index 938d2b0238..edad3fa348 100644 --- a/source/module_ri/exx_abfs-io.cpp +++ b/source/module_ri/exx_abfs-io.cpp @@ -18,26 +18,26 @@ std::vector>> Exx_Abfs::IO::constr { std::vector>> abfs( files_abfs.size() ); for( size_t T=0; T!=files_abfs.size(); ++T ) - abfs[T] = construct_abfs_T( + abfs[T] = construct_abfs_T( files_abfs[T], T, static_cast(orbs.get_kmesh() * kmesh_times) | 1, // Nk must be odd // orbs.get_dk() / kmesh_times, orbs.get_dk(), // Peize Lin change 2017-04-16 - orbs.get_dr_uniform() ); - + orbs.get_dr_uniform() ); + return abfs; } -std::vector>> Exx_Abfs::IO::construct_abfs( +std::vector>> Exx_Abfs::IO::construct_abfs( const std::vector>> & abfs_pre, const LCAO_Orbitals &orbs, const std::vector &files_abfs, const double kmesh_times ) { - std::vector>> + std::vector>> &&abfs = construct_abfs( orbs, files_abfs, kmesh_times ); - + assert( abfs.size() == abfs_pre.size() ); for( size_t T=0; T!=abfs.size(); ++T ) { @@ -47,12 +47,12 @@ std::vector>> Exx_Abfs::IO::constr { abfs[T][L].insert( abfs[T][L].begin(), abfs_pre[T][L].begin(), abfs_pre[T][L].end() ); } - } - + } + return abfs; } -std::vector> Exx_Abfs::IO::construct_abfs_T( +std::vector> Exx_Abfs::IO::construct_abfs_T( const std::string & file_name, const int &T, const int &nk, @@ -65,20 +65,20 @@ std::vector> Exx_Abfs::IO::construct_abfs_T( size_t meshr; double dr; std::map>> psis; - + /*---------------------- 1.read abfs ----------------------*/ std::string word; - + std::ifstream ifs( file_name.c_str() ); if(!ifs) - throw std::runtime_error(" Can't find the abfs ORBITAL file."); - + throw std::runtime_error(" Can't find the abfs ORBITAL file."); + while( ifs.good() ) { ifs >> word; - + if( "Element"==word ) { ModuleBase::GlobalFunc::READ_VALUE( ifs, label ); @@ -86,7 +86,7 @@ std::vector> Exx_Abfs::IO::construct_abfs_T( else if ( "Lmax"==word ) { ModuleBase::GlobalFunc::READ_VALUE( ifs, L_size ); - + if( L_size>=9 ) { std::stringstream ss; @@ -133,12 +133,12 @@ std::vector> Exx_Abfs::IO::construct_abfs_T( else if ( "END"==word ) { break; - } + } } - + ModuleBase::CHECK_NAME(ifs, "Mesh"); ifs >> meshr; - + ModuleBase::CHECK_NAME(ifs, "dr"); ifs >> dr; @@ -149,21 +149,21 @@ std::vector> Exx_Abfs::IO::construct_abfs_T( { std::string s_L, s_N; ifs >> s_L >> s_N; - + size_t T,L,N; ifs >> T >> L >> N; - + psis[L][N].resize(meshr); for(int ir=0; ir!=meshr; ir++) { ifs >> psis[L][N][ir]; - } + } } } ifs.close(); - + /*---------------------- 2.check L,N orbital ----------------------*/ @@ -183,12 +183,12 @@ std::vector> Exx_Abfs::IO::construct_abfs_T( throw std::domain_error(ss.str()); } - + /*---------------------- 3.rab, radial - ----------------------*/ + ----------------------*/ if(meshr%2==0) ++meshr; - + std::vector rab(meshr); std::vector radial(meshr); for( int ir=0; ir!=meshr; ++ir ) @@ -197,12 +197,12 @@ std::vector> Exx_Abfs::IO::construct_abfs_T( radial[ir] = ir*dr; //mohan 2010-04-19 } - + /*---------------------- 4.normalize psi ----------------------*/ for( size_t L=0; L<=L_size; ++L ) - { + { for( size_t N=0; N!=N_size[L]; ++N ) { std::vector psir(meshr); @@ -213,7 +213,7 @@ std::vector> Exx_Abfs::IO::construct_abfs_T( psir[ir] = psis[L][N][ir] * radial[ir]; inner[ir] = psir[ir] * psir[ir]; } - double unit = 0.0; + double unit = 0.0; ModuleBase::Integral::Simpson_Integral(meshr, ModuleBase::GlobalFunc::VECTOR_TO_PTR(inner), ModuleBase::GlobalFunc::VECTOR_TO_PTR(rab), unit); for( int ir=0; ir!=meshr; ++ir ) { @@ -222,10 +222,10 @@ std::vector> Exx_Abfs::IO::construct_abfs_T( } } - + /*---------------------- 5.construct abfs - ----------------------*/ + ----------------------*/ std::vector> abfs_T; abfs_T.resize(L_size+1); @@ -248,9 +248,9 @@ std::vector> Exx_Abfs::IO::construct_abfs_T( dk, dr_uniform, false, - true, PARAM.inp.cal_force); + true, PARAM.inp.cal_force); } } - + return abfs_T; } diff --git a/source/module_ri/exx_abfs-io.h b/source/module_ri/exx_abfs-io.h index 15dd9303a3..7fdf516b81 100644 --- a/source/module_ri/exx_abfs-io.h +++ b/source/module_ri/exx_abfs-io.h @@ -18,14 +18,14 @@ class LCAO_Orbitals; class Exx_Abfs::IO { public: - - static std::vector>> construct_abfs( + + static std::vector>> construct_abfs( const LCAO_Orbitals &orbs, const std::vector &files_abfs, - const double kmesh_times=1 ); // close dK, keep Kcut - - static std::vector>> construct_abfs( - const std::vector>> &abfs_pre, + const double kmesh_times=1 ); // close dK, keep Kcut + + static std::vector>> construct_abfs( + const std::vector>> &abfs_pre, const LCAO_Orbitals &orbs, const std::vector &files_abfs, const double kmesh_times=1 ); // close dK, keep Kcut diff --git a/source/module_ri/exx_abfs-jle.cpp b/source/module_ri/exx_abfs-jle.cpp index f9086fe777..d7957bc09b 100644 --- a/source/module_ri/exx_abfs-jle.cpp +++ b/source/module_ri/exx_abfs-jle.cpp @@ -10,9 +10,9 @@ bool Exx_Abfs::Jle::generate_matrix = false; int Exx_Abfs::Jle::Lmax = 2; double Exx_Abfs::Jle::Ecut_exx = 60; -double Exx_Abfs::Jle::tolerence = 1.0e-12; +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 UnitCell& ucell, const LCAO_Orbitals& orb) { @@ -23,7 +23,7 @@ void Exx_Abfs::Jle::init_jle(const double kmesh_times, jle[T].resize( Lmax+1 ); for (int L=0; L <= Lmax ; ++L) { - const size_t ecut_number + const size_t ecut_number = static_cast( sqrt( Ecut_exx ) * orb.Phi[T].getRcut() / ModuleBase::PI ); // Rydberg Unit. jle[T][L].resize( ecut_number ); @@ -35,10 +35,10 @@ void Exx_Abfs::Jle::init_jle(const double kmesh_times, { std::vector jle_r( orb.Phi[T].PhiLN(0,0).getNr() ); ModuleBase::Sphbes::Spherical_Bessel( - orb.Phi[T].PhiLN(0,0).getNr(), - orb.Phi[T].PhiLN(0,0).getRadial(), - en[E], - L, + 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( orb.Phi[T].PhiLN(0,0).getLabel(), diff --git a/source/module_ri/exx_abfs-jle.h b/source/module_ri/exx_abfs-jle.h index d41078248e..dec3178b89 100644 --- a/source/module_ri/exx_abfs-jle.h +++ b/source/module_ri/exx_abfs-jle.h @@ -10,13 +10,13 @@ class Exx_Abfs::Jle { public: - + std::vector< std::vector< std::vector< Numerical_Orbital_Lm>>> jle; - void init_jle(const double kmesh_times, + void init_jle(const double kmesh_times, const UnitCell& ucell, const LCAO_Orbitals& orb); diff --git a/source/module_ri/exx_abfs.h b/source/module_ri/exx_abfs.h index 845eb3d81b..e684807874 100644 --- a/source/module_ri/exx_abfs.h +++ b/source/module_ri/exx_abfs.h @@ -20,7 +20,7 @@ class Exx_Abfs class IO; class Construct_Orbs; class PCA; - + int rmesh_times = 5; // Peize Lin test int kmesh_times = 1; // Peize Lin test diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index 37cd742d7a..ec08cb09e3 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -29,12 +29,12 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co abfs = Exx_Abfs::Construct_Orbs::abfs_same_atom(ucell,orb, lcaos, this->kmesh_times, GlobalC::exx_info.info_ri.pca_threshold ); // ofs_mpi<<"memory:\t"<kmesh_times, ucell , orb); // ofs_mpi<<"memory:\t"<(abfs[T].size())-1 ); @@ -70,7 +70,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co // ofs_mpi<<"memory:\t"< - const auto ms_lcaoslcaos_lcaoslcaos = [&]() -> std::map>>>> + const auto ms_lcaoslcaos_lcaoslcaos = [&]() -> std::map>>>> { Matrix_Orbs22 m_lcaoslcaos_lcaoslcaos; m_lcaoslcaos_lcaoslcaos.init( 1, ucell,orb, this->kmesh_times, 1 ); @@ -82,7 +82,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co #endif return m_lcaoslcaos_lcaoslcaos.cal_overlap_matrix_all(ucell,index_lcaos, index_lcaos, index_lcaos, index_lcaos); }(); - + // ofs_mpi<<"memory:\t"< @@ -309,10 +309,10 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co } // m_big - m_left * m_middle * m_right.T -RI::Tensor Exx_Opt_Orb::cal_proj( - const RI::Tensor & m_big, - const std::vector> & m_left, - const std::vector>> & m_middle, +RI::Tensor Exx_Opt_Orb::cal_proj( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, const std::vector> & m_right ) const { ModuleBase::TITLE("Exx_Opt_Orb::cal_proj"); diff --git a/source/module_ri/exx_opt_orb.h b/source/module_ri/exx_opt_orb.h index 2b68dfa40a..094d444eeb 100644 --- a/source/module_ri/exx_opt_orb.h +++ b/source/module_ri/exx_opt_orb.h @@ -15,27 +15,27 @@ class Exx_Opt_Orb public: void generate_matrix(const K_Vectors &kv, const UnitCell& ucell, const LCAO_Orbitals& orb) const; private: - std::vector>> cal_I( - const std::map>>>> &ms, + std::vector>> cal_I( + const std::map>>>> &ms, const size_t TA, const size_t IA, const size_t TB, const size_t IB ) const; - RI::Tensor cal_proj( - const RI::Tensor & m_big, - const std::vector> & m_left, - const std::vector>> & m_middle, + RI::Tensor cal_proj( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, const std::vector> & m_right ) const; void print_matrix( const UnitCell& ucell, const K_Vectors &kv, const std::string& file_name, - const std::vector> &matrix_Q, + const std::vector> &matrix_Q, 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::Range &range_jles, const ModuleBase::Element_Basis_Index::IndexLNM &index_jles) const; std::map>> get_radial_R(const UnitCell& ucell) const; - + int kmesh_times = 4; }; #endif