Refactor LR_IO function and add test case 58_KP_LR_BSE

.gitignore

@@ -8,6 +8,7 @@ log.txt
 result.out
 test.sum
 *.dat
+!tests/08_EXX/58_KP_LR_BSE/librpa.d/KS_eigenvector_*.dat
 .DS_Store
 .cache
 .vscode

docs/advanced/input_files/input-main.md

@@ -2370,7 +2370,7 @@
 - **Type**: String
 - **Availability**: *[`basis_type`](#basis_type)==lcao*
 - **Description**: The directory to save files for LibRPA.
-- **Default**: "./OUT.librpa/"
+- **Default**: "OUT.librpa"
 
 ### out_pchg
 

docs/parameters.yaml

@@ -3617,7 +3617,7 @@ parameters:
     type: String
     description: |
       The directory to save files for LibRPA.
-    default_value: "\"./OUT.librpa/\""
+    default_value: "\"OUT.librpa\""
     unit: ""
     availability: basis_type==lcao
   - name: out_pchg

source/source_cell/atom_spec.h

@@ -49,7 +49,7 @@ class Atom
     int stapos_wf = 0;                           ///< start position of wave functions
 
     std::string label;                            ///< atomic symbol
-    std::vector<ModuleBase::Vector3<double>> tau; ///< Cartesian coordinates of each atom in this type
+    std::vector<ModuleBase::Vector3<double>> tau; ///< Cartesian coordinates of each atom in this type, in unit of lat0
     std::vector<ModuleBase::Vector3<double>> dis; ///< direct displacements of each atom in this type in current step (liuyu modify 2023-03-22)
     std::vector<ModuleBase::Vector3<double>> taud; ///< Direct coordinates of each atom in this type
     std::vector<ModuleBase::Vector3<int>> boundary_shift; ///< record for periodic boundary adjustment

source/source_cell/klist.h

@@ -14,9 +14,9 @@
 class K_Vectors
 {
 public:
-    std::vector<ModuleBase::Vector3<double>> kvec_c; ///< Cartesian coordinates of k points
+    std::vector<ModuleBase::Vector3<double>> kvec_c; ///< Cartesian coordinates of k points, in unit of 2*pi/lat0
     std::vector<ModuleBase::Vector3<double>> kvec_d; ///< Direct coordinates of k points
-    std::vector<ModuleBase::Vector3<double>> kvec_c_full; ///< Cartesian coordinates of full k mesh match with nkstot_full
+    std::vector<ModuleBase::Vector3<double>> kvec_c_full; ///< Cartesian coordinates of full k mesh match with nkstot_full, in unit of 2*pi/lat0
 
     std::vector<double> wk; ///< wk, weight of k points
 
@@ -328,4 +328,4 @@ class K_Vectors
     friend void KVectorUtils::kvec_mpi_k(K_Vectors& kvec);
 #endif
 };
-#endif // KVECT_H
+#endif // KVECT_H

source/source_cell/unitcell_data.h

@@ -24,7 +24,7 @@ struct Lattice
     double omega = 0.0;                ///< the volume of the unit cell
     std::vector<int> lat_axis_free{0, 0, 0};      ///< whether each lattice axis (a,b,c) is allowed to relax (0=fixed, 1=free)
 
-    ModuleBase::Matrix3 latvec = ModuleBase::Matrix3();           ///< Unitcell lattice vectors
+    ModuleBase::Matrix3 latvec = ModuleBase::Matrix3();           ///< Unitcell lattice vectors, in unit of lat0
     ModuleBase::Vector3<double> a1, a2, a3;                       ///< Same as latvec, just at another form
     ModuleBase::Vector3<double> latcenter;                        ///< (a1+a2+a3)/2 the center of vector
     ModuleBase::Matrix3 latvec_supercell = ModuleBase::Matrix3(); ///< Supercell lattice vectors

source/source_esolver/esolver_lr_lcao_bse.cpp

@@ -8,6 +8,7 @@
 #include "source_lcao/module_lr/lr_spectrum.h"
 #include "source_lcao/module_lr/utils/exciton_plotter.h"
 #include "source_lcao/module_lr/utils/lr_io.h"
+#include "source_lcao/module_lr/utils/lr_util.hpp"
 
 namespace ModuleESolver
 {
@@ -28,11 +29,8 @@ void ESolver_BSE<T, TR>::before_all_runners(BaseCell& basecell, const Input_para
     ModuleESolver::ESolver_FP::before_all_runners(basecell, inp);
     this->pelec = new elecstate::ElecStateLCAO<T>();
 
-    this->kRlist = LR_IO::RI_kRlist(*this->ucell_,
-                                    &this->kv,
-                                    this->rpa_dir,
-                                    inp.bse_use_fine_kgrid,
-                                    this->out_dir);
+    this->kRlist = LR_IO::RI_kRlist(*this->ucell_, &this->kv, this->nspin,
+                                    this->rpa_dir, this->out_dir, inp.bse_use_fine_kgrid);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "Set K-POINTS and R-list for RI");
     ModuleIO::print_parameters(ucell, this->kv, inp);
 
@@ -609,6 +607,11 @@ void ESolver_BSE<T, TR>::lri_init()
             Vs_in = LR_IO::read_coulomb_mat_k<T, T>(this->rpa_dir, Cs_in, this->kRlist);
         }
         Ws_in = LR_IO::read_Ws<T, T>(Vs_in, this->kRlist.Rlist);
+        // if (GlobalV::MY_RANK == 0)
+        // {
+        //     LR_IO::write_lri_R_max_norm(Vs_in, *this->ucell_, this->out_dir + "V_R_max_norm.dat");
+        //     LR_IO::write_lri_R_max_norm(Ws_in, *this->ucell_, this->out_dir + "W_R_max_norm.dat");
+        // }
     // }
 #ifdef __MPI
     MPI_Barrier(MPI_COMM_WORLD);

source/source_esolver/esolver_lr_lcao_tddft.cpp

@@ -587,6 +587,7 @@ void ModuleESolver::ESolver_LR<T, TR>::runner(BaseCell& basecell, const int iste
                                 this->paraC_,
                                 this->paraMat_,
                                 spin_types[is],
+                                this->in_dir,
                                 this->out_dir,
                                 this->inp_->ri_hartree_benchmark,
                                 (this->inp_->ri_hartree_benchmark == "aims" ? this->inp_->aims_nbasis : std::vector<int>({})));

source/source_io/module_parameter/input_parameter.h

@@ -465,7 +465,7 @@ struct Input_para
     bool restart_save = false;               ///< restart //Peize Lin add 2020-04-04
     bool rpa = false;                        ///< rpa calculation
     bool rpa_out_vel = false;                ///< whether to output velocity matrix for librpa
-    std::string rpa_outdir = "./OUT.librpa/";///< output directory for librpa
+    std::string rpa_outdir = "OUT.librpa";   ///< output directory for librpa
     std::vector<int> out_pchg = {};          ///< specify the bands to be calculated for partial charge
     std::vector<int> out_wfc_norm = {};      ///< specify the bands to be calculated for norm of wfc
     std::vector<int> out_wfc_re_im = {};     ///< specify the bands to be calculated for real and imaginary parts of wfc

source/source_io/module_parameter/read_inp_out.cpp

@@ -1481,7 +1481,7 @@ If EXX(exact exchange) is calculated (i.e. dft_fuctional==hse/hf/pbe0/scan0 or r
         item.category = "Output information";
         item.type = "String";
         item.description = "The directory to save files for LibRPA.";
-        item.default_value = "\"./OUT.librpa/\"";
+        item.default_value = "\"OUT.librpa\"";
         item.unit = "";
         item.set_availability("basis_type==lcao");
         read_sync_string(input.rpa_outdir);

source/source_io/test/read_input_ptest.cpp

@@ -362,7 +362,7 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_FALSE(param.inp.dft_plus_dmft);
     EXPECT_FALSE(param.inp.rpa);
     EXPECT_FALSE(param.inp.rpa_out_vel);
-    EXPECT_EQ(param.inp.rpa_outdir, "./OUT.librpa/");
+    EXPECT_EQ(param.inp.rpa_outdir, "OUT.librpa/");
     EXPECT_EQ(param.inp.imp_sol, 0);
     EXPECT_DOUBLE_EQ(param.inp.eb_k, 80.0);
     EXPECT_DOUBLE_EQ(param.inp.tau, 1.0798 * 1e-5);

source/source_lcao/module_bse/hamilt_bse.cpp

@@ -142,13 +142,13 @@ HamiltBSE<T>::HamiltBSE(const int& nspin,
 
 template <typename T>
 void HamiltBSE<T>::cal_V_for_A(){
-    ModuleBase::TITLE("HamiltBSE", "cal_V_for_A");
-    ModuleBase::timer::start("HamiltBSE", "cal_V_for_A");
-    std::cout<<"in cal_V_for_A"<<std::endl;
     if (!this->VA_local.empty()) {
         std::cout<< "V for A has been calculated, skip." <<std::endl;
         return;
     }
+    ModuleBase::TITLE("HamiltBSE", "cal_V_for_A");
+    ModuleBase::timer::start("HamiltBSE", "cal_V_for_A");
+    std::cout<<"in cal_V_for_A"<<std::endl;
     BSE_Util::print_mem_estimate("V matrix of A", this->pA.get_local_size(), sizeof(T));
     this->VA_local.resize(this->pA.get_local_size(), 0.0);
     if (this->ri_hartree_benchmark == "aims" || this->ri_hartree_benchmark == "abacus") {
@@ -171,13 +171,13 @@ void HamiltBSE<T>::cal_V_for_A(){
 
 template <typename T>
 void HamiltBSE<T>::cal_V_for_B(){
-    ModuleBase::TITLE("HamiltBSE", "cal_V_for_B");
-    ModuleBase::timer::start("HamiltBSE", "cal_V_for_B");
-    std::cout<<"in cal_V_for_B"<<std::endl;
     if (!this->VB_local.empty()) {
         std::cout<< "V for B has been calculated, skip." <<std::endl;
         return;
     }
+    ModuleBase::TITLE("HamiltBSE", "cal_V_for_B");
+    ModuleBase::timer::start("HamiltBSE", "cal_V_for_B");
+    std::cout<<"in cal_V_for_B"<<std::endl;
     BSE_Util::print_mem_estimate("V matrix of B", this->pA.get_local_size(), sizeof(T));
     this->VB_local.resize(this->pA.get_local_size(), 0.0);
     if (this->ri_hartree_benchmark == "aims" || this->ri_hartree_benchmark == "abacus") {
@@ -200,13 +200,13 @@ void HamiltBSE<T>::cal_V_for_B(){
 
 template <typename T>
 void HamiltBSE<T>::cal_W_for_A(){
-    ModuleBase::TITLE("HamiltBSE", "cal_W_for_A");
-    ModuleBase::timer::start("HamiltBSE", "cal_W_for_A");
-    std::cout<<"in cal_W_for_A"<<std::endl;
     if (!this->WA_local.empty()) {
         std::cout<< "W for A has been calculated, skip." <<std::endl;
         return;
     }
+    ModuleBase::TITLE("HamiltBSE", "cal_W_for_A");
+    ModuleBase::timer::start("HamiltBSE", "cal_W_for_A");
+    std::cout<<"in cal_W_for_A"<<std::endl;
     BSE_Util::print_mem_estimate("W matrix of A", this->pA.get_local_size(), sizeof(T));
     this->WA_local.resize(this->pA.get_local_size(), 0.0);
     this->mo_lri.cal_W_for_A(this->WA_local, this->pA);    
@@ -219,13 +219,13 @@ void HamiltBSE<T>::cal_W_for_A(){
 
 template <typename T>
 void HamiltBSE<T>::cal_W_for_B(){
-    ModuleBase::TITLE("HamiltBSE", "cal_W_for_B");
-    ModuleBase::timer::start("HamiltBSE", "cal_W_for_B");
-    std::cout<<"in cal_W_for_B"<<std::endl;
     if (!this->WB_local.empty()) {
         std::cout<< "W for B has been calculated, skip." <<std::endl;
         return;
     }
+    ModuleBase::TITLE("HamiltBSE", "cal_W_for_B");
+    ModuleBase::timer::start("HamiltBSE", "cal_W_for_B");
+    std::cout<<"in cal_W_for_B"<<std::endl;
     BSE_Util::print_mem_estimate("W matrix of B", this->pA.get_local_size(), sizeof(T));
     this->WB_local.resize(this->pA.get_local_size(), 0.0);
     this->mo_lri.cal_W_for_B(this->WB_local, this->pA);
@@ -336,7 +336,7 @@ void HamiltBSE<T>::init_bse_matrix(const bool is_full, const int & st_index){
     }
     if (this->out_bse_ab)
     {
-        this->write_AB_matrix("A_matrix_"+std::to_string(this->my_rank)+".dat", 6, this->BSE_A_local.data(), this->ndim, this->ndim);
+        this->write_AB_matrix(this->out_dir+"A_matrix_"+std::to_string(this->my_rank)+".dat", 6, this->BSE_A_local.data(), this->ndim, this->ndim);
     }
 
     if (is_full)
@@ -351,7 +351,7 @@ void HamiltBSE<T>::init_bse_matrix(const bool is_full, const int & st_index){
         }
         if (this->out_bse_ab)
         {
-            this->write_AB_matrix("B_matrix_"+std::to_string(this->my_rank)+".dat", 6, this->BSE_B_local.data(), this->ndim, this->ndim);
+            this->write_AB_matrix(this->out_dir+"B_matrix_"+std::to_string(this->my_rank)+".dat", 6, this->BSE_B_local.data(), this->ndim, this->ndim);
         }
     }
 

source/source_lcao/module_lr/CMakeLists.txt

@@ -8,6 +8,7 @@ if(ENABLE_LCAO)
     utils/lr_util.cpp
     utils/lr_util_hcontainer.cpp
     utils/lr_io.cpp
+    utils/lr_io_krlist.cpp
     utils/exciton_plotter.cpp
     ao_to_mo_transformer/ao_to_mo_parallel.cpp
     ao_to_mo_transformer/ao_to_mo_serial.cpp

source/source_lcao/module_lr/hamilt_casida.h

@@ -38,6 +38,7 @@ namespace LR
                const Parallel_2D& pc_in,
                const Parallel_Orbitals& pmat_in,
                const std::string& spin_type,
+               const std::string& in_dir,
                const std::string& out_dir,
                const std::string& ri_hartree_benchmark = "none",
                const std::vector<int>& aims_nbasis = {})
@@ -56,11 +57,10 @@ namespace LR
 #ifdef __EXX
             using TAC = std::pair<int, std::array<int, 3>>;
             using TLRI = std::map<int, std::map<TAC, RI::Tensor<T>>>;
-            const std::string& dir = PARAM.globalv.global_readin_dir;
             TLRI Cs_read; 
             TLRI Vs_read; 
 #ifdef __DEBUG
-            // TLRI Vs_compare = LRI_CV_Tools::read_Vs_abf<T>(dir + "Vs");
+            // TLRI Vs_compare = LRI_CV_Tools::read_Vs_abf<T>(in_dir + "Vs");
             // LRI_CV_Tools::write_Vs_abf(Vs_read, "Vs_read_from_coulomb");
             // LRI_CV_Tools::write_Cs_ao(Cs_read, "Cs_ao_read"); // ensure Cs_ao is read correctly
             // assert(RI_Benchmark::compare_Vs(Vs_read, Vs_compare));
@@ -75,27 +75,23 @@ namespace LR
                     {
                         LR_IO::RI_kRlist kRlist (ucell_in,
                                                 const_cast<K_Vectors*>(&kv_in),
-                                                dir,
-                                                use_fine_kgrid,
-                                                out_dir);
+                                                nspin, in_dir, out_dir, use_fine_kgrid);
                         // though C and V are real, here still use <T> to multiply with psi
-                        Cs_read = LRI_CV_Tools::read_Cs_ao_all<T>(dir);
-                        Vs_read = LR_IO::read_coulomb_mat_general_k<T,T>(dir, Cs_read, kRlist);
+                        Cs_read = LRI_CV_Tools::read_Cs_ao_all<T>(in_dir);
+                        Vs_read = LR_IO::read_coulomb_mat_general_k<T,T>(in_dir, Cs_read, kRlist);
                     }
                     else if (ri_hartree_benchmark == "abacus")
                     {
-                        Cs_read = LRI_CV_Tools::read_Cs_ao<T>(dir + "Cs");
-                        Vs_read = LRI_CV_Tools::read_Vs_abf<T>(dir + "Vs");
+                        Cs_read = LRI_CV_Tools::read_Cs_ao<T>(in_dir + "Cs");
+                        Vs_read = LRI_CV_Tools::read_Vs_abf<T>(in_dir + "Vs");
                     }
                     else if (ri_hartree_benchmark == "abacus-librpa")
                     {
                         LR_IO::RI_kRlist kRlist (ucell_in,
                                                 const_cast<K_Vectors*>(&kv_in),
-                                                dir,
-                                                use_fine_kgrid,
-                                                out_dir);
-                        Cs_read = LRI_CV_Tools::read_Cs_ao_all<T>(dir);
-                        Vs_read = LR_IO::read_coulomb_mat_k<T,T>(dir, Cs_read, kRlist);
+                                                nspin, in_dir, out_dir, use_fine_kgrid);
+                        Cs_read = LRI_CV_Tools::read_Cs_ao_all<T>(in_dir);
+                        Vs_read = LR_IO::read_coulomb_mat_k<T,T>(in_dir, Cs_read, kRlist);
                     }
                     if (!std::set<std::string>({ "rpa", "hf"}).count(xc_kernel)) {
                         throw std::runtime_error("ri_hartree_benchmark is only supported for xc_kernel = rpa, hf");