diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 40abf00c3b..ec765febe2 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -512,6 +512,7 @@ These variables are used to control general system parameters. - get_wf: obtain real space wave functions (for LCAO basis only). See `out_wfc_norm` and `out_wfc_re_im` for more information - get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is `SR.csr` with file format being the same as that generated by [out_mat_hs2](#out_mat_hs2). Note: in the 3.10-LTS version, the command was named `get_S` - gen_bessel: generates projectors, i.e., a series of Bessel functions, for the DeePKS method (for LCAO basis only); see also keywords `bessel_descriptor_lmax`, `bessel_descriptor_rcut` and `bessel_descriptor_tolerence`. A file named `jle.orb` will be generated which contains the projectors. An example is provided in examples/H2O-deepks-pw + - gen_opt_abfs: generate opt-ABFs as discussed in this [article](https://pubs.acs.org/doi/abs/10.1021/acs.jpclett.0c00481). - test_memory: obtain a rough estimation of memory consuption for the calculation - test_neighbour: obtain information of neighboring atoms (for LCAO basis only), please specify a positive [search_radius](#search_radius) manually - **Default**: scf @@ -540,7 +541,7 @@ These variables are used to control general system parameters. - **Default**: - 0: - if [calculation](#calculation)==md/nscf/get_pchg/get_wf/get_s or [gamma_only](#gamma_only)==True; - - If ([dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb or [rpa](#rpa)==True). + - If ([dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0 or [rpa](#rpa)==True). - If [efield_flag](#efield_flag)==1 - 1: else @@ -763,7 +764,7 @@ These variables are used to control parameters related to input files. - **Availability**: Used only when numerical atomic orbitals are employed as basis set. - **Description**: If [restart_save](#restart_save) is set to true and an electronic iteration is finished, calculations can be restarted from the charge density file, which are saved in the former calculation. Please ensure [read_file_dir](#read_file_dir) is correct, and the charge density file exist. - If EXX(exact exchange) is calculated (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*), the Hexx(R) files in the same folder for each processor will also be read. + If EXX(exact exchange) is calculated (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0* or *[rpa](#rpa)==True*), the Hexx(R) files in the same folder for each processor will also be read. - **Default**: False ### wannier_card @@ -1057,7 +1058,6 @@ calculations. Furthermore, the old INPUT parameter exx_hybrid_type for hybrid functionals has been absorbed into dft_functional. Options are `hf` (pure Hartree-Fock), `pbe0`(PBE0), `hse` (Note: in order to use HSE functional, LIBXC is required). Note also that HSE has been tested while PBE0 has NOT been fully tested yet, and the maximum CPU cores for running exx in parallel is $N(N+1)/2$, with N being the number of atoms. - If set to `opt_orb`, the program will not perform hybrid functional calculation. Instead, it is going to generate opt-ABFs as discussed in this [article](https://pubs.acs.org/doi/abs/10.1021/acs.jpclett.0c00481). - **Default**: Used the same as DFT functional as specified in the pseudopotential files. ### xc_temperature @@ -1645,7 +1645,7 @@ These variables are used to control the output of properties. --- The circle order of the charge density on real space grids is: x is the outer loop, then y and finally z (z is moving fastest). - In EXX(exact exchange) calculations, (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*), the Hexx(R) files will be output in the folder `OUT.${suffix}` too, which can be read in NSCF calculation. + In EXX(exact exchange) calculations, (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0* or *[rpa](#rpa)==True*), the Hexx(R) files will be output in the folder `OUT.${suffix}` too, which can be read in NSCF calculation. In molecular dynamics simulations, the output frequency is controlled by [out_interval](#out_interval). - **Default**: 0 3 @@ -1965,7 +1965,7 @@ These variables are used to control the output of properties. - auto: These files are saved in folder `OUT.${suffix}/restart/`; - other: These files are saved in folder `${read_file_dir}/restart/`. - If EXX(exact exchange) is calculated (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*), the Hexx(R) files for each processor will also be saved in the above folder, which can be read in EXX calculation with *[restart_load](#restart_load)==True*. + If EXX(exact exchange) is calculated (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0* or *[rpa](#rpa)==True*), the Hexx(R) files for each processor will also be saved in the above folder, which can be read in EXX calculation with *[restart_load](#restart_load)==True*. - **Default**: False ### rpa (Under Development Feature) @@ -2927,14 +2927,14 @@ These variables are relevant when using hybrid functionals with *[basis_type](#b ### exx_opt_orb_lmax - **Type**: Integer -- **Availability**: *[dft_functional](#dft_functional)==opt_orb* +- **Availability**: *[calculation](#calculation)==gen_opt_abfs* - **Description**: The maximum l of the spherical Bessel functions, when the radial part of opt-ABFs are generated as linear combinations of spherical Bessel functions. A reasonable choice is 2. - **Default**: 0 ### exx_opt_orb_ecut - **Type**: Real -- **Availability**: *[dft_functional](#dft_functional)==opt_orb* +- **Availability**: *[calculation](#calculation)==gen_opt_abfs* - **Description**: The cut-off of plane wave expansion, when the plane wave basis is used to optimize the radial ABFs. A reasonable choice is 60. - **Default**: 0 - **Unit**: Ry @@ -2942,7 +2942,7 @@ These variables are relevant when using hybrid functionals with *[basis_type](#b ### exx_opt_orb_tolerence - **Type**: Real -- **Availability**: *[dft_functional](#dft_functional)==opt_orb* +- **Availability**: *[calculation](#calculation)==gen_opt_abfs* - **Description**: The threshold when solving for the zeros of spherical Bessel functions. A reasonable choice is 1e-12. - **Default**: 0 @@ -2975,7 +2975,7 @@ These variables are relevant when using hybrid functionals with *[basis_type](#b ### exx_symmetry_realspace - **Type**: Boolean -- **Availability**: *[symmetry](#symmetry)==1* and exx calculation (*[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*) +- **Availability**: *[symmetry](#symmetry)==1* and exx calculation (*[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0* or *[rpa](#rpa)==True*) - **Description**: - False: only rotate k-space density matrix D(k) from irreducible k-points to accelerate diagonalization - True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index cff799a8eb..84844ab720 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -58,6 +58,10 @@ #include "source_lcao/hamilt_lcao.h" #include "source_hsolver/hsolver_lcao.h" +#ifdef __EXX +#include "../source_lcao/module_ri/exx_opt_orb.h" +#endif + // test RDMFT #include "source_lcao/module_rdmft/rdmft.h" @@ -132,6 +136,17 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa two_center_bundle_, orb_); + if (PARAM.inp.calculation == "gen_opt_abfs") + { + #ifdef __EXX + Exx_Opt_Orb exx_opt_orb; + exx_opt_orb.generate_matrix(GlobalC::exx_info.info_opt_abfs, this->kv, ucell, this->orb_); + #else + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_all_runners", "calculation=gen_opt_abfs must compile __EXX"); + #endif + return; + } + // 4) initialize electronic wave function psi if (this->psi == nullptr) { diff --git a/source/source_esolver/lcao_others.cpp b/source/source_esolver/lcao_others.cpp index 850c74d393..2cc85ca253 100644 --- a/source/source_esolver/lcao_others.cpp +++ b/source/source_esolver/lcao_others.cpp @@ -74,6 +74,10 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "testing neighbour"); return; } + else if (cal_type == "gen_opt_abfs") + { + return; + } // 1. prepare HS matrices, prepare grid integral // (1) Find adjacent atoms for each atom. diff --git a/source/source_hamilt/module_xc/exx_info.h b/source/source_hamilt/module_xc/exx_info.h index 2a8b4269b2..0b427b438d 100644 --- a/source/source_hamilt/module_xc/exx_info.h +++ b/source/source_hamilt/module_xc/exx_info.h @@ -77,6 +77,15 @@ struct Exx_Info }; Exx_Info_RI info_ri; + struct Exx_Info_Opt_ABFs + { + int abfs_Lmax = 0; // tmp + double ecut_exx = 60; + double tolerence = 1E-2; + double kmesh_times = 4; + }; + Exx_Info_Opt_ABFs info_opt_abfs; + Exx_Info() : info_lip(this->info_global), info_ri(this->info_global) { } diff --git a/source/source_io/input_conv.cpp b/source/source_io/input_conv.cpp index 330c5514f9..37cc906c6d 100644 --- a/source/source_io/input_conv.cpp +++ b/source/source_io/input_conv.cpp @@ -271,7 +271,6 @@ void Input_Conv::Convert() ModuleBase::GlobalFunc::MAKE_DIR(GlobalC::restart.folder); if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0" || dft_functional_lower == "hse" - || dft_functional_lower == "opt_orb" || dft_functional_lower == "scan0") { GlobalC::restart.info_save.save_charge = true; GlobalC::restart.info_save.save_H = true; @@ -297,7 +296,6 @@ void Input_Conv::Convert() GlobalC::restart.folder = PARAM.globalv.global_readin_dir + "restart/"; if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0" || dft_functional_lower == "hse" - || dft_functional_lower == "opt_orb" || dft_functional_lower == "scan0" || dft_functional_lower == "lc_pbe" || dft_functional_lower == "lc_wpbe" @@ -422,13 +420,6 @@ void Input_Conv::Convert() } } } -#ifdef __EXX - else if (dft_functional_lower == "opt_orb") - { - GlobalC::exx_info.info_global.cal_exx = false; - Exx_Abfs::Jle::generate_matrix = true; - } -#endif else { GlobalC::exx_info.info_global.cal_exx = false; @@ -456,7 +447,6 @@ void Input_Conv::Convert() if (GlobalC::exx_info.info_global.cal_exx #ifdef __EXX - || Exx_Abfs::Jle::generate_matrix || PARAM.inp.rpa #endif ) @@ -483,11 +473,9 @@ void Input_Conv::Convert() GlobalC::exx_info.info_ri.ccp_rmesh_times = std::stod(PARAM.inp.exx_ccp_rmesh_times); GlobalC::exx_info.info_ri.exx_symmetry_realspace = PARAM.inp.exx_symmetry_realspace; -#ifdef __EXX - Exx_Abfs::Jle::Lmax = PARAM.inp.exx_opt_orb_lmax; - Exx_Abfs::Jle::Ecut_exx = PARAM.inp.exx_opt_orb_ecut; - Exx_Abfs::Jle::tolerence = PARAM.inp.exx_opt_orb_tolerence; -#endif + GlobalC::exx_info.info_opt_abfs.abfs_Lmax = PARAM.inp.exx_opt_orb_lmax; + GlobalC::exx_info.info_opt_abfs.ecut_exx = PARAM.inp.exx_opt_orb_ecut; + GlobalC::exx_info.info_opt_abfs.tolerence = PARAM.inp.exx_opt_orb_tolerence; // EXX does not support symmetry for nspin==4 if (PARAM.inp.calculation != "nscf" && PARAM.inp.symmetry == "1" && PARAM.inp.nspin == 4 && PARAM.inp.basis_type == "lcao") diff --git a/source/source_io/read_input_item_system.cpp b/source/source_io/read_input_item_system.cpp index ed8cd2f778..378d77fcd1 100644 --- a/source/source_io/read_input_item_system.cpp +++ b/source/source_io/read_input_item_system.cpp @@ -62,7 +62,7 @@ void ReadInput::item_system() } { Input_Item item("calculation"); - item.annotation = "test; scf; relax; nscf; get_wf; get_pchg"; + item.annotation = "scf; relax; md; cell-relax; nscf; get_s; get_wf; get_pchg; gen_bessel; gen_opt_abfs; test_memory; test_neighbour"; item.read_value = [](const Input_Item& item, Parameter& para) { para.input.calculation = strvalue; std::string& calculation = para.input.calculation; @@ -73,13 +73,14 @@ void ReadInput::item_system() "relax", "md", "cell-relax", - "test_memory", - "test_neighbour", "nscf", "get_s", "get_wf", "get_pchg", - "gen_bessel"}; + "gen_bessel", + "gen_opt_abfs", + "test_memory", + "test_neighbour"}; if (std::find(callist.begin(), callist.end(), calculation) == callist.end()) { const std::string warningstr = nofound_str(callist, "calculation"); diff --git a/source/source_lcao/module_ri/Exx_LRI_interface.hpp b/source/source_lcao/module_ri/Exx_LRI_interface.hpp index ff16edfac1..0c78872433 100644 --- a/source/source_lcao/module_ri/Exx_LRI_interface.hpp +++ b/source/source_lcao/module_ri/Exx_LRI_interface.hpp @@ -4,7 +4,6 @@ #include "Exx_LRI_interface.h" #include "source_lcao/module_ri/exx_abfs-jle.h" -#include "source_lcao/module_ri/exx_opt_orb.h" #include "source_lcao/hamilt_lcao.h" #include "source_lcao/module_operator_lcao/op_exx_lcao.h" #include "source_base/parallel_common.h" @@ -183,15 +182,6 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, this->cal_exx_ions(ucell,PARAM.inp.out_ri_cv); } - if (Exx_Abfs::Jle::generate_matrix) - { - //program should be stopped after this judgement - Exx_Opt_Orb exx_opt_orb; - exx_opt_orb.generate_matrix(kv, ucell,orb); - ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf"); - return; - } - // set initial parameter for mix_DMk_2D if(GlobalC::exx_info.info_global.cal_exx) { diff --git a/source/source_lcao/module_ri/Matrix_Orbs22.cpp b/source/source_lcao/module_ri/Matrix_Orbs22.cpp index 9fddea6a10..01d0330406 100644 --- a/source/source_lcao/module_ri/Matrix_Orbs22.cpp +++ b/source/source_lcao/module_ri/Matrix_Orbs22.cpp @@ -55,7 +55,6 @@ void Matrix_Orbs22::init(const int mode, Lmax = 2 * Lmax + 1; ModuleBase::timer::tick("Matrix_Orbs22", "init"); - std::cout << "Matrix_Orbs22::init()::done" << std::endl; } void Matrix_Orbs22::init_radial(const std::vector>>& orb_A1, diff --git a/source/source_lcao/module_ri/exx_abfs-construct_orbs.cpp b/source/source_lcao/module_ri/exx_abfs-construct_orbs.cpp index dff90d6de1..1556335572 100644 --- a/source/source_lcao/module_ri/exx_abfs-construct_orbs.cpp +++ b/source/source_lcao/module_ri/exx_abfs-construct_orbs.cpp @@ -45,16 +45,7 @@ std::vector>> Exx_Abfs::Construct_ } } - for (int T = 0; T < orbs.size() ; T++) - { - for (int L=orbs[T].size()-1; L >= 0 ; L--) - { - if (orbs[T][L].size()>0) - break; - else - orbs[T].resize(L); - } - } + Exx_Abfs::Construct_Orbs::filter_empty_orbs(orbs); return orbs; } @@ -548,4 +539,18 @@ std::vector Exx_Abfs::Construct_Orbs::get_Rcut(const std::vector>> &orbs) +{ + for (int T=0; T=0 ; --L) + { + if (orbs[T][L].size()>0) + { break; } + else + { orbs[T].resize(L); } + } + } +} diff --git a/source/source_lcao/module_ri/exx_abfs-construct_orbs.h b/source/source_lcao/module_ri/exx_abfs-construct_orbs.h index c8079b13dd..ae698b33bd 100644 --- a/source/source_lcao/module_ri/exx_abfs-construct_orbs.h +++ b/source/source_lcao/module_ri/exx_abfs-construct_orbs.h @@ -57,6 +57,9 @@ class Exx_Abfs::Construct_Orbs return get_Rmax(rcut); } + static void filter_empty_orbs( + std::vector>> &orbs); + private: static std::vector>>> psi_mult_psi( const std::vector>> &lcaos ); diff --git a/source/source_lcao/module_ri/exx_abfs-jle.cpp b/source/source_lcao/module_ri/exx_abfs-jle.cpp index e9de7fabe2..97dd5e9153 100644 --- a/source/source_lcao/module_ri/exx_abfs-jle.cpp +++ b/source/source_lcao/module_ri/exx_abfs-jle.cpp @@ -1,37 +1,34 @@ #include "exx_abfs-jle.h" #include "source_io/module_parameter/parameter.h" -#include "../../source_pw/module_pwdft/global.h" #include "../../source_basis/module_ao/ORB_read.h" -#include "../../source_base/global_function.h" +#include "../../source_cell/unitcell.h" #include "../../source_base/mathzone.h" #include "../../source_base/math_sphbes.h" // mohan add 2021-05-06 -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; - -void Exx_Abfs::Jle::init_jle(const double kmesh_times, - const UnitCell& ucell, - const LCAO_Orbitals& orb) +std::vector>> +Exx_Abfs::Jle::init_jle( + const Exx_Info::Exx_Info_Opt_ABFs &info, + const double kmesh_times, + const UnitCell& ucell, + const LCAO_Orbitals& orb) { - jle.resize( ucell.ntype ); + std::vector>> jle( ucell.ntype ); - for (int T = 0; T < ucell.ntype ; T++) + for(int T=0; T( sqrt( Ecut_exx ) * orb.Phi[T].getRcut() / ModuleBase::PI ); // Rydberg Unit. + = static_cast( std::sqrt( info.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), orb.Phi[T].getRcut()); + ModuleBase::Sphbes::Spherical_Bessel_Roots(ecut_number, L, info.tolerence, en.data(), orb.Phi[T].getRcut()); - for(size_t E=0; E!=ecut_number; ++E) + for(size_t E=0; E jle_r( orb.Phi[T].PhiLN(0,0).getNr() ); ModuleBase::Sphbes::Spherical_Bessel( @@ -39,7 +36,7 @@ void Exx_Abfs::Jle::init_jle(const double kmesh_times, orb.Phi[T].PhiLN(0,0).getRadial(), en[E], L, - ModuleBase::GlobalFunc::VECTOR_TO_PTR(jle_r)); + jle_r.data()); jle[T][L][E].set_orbital_info( orb.Phi[T].PhiLN(0,0).getLabel(), orb.Phi[T].PhiLN(0,0).getType(), @@ -49,13 +46,16 @@ void Exx_Abfs::Jle::init_jle(const double kmesh_times, 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), + jle_r.data(), 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); + true, + PARAM.inp.cal_force); } } } + + return jle; } diff --git a/source/source_lcao/module_ri/exx_abfs-jle.h b/source/source_lcao/module_ri/exx_abfs-jle.h index d1ea211eaf..66e9b3975d 100644 --- a/source/source_lcao/module_ri/exx_abfs-jle.h +++ b/source/source_lcao/module_ri/exx_abfs-jle.h @@ -2,28 +2,23 @@ #define EXX_ABFS_JLE_H #include "exx_abfs.h" +#include "../../source_hamilt/module_xc/exx_info.h" +#include "../../source_basis/module_ao/ORB_atomic_lm.h" #include -#include "source_cell/unitcell.h" -#include "../../source_basis/module_ao/ORB_read.h" + + class LCAO_Orbitals; + class UnitCell; class Exx_Abfs::Jle { public: - - std::vector< - std::vector< - std::vector< - Numerical_Orbital_Lm>>> jle; - - void init_jle(const double kmesh_times, - const UnitCell& ucell, - const LCAO_Orbitals& orb); - - static bool generate_matrix; - static int Lmax; - static double Ecut_exx; - static double tolerence; + static std::vector>> + init_jle( + const Exx_Info::Exx_Info_Opt_ABFs &info, + const double kmesh_times, + const UnitCell& ucell, + const LCAO_Orbitals& orb); }; #endif // EXX_ABFS_JLE_H diff --git a/source/source_lcao/module_ri/exx_opt_orb-print.cpp b/source/source_lcao/module_ri/exx_opt_orb-print.cpp index 0598eb9aa1..a851d2c579 100644 --- a/source/source_lcao/module_ri/exx_opt_orb-print.cpp +++ b/source/source_lcao/module_ri/exx_opt_orb-print.cpp @@ -3,6 +3,7 @@ #include "exx_abfs-jle.h" void Exx_Opt_Orb::print_matrix( + const Exx_Info::Exx_Info_Opt_ABFs &info, const UnitCell& ucell, const K_Vectors &kv, const std::string& file_name, @@ -60,47 +61,44 @@ void Exx_Opt_Orb::print_matrix( } // ecutwfc_jlq determine the jlq corresponding to plane wave calculation. - ofs << Exx_Abfs::Jle::Ecut_exx << " ecutwfc" << std::endl; // mohan add 2009-09-08 + ofs << info.ecut_exx << " ecutwfc" << std::endl; // mohan add 2009-09-08 // this parameter determine the total number of jlq. - ofs << Exx_Abfs::Jle::Ecut_exx << " ecutwfc_jlq" << std::endl;//mohan modify 2009-09-08 + ofs << info.ecut_exx << " ecutwfc_jlq" << std::endl;//mohan modify 2009-09-08 - if(TA==TB) { - ofs << orb_cutoff[TA] << " rcut_Jlq" << std::endl; - } else { - ofs << orb_cutoff[TA] << " " << orb_cutoff[TB] << " rcut_Jlq" << std::endl; -} + if(TA==TB) + { ofs << orb_cutoff[TA] << " rcut_Jlq" << std::endl; } + else + { 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; ofs << 0 << " smearing_sigma" << std::endl; - ofs << Exx_Abfs::Jle::tolerence << " tolerence" << std::endl; + ofs << info.tolerence << " tolerence" << std::endl; - ofs << Exx_Abfs::Jle::Lmax << " lmax" << std::endl; + ofs << info.abfs_Lmax << " lmax" << std::endl; ofs << kv.get_nkstot() << " nks" << std::endl; - assert( matrix_V.shape[0] == matrix_V.shape[1] ); - ofs << matrix_V.shape[0] << " nbands" << std::endl; + assert( matrix_V.shape[0]*matrix_V.shape[1] == matrix_V.shape[2]*matrix_V.shape[3] ); + ofs << matrix_V.shape[0]*matrix_V.shape[1] << " nbands" << std::endl; auto cal_sum_M = [&range_jles](size_t T) -> size_t { size_t sum_M = 0; - for( size_t L = 0; L!=range_jles[T].size(); ++L ) { - sum_M += range_jles[T][L].M; -} + for( size_t L = 0; L!=range_jles[T].size(); ++L ) + { sum_M += range_jles[T][L].M; } return sum_M; }; 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 ) * 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 { - ofs << ecut_numberA << " " << ecut_numberB << " ne" << std::endl; -} + const size_t ecut_numberA = static_cast( std::sqrt( info.ecut_exx ) * orb_cutoff[TA] / ModuleBase::PI ); // Rydberg Unit + const size_t ecut_numberB = static_cast( std::sqrt( info.ecut_exx ) * orb_cutoff[TB] / ModuleBase::PI ); // Rydberg Unit + if(TA==TB) + { ofs << ecut_numberA << " ne" << std::endl; } + else + { ofs << ecut_numberA << " " << ecut_numberB << " ne" << std::endl; } ofs << "" << std::endl; for( int ik=0; ik!=kv.get_nkstot(); ++ik ) @@ -119,26 +117,27 @@ void Exx_Opt_Orb::print_matrix( //--------------------- // < Psi | jY > //--------------------- - ofs<< "" << std::endl; - - for( size_t ib=0; ib!=matrix_V.shape[0]; ++ib ) + ofs<< "" << std::endl; + for( size_t iw0=0; iw0!=matrix_V.shape[0]; ++iw0 ) { - for( size_t iat=0; iat!=matrix_Q.size(); ++iat ) + for( size_t iw1=0; iw1!=matrix_V.shape[1]; ++iw1 ) { - const size_t it = (iat==0) ? TA : TB; - for( size_t il=0; il!=range_jles[it].size(); ++il ) + for( size_t iat=0; iat!=matrix_Q.size(); ++iat ) { - for( size_t im=0; im!=range_jles[it][il].M; ++im ) + const size_t it = (iat==0) ? TA : TB; + for( size_t il=0; il!=range_jles[it].size(); ++il ) { - for( size_t iq=0; iq!=range_jles[it][il].N; ++iq ) + for( size_t im=0; im!=range_jles[it][il].M; ++im ) { - ofs<" << std::endl << std::endl; }; @@ -149,7 +148,6 @@ void Exx_Opt_Orb::print_matrix( // < jY | jY > //--------------------- ofs<< "" <" << std::endl << std::endl; }; @@ -188,23 +185,26 @@ void Exx_Opt_Orb::print_matrix( // < Psi | Psi > //--------------------- ofs << "" << std::endl; - - for( size_t ib1=0; ib1!=matrix_V.shape[0]; ++ib1 ) + for( size_t iw0=0; iw0!=matrix_V.shape[0]; ++iw0 ) { - for( size_t ib2=0; ib2!=matrix_V.shape[1]; ++ib2 ) + for( size_t iw1=0; iw1!=matrix_V.shape[1]; ++iw1 ) { - ofs<" << std::endl << std::endl; }; - std::ofstream ofs(file_name+"_"+ModuleBase::GlobalFunc::TO_STRING(TA)+"_"+ModuleBase::GlobalFunc::TO_STRING(IA)+"_"+ModuleBase::GlobalFunc::TO_STRING(TB)+"_"+ModuleBase::GlobalFunc::TO_STRING(IB)); + std::ofstream ofs(file_name+"_"+std::to_string(TA)+"_"+std::to_string(IA)+"_"+std::to_string(TB)+"_"+std::to_string(IB)); print_header(ofs); print_Q(ofs); print_S(ofs); print_V(ofs); - ofs.close(); } diff --git a/source/source_lcao/module_ri/exx_opt_orb.cpp b/source/source_lcao/module_ri/exx_opt_orb.cpp index 1231552555..c22c2f8939 100644 --- a/source/source_lcao/module_ri/exx_opt_orb.cpp +++ b/source/source_lcao/module_ri/exx_opt_orb.cpp @@ -10,35 +10,31 @@ #include "source_lcao/module_ri/Matrix_Orbs21.h" #include "source_lcao/module_ri/Matrix_Orbs22.h" #include "source_lcao/module_ri/LRI_CV_Tools.h" +#include -#include "../../source_lcao/module_ri/test_code/element_basis_index-test.h" -#include "../../source_lcao/module_ri/test_code/test_function.h" -#include - -void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, const LCAO_Orbitals& orb) const +void Exx_Opt_Orb::generate_matrix( + const Exx_Info::Exx_Info_Opt_ABFs &info, + const K_Vectors &kv, + const UnitCell &ucell, + 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); - ModuleBase::TITLE("Exx_Opt_Orb::generate_matrix"); -// ofs_mpi<<"memory:\t"<>> - lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->kmesh_times ); + std::vector>> + lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, info.kmesh_times ); + Exx_Abfs::Construct_Orbs::filter_empty_orbs(lcaos); - const std::vector>> - abfs = Exx_Abfs::Construct_Orbs::abfs_same_atom(ucell,orb, lcaos, this->kmesh_times, GlobalC::exx_info.info_ri.pca_threshold ); + std::vector>> + abfs = Exx_Abfs::Construct_Orbs::abfs_same_atom(ucell,orb, lcaos, info.kmesh_times, GlobalC::exx_info.info_ri.pca_threshold ); + Exx_Abfs::Construct_Orbs::filter_empty_orbs(abfs); -// ofs_mpi<<"memory:\t"<kmesh_times, ucell , orb); + std::vector< std::vector< std::vector< Numerical_Orbital_Lm>>> + jle = Exx_Abfs::Jle::init_jle(info, info.kmesh_times, ucell , orb); + Exx_Abfs::Construct_Orbs::filter_empty_orbs(jle); -// ofs_mpi<<"memory:\t"<(abfs[T].size())-1 ); -} + GlobalC::exx_info.info_ri.abfs_Lmax = info.abfs_Lmax; + for( 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 ); } 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 ); @@ -46,28 +42,12 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co const ModuleBase::Element_Basis_Index::Range range_abfs = Exx_Abfs::Abfs_Index::construct_range( abfs ); const ModuleBase::Element_Basis_Index::IndexLNM index_abfs = ModuleBase::Element_Basis_Index::construct_index( range_abfs ); - const ModuleBase::Element_Basis_Index::Range range_jys = Exx_Abfs::Abfs_Index::construct_range( jle.jle ); + const ModuleBase::Element_Basis_Index::Range range_jys = Exx_Abfs::Abfs_Index::construct_range( jle ); const ModuleBase::Element_Basis_Index::IndexLNM index_jys = ModuleBase::Element_Basis_Index::construct_index( range_jys ); -// ofs_mpi<>> radial_R = get_radial_R(ucell); -#if TEST_EXX_RADIAL==2 - { - for(const auto & rA : radial_R) - for(const auto & rB : rA.second) - { - ofs_mpi<>> radial_R = get_radial_R(ucell); // < lcaos lcaos | lcaos lcaos > const auto ms_lcaoslcaos_lcaoslcaos = [&]() -> std::map>>>> @@ -75,7 +55,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co Matrix_Orbs22 m_lcaoslcaos_lcaoslcaos; ORB_gaunt_table MGT; int Lmax; - m_lcaoslcaos_lcaoslcaos.init( 1, ucell,orb, this->kmesh_times, orb.get_Rmax(), Lmax ); + m_lcaoslcaos_lcaoslcaos.init( 1, ucell,orb, info.kmesh_times, orb.get_Rmax(), Lmax ); MGT.init_Gaunt_CH(Lmax); MGT.init_Gaunt(Lmax); m_lcaoslcaos_lcaoslcaos.init_radial( lcaos, lcaos, lcaos, lcaos, MGT ); @@ -86,8 +66,6 @@ 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"< const auto ms_lcaoslcaos_jys = [&]() -> std::map>>>>> @@ -95,10 +73,10 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co Matrix_Orbs21 m_jyslcaos_lcaos; ORB_gaunt_table MGT; int Lmax; - m_jyslcaos_lcaos.init( 1, ucell , orb, this->kmesh_times, orb.get_Rmax(), Lmax ); + m_jyslcaos_lcaos.init( 1, ucell , orb, info.kmesh_times, orb.get_Rmax(), Lmax ); MGT.init_Gaunt_CH(Lmax); MGT.init_Gaunt(Lmax); - m_jyslcaos_lcaos.init_radial( jle.jle, lcaos, lcaos, MGT); + m_jyslcaos_lcaos.init_radial( jle, lcaos, lcaos, MGT); #if TEST_EXX_RADIAL>=1 m_jyslcaos_lcaos.init_radial_table( radial_R); #else @@ -107,18 +85,16 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co return m_jyslcaos_lcaos.cal_overlap_matrix_all(ucell,index_jys, index_lcaos, index_lcaos ); }(); -// ofs_mpi<<"memory:\t"< const auto ms_jys_jys = [&]() -> std::map>>>> { Matrix_Orbs11 m_jys_jys; ORB_gaunt_table MGT; int Lmax; - m_jys_jys.init( 2,ucell,orb, this->kmesh_times, orb.get_Rmax(), Lmax ); + m_jys_jys.init( 2,ucell,orb, info.kmesh_times, orb.get_Rmax(), Lmax ); MGT.init_Gaunt_CH(Lmax); MGT.init_Gaunt(Lmax); - m_jys_jys.init_radial( jle.jle, jle.jle, MGT ); + m_jys_jys.init_radial( jle, jle, MGT ); #if TEST_EXX_RADIAL>=1 m_jys_jys.init_radial_table(radial_R); #else @@ -127,15 +103,13 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co return m_jys_jys.cal_overlap_matrix_all(ucell,index_jys, index_jys ); }(); -// ofs_mpi<<"memory:\t"< const auto ms_abfs_abfs = [&]() -> std::map>>>> { Matrix_Orbs11 m_abfs_abfs; ORB_gaunt_table MGT; int Lmax; - m_abfs_abfs.init( 2, ucell, orb, this->kmesh_times, orb.get_Rmax(), Lmax ); + m_abfs_abfs.init( 2, ucell, orb, info.kmesh_times, orb.get_Rmax(), Lmax ); MGT.init_Gaunt_CH(Lmax); MGT.init_Gaunt(Lmax); m_abfs_abfs.init_radial( abfs, abfs, MGT ); @@ -147,15 +121,13 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co return m_abfs_abfs.cal_overlap_matrix_all(ucell,index_abfs, index_abfs ); }(); -// ofs_mpi<<"memory:\t"< const auto ms_lcaoslcaos_abfs = [&]() -> std::map>>>>> { Matrix_Orbs21 m_abfslcaos_lcaos; ORB_gaunt_table MGT; int Lmax; - m_abfslcaos_lcaos.init( 1, ucell , orb, this->kmesh_times, orb.get_Rmax(), Lmax ); + m_abfslcaos_lcaos.init( 1, ucell , orb, info.kmesh_times, orb.get_Rmax(), Lmax ); MGT.init_Gaunt_CH(Lmax); MGT.init_Gaunt(Lmax); m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos, MGT ); @@ -167,18 +139,16 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co return m_abfslcaos_lcaos.cal_overlap_matrix_all(ucell,index_abfs, index_lcaos, index_lcaos ); }(); -// ofs_mpi<<"memory:\t"< const auto ms_jys_abfs = [&]() -> std::map>>>> { Matrix_Orbs11 m_jys_abfs; ORB_gaunt_table MGT; int Lmax; - m_jys_abfs.init( 2, ucell,orb, this->kmesh_times, orb.get_Rmax(), Lmax ); + m_jys_abfs.init( 2, ucell,orb, info.kmesh_times, orb.get_Rmax(), Lmax ); MGT.init_Gaunt_CH(Lmax); MGT.init_Gaunt(Lmax); - m_jys_abfs.init_radial( jle.jle, abfs, MGT ); + m_jys_abfs.init_radial( jle, abfs, MGT ); #if TEST_EXX_RADIAL>=1 m_jys_abfs.init_radial_table(radial_R); #else @@ -187,15 +157,6 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co return m_jys_abfs.cal_overlap_matrix_all(ucell,index_jys, index_abfs ); }(); -// ofs_mpi<<"memory:\t"<>> ms_abfs_abfs_I = cal_I( ms_abfs_abfs, T,I,T,I ); // < lcaos lcaos | lcaos lcaos > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | lcaos lcaos > const RI::Tensor m_lcaoslcaos_lcaoslcaos_proj = - cal_proj( + cal_proj_22( ms_lcaoslcaos_lcaoslcaos.at(T).at(I).at(T).at(I), ms_lcaoslcaos_abfs.at(T).at(I).at(T).at(I), ms_abfs_abfs_I, ms_lcaoslcaos_abfs.at(T).at(I).at(T).at(I)); // < lcaos lcaos | jys > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector> m_lcaoslcaos_jys_proj = - {cal_proj( + {cal_proj_21( ms_lcaoslcaos_jys.at(T).at(I).at(T).at(I)[0], ms_lcaoslcaos_abfs.at(T).at(I).at(T).at(I), ms_abfs_abfs_I, {ms_jys_abfs.at(T).at(I).at(T).at(I)})}; // < jys | jys > - < jys | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector>> m_jys_jys_proj = - {{cal_proj( + {{cal_proj_11( ms_jys_jys.at(T).at(I).at(T).at(I), {ms_jys_abfs.at(T).at(I).at(T).at(I)}, ms_abfs_abfs_I, {ms_jys_abfs.at(T).at(I).at(T).at(I)})}}; - print_matrix(ucell, + print_matrix( + info, + ucell, kv, - "matrix", + PARAM.globalv.global_out_dir+"/matrix-opt-abfs", m_lcaoslcaos_jys_proj, m_jys_jys_proj, m_lcaoslcaos_lcaoslcaos_proj, @@ -244,9 +207,11 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co } else { - print_matrix(ucell, + print_matrix( + info, + ucell, kv, - "matrix", + PARAM.globalv.global_out_dir+"/matrix-opt-abfs", ms_lcaoslcaos_jys.at(T).at(I).at(T).at(I), {{ms_jys_jys.at(T).at(I).at(T).at(I)}}, ms_lcaoslcaos_lcaoslcaos.at(T).at(I).at(T).at(I), @@ -263,48 +228,50 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co const std::vector>> ms_abfs_abfs_I = cal_I( ms_abfs_abfs, TA,IA,TB,IB ); // < lcaos lcaos | lcaos lcaos > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | lcaos lcaos > const RI::Tensor m_lcaoslcaos_lcaoslcaos_proj = - cal_proj( + cal_proj_22( ms_lcaoslcaos_lcaoslcaos.at(TA).at(IA).at(TB).at(IB), ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB), ms_abfs_abfs_I, ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB)); // < lcaos lcaos | jys > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector> m_lcaoslcaos_jys_proj = - {cal_proj( + {cal_proj_21( ms_lcaoslcaos_jys.at(TA).at(IA).at(TB).at(IB)[0], ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB), ms_abfs_abfs_I, { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }), - cal_proj( + cal_proj_21( ms_lcaoslcaos_jys.at(TA).at(IA).at(TB).at(IB)[1], ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB), ms_abfs_abfs_I, { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) })}; // < jys | jys > - < jys | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector>> m_jys_jys_proj = - {{cal_proj( + {{cal_proj_11( ms_jys_jys.at(TA).at(IA).at(TA).at(IA), { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }), - cal_proj( + cal_proj_11( ms_jys_jys.at(TA).at(IA).at(TB).at(IB), { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }) }, - {cal_proj( + {cal_proj_11( ms_jys_jys.at(TB).at(IB).at(TA).at(IA), { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }), - cal_proj( + cal_proj_11( ms_jys_jys.at(TB).at(IB).at(TB).at(IB), { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }) }}; - print_matrix(ucell, + print_matrix( + info, + ucell, kv, - "matrix", + PARAM.globalv.global_out_dir+"/matrix-opt-abfs", m_lcaoslcaos_jys_proj, m_jys_jys_proj, m_lcaoslcaos_lcaoslcaos_proj, @@ -314,9 +281,11 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co } else { - print_matrix(ucell, + print_matrix( + info, + ucell, kv, - "matrix", + PARAM.globalv.global_out_dir+"/matrix-opt-abfs", ms_lcaoslcaos_jys.at(TA).at(IA).at(TB).at(IB), {{ms_jys_jys.at(TA).at(IA).at(TA).at(IA), ms_jys_jys.at(TA).at(IA).at(TB).at(IB)}, {ms_jys_jys.at(TB).at(IB).at(TA).at(IA), ms_jys_jys.at(TB).at(IB).at(TB).at(IB)}}, @@ -333,27 +302,82 @@ 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_22( + 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"); - -//auto print_nrc = [](const matrix & m){ std::cout<<"\t"< m_proj = m_big; -//print_nrc(m_proj); + ModuleBase::TITLE("Exx_Opt_Orb::cal_proj_22"); + RI::Tensor m_proj = m_big.copy(); + for( size_t il=0; il!=m_left.size(); ++il ) + { + for( size_t ir=0; ir!=m_right.size(); ++ir ) + { + // m_proj = m_proj - m_left[il] * m_middle[il][ir] * m_right[ir].T; + const RI::Tensor m_lm = RI::Tensor_Multiply::x0x1y1_x0x1a_ay1(m_left[il], m_middle[il][ir]); + const RI::Tensor m_lmr = RI::Tensor_Multiply::x0x1y0y1_x0x1a_y0y1a(m_lm, m_right[ir]); + m_proj -= m_lmr; + } + } + return m_proj; +} +RI::Tensor Exx_Opt_Orb::cal_proj_21( + 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_21"); + RI::Tensor m_proj = m_big.copy(); + for( size_t il=0; il!=m_left.size(); ++il ) + { + for( size_t ir=0; ir!=m_right.size(); ++ir ) + { + // m_proj = m_proj - m_left[il] * m_middle[il][ir] * m_right[ir].T; + const RI::Tensor m_lm = RI::Tensor_Multiply::x0x1y1_x0x1a_ay1(m_left[il], m_middle[il][ir]); + const RI::Tensor m_lmr = RI::Tensor_Multiply::x0x1y0_x0x1a_y0a(m_lm, m_right[ir]); + m_proj -= m_lmr; + } + } + return m_proj; +} +RI::Tensor Exx_Opt_Orb::cal_proj_12( + 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_12"); + RI::Tensor m_proj = m_big.copy(); + for( size_t il=0; il!=m_left.size(); ++il ) + { + for( size_t ir=0; ir!=m_right.size(); ++ir ) + { + // m_proj = m_proj - m_left[il] * m_middle[il][ir] * m_right[ir].T; + const RI::Tensor m_lm = RI::Tensor_Multiply::x0y1_x0a_ay1(m_left[il], m_middle[il][ir]); + const RI::Tensor m_lmr = RI::Tensor_Multiply::x0y0y1_x0a_y0y1a(m_lm, m_right[ir]); + m_proj -= m_lmr; + } + } + return m_proj; +} +RI::Tensor Exx_Opt_Orb::cal_proj_11( + 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_11"); + RI::Tensor m_proj = m_big.copy(); for( size_t il=0; il!=m_left.size(); ++il ) { for( size_t ir=0; ir!=m_right.size(); ++ir ) { -//std::cout< m_lm = RI::Tensor_Multiply::x0y1_x0a_ay1(m_left[il], m_middle[il][ir]); + const RI::Tensor m_lmr = RI::Tensor_Multiply::x0y0_x0a_y0a(m_lm, m_right[ir]); + m_proj -= m_lmr; } } return m_proj; @@ -367,18 +391,16 @@ std::vector>> Exx_Opt_Orb::cal_I( if( TA==TB && IA==IB ) { - std::vector>> m_I - {{ RI::Tensor(ms.at(TA).at(IA).at(TA).at(IA).shape) }}; - return LRI_CV_Tools::cal_I(m_I); + return {{LRI_CV_Tools::cal_I(RI::Tensor(ms.at(TA).at(IA).at(TA).at(IA)))}}; } else { - std::vector>> m_I - {{ RI::Tensor(ms.at(TA).at(IA).at(TA).at(IA).shape), - RI::Tensor(ms.at(TA).at(IA).at(TB).at(IB).shape) }, - { RI::Tensor(ms.at(TB).at(IB).at(TA).at(IA).shape), - RI::Tensor(ms.at(TB).at(IB).at(TB).at(IB).shape) }}; - return LRI_CV_Tools::cal_I(m_I); + std::vector>> m_in + {{ ms.at(TA).at(IA).at(TA).at(IA), + ms.at(TA).at(IA).at(TB).at(IB) }, + { ms.at(TB).at(IB).at(TA).at(IA), + ms.at(TB).at(IB).at(TB).at(IB) }}; + return LRI_CV_Tools::cal_I(m_in); } } diff --git a/source/source_lcao/module_ri/exx_opt_orb.h b/source/source_lcao/module_ri/exx_opt_orb.h index 5bebe91771..2520a584ba 100644 --- a/source/source_lcao/module_ri/exx_opt_orb.h +++ b/source/source_lcao/module_ri/exx_opt_orb.h @@ -1,29 +1,50 @@ #ifndef EXX_OPT_ORB_H #define EXX_OPT_ORB_H +#include "../../source_hamilt/module_xc/exx_info.h" #include "../../source_base/matrix.h" #include "../../source_base/element_basis_index.h" #include "source_cell/klist.h" #include "source_basis/module_ao/ORB_read.h" +#include #include #include #include -#include class Exx_Opt_Orb { public: - void generate_matrix(const K_Vectors &kv, const UnitCell& ucell, const LCAO_Orbitals& orb) const; + void generate_matrix( + const Exx_Info::Exx_Info_Opt_ABFs &info, + const K_Vectors &kv, + const UnitCell &ucell, + const LCAO_Orbitals &orb) const; private: 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_22( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ) const; + RI::Tensor cal_proj_21( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ) const; + RI::Tensor cal_proj_12( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ) const; + RI::Tensor cal_proj_11( + 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 Exx_Info::Exx_Info_Opt_ABFs &info, const UnitCell& ucell, const K_Vectors &kv, const std::string& file_name, @@ -35,7 +56,5 @@ class Exx_Opt_Orb 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 diff --git a/source/source_main/driver_run.cpp b/source/source_main/driver_run.cpp index 1666498bb3..d628749e24 100644 --- a/source/source_main/driver_run.cpp +++ b/source/source_main/driver_run.cpp @@ -75,7 +75,7 @@ void Driver::driver_run() { p_esolver->runner(ucell, 0); } - else if (cal == "get_pchg" || cal == "get_wf" || cal == "gen_bessel" || + else if (cal == "get_pchg" || cal == "get_wf" || cal == "gen_bessel" || cal == "gen_opt_abfs" || cal == "test_memory" || cal == "test_neighbour") { //! supported "other" functions: