diff --git a/driver/driver_utils.cpp b/driver/driver_utils.cpp index 945c935a..46fbaea8 100644 --- a/driver/driver_utils.cpp +++ b/driver/driver_utils.cpp @@ -48,30 +48,27 @@ std::vector interpolate_dielec_func(int option, const std::vector> kfrac_band = - read_band_kpath_info(driver_params.input_dir + "pyatb_librpa_df/k_path_info", - n_basis, n_states, n_spin, flag); - if (Params::use_soc) - { - assert(n_basis % 2 == 0 && "Error: nbasis is not even when SOC!"); - n_basis = n_basis / 2; - } - df_headwing.set(pyatb_meanfield, kfrac_band, frequencies_target, n_basis, n_states, + read_velocity(file_abacus, meanfield); + n_basis = meanfield.get_n_aos(); + n_states = meanfield.get_n_bands(); + n_spin = meanfield.get_n_spins(); + // int flag; + // std::vector> kfrac_band = + // read_band_kpath_info(driver_params.input_dir + "pyatb_librpa_df/k_path_info", + // n_basis, n_states, n_spin, flag); + df_headwing.set(meanfield, kfrac_list, frequencies_target, n_basis, n_states, n_spin); } else if (infile_aims.is_open()) @@ -104,9 +101,8 @@ std::vector interpolate_dielec_func(int option, const std::vector interpolate_dielec_func(int option, const std::vector> kfrac_band = - read_band_kpath_info(driver_params.input_dir + "pyatb_librpa_df/k_path_info", - n_basis, n_states, n_spin, flag); - df_headwing.set(meanfield, kfrac_band, frequencies_target, n_basis, n_states, + n_basis = meanfield.get_n_aos(); + n_states = meanfield.get_n_bands(); + n_spin = meanfield.get_n_spins(); + // int flag; + // std::vector> kfrac_band = + // read_band_kpath_info(driver_params.input_dir + "pyatb_librpa_df/k_path_info", + // n_basis, n_states, n_spin, flag); + df_headwing.set(meanfield, kfrac_list, frequencies_target, n_basis, n_states, n_spin); } else if (infile_aims.is_open()) diff --git a/driver/inputfile.cpp b/driver/inputfile.cpp index 00d7e38d..d308edd3 100644 --- a/driver/inputfile.cpp +++ b/driver/inputfile.cpp @@ -189,6 +189,10 @@ void parse_inputfile_to_params(const std::string &fn) parser.parse_bool("replace_w_head", Params::replace_w_head, true, flag); parser.parse_int("option_dielect_func", Params::option_dielect_func, 2, flag); + parser.parse_bool("band_continue", Params::band_continue, false, flag); + + parser.parse_int("output_Wc_Rf_mat", Params::output_Wc_Rf_mat, false, flag); + parser.parse_bool("output_energy_qp", Params::output_energy_qp, false, flag); parser.parse_bool("output_gw_sigc_mat", Params::output_gw_sigc_mat, false, flag); parser.parse_bool("output_gw_sigc_mat_rt", Params::output_gw_sigc_mat_rt, false, flag); parser.parse_bool("output_gw_sigc_mat_rf", Params::output_gw_sigc_mat_rf, false, flag); diff --git a/driver/task_gw.cpp b/driver/task_gw.cpp index 50cadffd..8ac3c7d3 100644 --- a/driver/task_gw.cpp +++ b/driver/task_gw.cpp @@ -408,6 +408,37 @@ void task_g0w0(std::map, ComplexMatrix> &sinvS) printf("CBM: k-point %4d: (%.5f, %.5f, %.5f) \n", ik_cond + 1, k_cond.x, k_cond.y, k_cond.z); lib_printf("Bandgap(eV): %12.7f \n", bandgap); + + if (Params::output_energy_qp) { + std::ofstream ofs("energy_qp"); + ofs << " state occ_num e_gs(Ha) e_qp(Ha)"<, ComplexMatrix> &sinvS) qlist.push_back(q_weight.first); } - // Prepare time-frequency grids - auto tfg = - LIBRPA::utils::generate_timefreq_grids(Params::nfreq, Params::tfgrids_type, meanfield); - - Chi0 chi0(meanfield, klist, tfg); - chi0.gf_R_threshold = Params::gf_R_threshold; - Profiler::start("read_vq_cut", "Load truncated Coulomb"); if (LIBRPA::parallel_routing == LIBRPA::ParallelRouting::R_TAU) { @@ -64,61 +57,21 @@ void task_g0w0_band(std::map, ComplexMatrix> &sinvS) } Profiler::cease("read_vq_cut"); - std::vector epsmac_LF_imagfreq_re; - if (Params::replace_w_head) - { - std::vector omegas_dielect; - std::vector dielect_func; - if (Params::option_dielect_func != 3 && Params::option_dielect_func != 4) - read_dielec_func(driver_params.input_dir + "dielecfunc_out", omegas_dielect, - dielect_func); - - epsmac_LF_imagfreq_re = interpolate_dielec_func(Params::option_dielect_func, omegas_dielect, - dielect_func, chi0.tfg.get_freq_nodes()); - } + Profiler::start("read_vxc", "Load DFT xc potential"); + std::vector vxc; + int flag_read_vxc = read_vxc(driver_params.input_dir + "vxc_out", vxc); + Profiler::stop("read_vxc"); - chi0.set_input_dir(driver_params.input_dir); - Profiler::start("chi0_build", "Build response function chi0"); - if (Params::use_shrink_chi) - { - chi0.build(Cs_data, Rlist, period, local_atpair, qlist, sinvS); - } - else + if (flag_read_vxc != 0) { - chi0.build(Cs_shrinked_data, Rlist, period, local_atpair, qlist, sinvS); - } - - Profiler::stop("chi0_build"); - - std::flush(ofs_myid); - mpi_comm_global_h.barrier(); - - if (Params::debug) - { // debug, check chi0 - char fn[80]; - for (const auto &chi0q : chi0.get_chi0_q()) + if (mpi_comm_global_h.myid == 0) { - const int ifreq = chi0.tfg.get_freq_index(chi0q.first); - for (const auto &q_IJchi0 : chi0q.second) - { - const int iq = std::distance(klist.begin(), - std::find(klist.begin(), klist.end(), q_IJchi0.first)); - for (const auto &I_Jchi0 : q_IJchi0.second) - { - const auto &I = I_Jchi0.first; - for (const auto &J_chi0 : I_Jchi0.second) - { - const auto &J = J_chi0.first; - sprintf(fn, "chi0fq_ifreq_%d_iq_%d_I_%zu_J_%zu_id_%d.mtx", ifreq, iq, I, J, - mpi_comm_global_h.myid); - print_complex_matrix_mm(J_chi0.second, Params::output_dir + "/" + fn, - 1e-15); - } - } - } + lib_printf("Error in reading Vxc on kgrid, task failed!\n"); } + Profiler::stop("g0w0_band"); + return; } - + // ================== EXX part: Compute Hexx(R) ================== Profiler::start("g0w0_exx", "Build exchange self-energy"); auto exx = LIBRPA::Exx(meanfield, kfrac_list, period); { @@ -156,97 +109,151 @@ void task_g0w0_band(std::map, ComplexMatrix> &sinvS) Profiler::stop("g0w0_exx"); std::flush(ofs_myid); - Profiler::start("g0w0_wc", "Build screened interaction"); - vector> epsmac_LF_imagfreq(epsmac_LF_imagfreq_re.cbegin(), - epsmac_LF_imagfreq_re.cend()); - map, matrix_m>>>::pair_t_old> - Wc_freq_q; - if (Params::use_scalapack_gw_wc) + // ================== GW part: Compute Sigma_c(R, iw) ================== + // Prepare time-frequency grids + auto tfg = + LIBRPA::utils::generate_timefreq_grids(Params::nfreq, Params::tfgrids_type, meanfield); + + LIBRPA::G0W0 s_g0w0(meanfield, kfrac_list, tfg, period); + + if (Params::band_continue) { - if (Params::use_fullcoul_wc) + s_g0w0.read_sigc(driver_params.input_dir + "librpa.d/", Rlist); + } + else + { + Chi0 chi0(meanfield, klist, tfg); + chi0.gf_R_threshold = Params::gf_R_threshold; + + std::vector epsmac_LF_imagfreq_re; + if (Params::replace_w_head) + { + std::vector omegas_dielect; + std::vector dielect_func; + if (Params::option_dielect_func != 3 && Params::option_dielect_func != 4) + read_dielec_func(driver_params.input_dir + "dielecfunc_out", omegas_dielect, + dielect_func); + + epsmac_LF_imagfreq_re = interpolate_dielec_func(Params::option_dielect_func, omegas_dielect, + dielect_func, chi0.tfg.get_freq_nodes()); + } + + chi0.set_input_dir(driver_params.input_dir); + Profiler::start("chi0_build", "Build response function chi0"); + if (Params::use_shrink_chi) { - Wc_freq_q = compute_Wc_freq_q_blacs(chi0, Vq, Vq, epsmac_LF_imagfreq); + chi0.build(Cs_data, Rlist, period, local_atpair, qlist, sinvS); } else { - Wc_freq_q = compute_Wc_freq_q_blacs(chi0, Vq, Vq_cut, epsmac_LF_imagfreq); + chi0.build(Cs_shrinked_data, Rlist, period, local_atpair, qlist, sinvS); } - } - else - { - if (Params::use_fullcoul_wc) + + Profiler::stop("chi0_build"); + + std::flush(ofs_myid); + mpi_comm_global_h.barrier(); + + if (Params::debug) + { // debug, check chi0 + char fn[80]; + for (const auto &chi0q : chi0.get_chi0_q()) + { + const int ifreq = chi0.tfg.get_freq_index(chi0q.first); + for (const auto &q_IJchi0 : chi0q.second) + { + const int iq = std::distance(klist.begin(), + std::find(klist.begin(), klist.end(), q_IJchi0.first)); + for (const auto &I_Jchi0 : q_IJchi0.second) + { + const auto &I = I_Jchi0.first; + for (const auto &J_chi0 : I_Jchi0.second) + { + const auto &J = J_chi0.first; + sprintf(fn, "chi0fq_ifreq_%d_iq_%d_I_%zu_J_%zu_id_%d.mtx", ifreq, iq, I, J, + mpi_comm_global_h.myid); + print_complex_matrix_mm(J_chi0.second, Params::output_dir + "/" + fn, + 1e-15); + } + } + } + } + } + + Profiler::start("g0w0_wc", "Build screened interaction"); + vector> epsmac_LF_imagfreq(epsmac_LF_imagfreq_re.cbegin(), + epsmac_LF_imagfreq_re.cend()); + map, matrix_m>>>::pair_t_old> + Wc_freq_q; + if (Params::use_scalapack_gw_wc) { - Wc_freq_q = compute_Wc_freq_q(chi0, Vq, Vq, epsmac_LF_imagfreq); + if (Params::use_fullcoul_wc) + { + Wc_freq_q = compute_Wc_freq_q_blacs(chi0, Vq, Vq, epsmac_LF_imagfreq); + } + else + { + Wc_freq_q = compute_Wc_freq_q_blacs(chi0, Vq, Vq_cut, epsmac_LF_imagfreq); + } } else { - Wc_freq_q = compute_Wc_freq_q(chi0, Vq, Vq_cut, epsmac_LF_imagfreq); + if (Params::use_fullcoul_wc) + { + Wc_freq_q = compute_Wc_freq_q(chi0, Vq, Vq, epsmac_LF_imagfreq); + } + else + { + Wc_freq_q = compute_Wc_freq_q(chi0, Vq, Vq_cut, epsmac_LF_imagfreq); + } } - } - Profiler::stop("g0w0_wc"); - if (Params::debug) - { // debug, check Wc - char fn[80]; - for (const auto &Wc : Wc_freq_q) - { - const int ifreq = chi0.tfg.get_freq_index(Wc.first); - for (const auto &I_JqWc : Wc.second) + Profiler::stop("g0w0_wc"); + if (Params::debug) + { // debug, check Wc(q, iw) + char fn[80]; + for (const auto &Wc : Wc_freq_q) { - const auto &I = I_JqWc.first; - for (const auto &J_qWc : I_JqWc.second) + const int ifreq = chi0.tfg.get_freq_index(Wc.first); + for (const auto &I_JqWc : Wc.second) { - const auto &J = J_qWc.first; - for (const auto &q_Wc : J_qWc.second) + const auto &I = I_JqWc.first; + for (const auto &J_qWc : I_JqWc.second) { - const int iq = std::distance( - klist.begin(), std::find(klist.begin(), klist.end(), q_Wc.first)); - sprintf(fn, "Wcfq_ifreq_%d_iq_%d_I_%zu_J_%zu_id_%d.mtx", ifreq, iq, I, J, - mpi_comm_global_h.myid); - print_matrix_mm_file(q_Wc.second, Params::output_dir + "/" + fn, 1e-15); + const auto &J = J_qWc.first; + for (const auto &q_Wc : J_qWc.second) + { + const int iq = std::distance( + klist.begin(), std::find(klist.begin(), klist.end(), q_Wc.first)); + sprintf(fn, "Wcfq_ifreq_%d_iq_%d_I_%zu_J_%zu_id_%d.mtx", ifreq, iq, I, J, + mpi_comm_global_h.myid); + print_matrix_mm_file(q_Wc.second, Params::output_dir + "/" + fn, 1e-15); + } } } } } - } - if (Params::use_shrink_abfs) - { - Profiler::start("read_shrink_sinvS_fold", "Load shrink transformation"); - // change atom_mu: number of {Mu,mu} in the later calculations - read_shrink_sinvS(driver_params.input_dir, "shrink_sinvS_", sinvS); - Profiler::stop("read_shrink_sinvS_fold"); - } - - LIBRPA::G0W0 s_g0w0(meanfield, kfrac_list, chi0.tfg, period); - Profiler::start("g0w0_sigc_IJ", "Build real-space correlation self-energy"); + if (Params::use_shrink_abfs) + { + Profiler::start("read_shrink_sinvS_fold", "Load shrink transformation"); + // change atom_mu: number of {Mu,mu} in the later calculations + read_shrink_sinvS(driver_params.input_dir, "shrink_sinvS_", sinvS); + Profiler::stop("read_shrink_sinvS_fold"); + } - if (Params::use_soc) - s_g0w0.build_spacetime>(Cs_data, Wc_freq_q, Rlist, qlist, sinvS); - else - s_g0w0.build_spacetime(Cs_data, Wc_freq_q, Rlist, qlist, sinvS); + Profiler::start("g0w0_sigc_IJ", "Build real-space correlation self-energy"); - Profiler::stop("g0w0_sigc_IJ"); - std::flush(ofs_myid); - - /* - * Compute the QP energies on k-grid - */ - Profiler::start("read_vxc", "Load DFT xc potential"); - std::vector vxc; - int flag_read_vxc = read_vxc(driver_params.input_dir + "vxc_out", vxc); - Profiler::stop("read_vxc"); + if (Params::use_soc) + s_g0w0.build_spacetime>(Cs_data, Wc_freq_q, Rlist, qlist, sinvS); + else + s_g0w0.build_spacetime(Cs_data, Wc_freq_q, Rlist, qlist, sinvS); - if (flag_read_vxc != 0) - { - if (mpi_comm_global_h.myid == 0) - { - lib_printf("Error in reading Vxc on kgrid, task failed!\n"); - } - Profiler::stop("g0w0_band"); - return; + Profiler::stop("g0w0_sigc_IJ"); + std::flush(ofs_myid); } + // ================== Compute the QP energies on k-grid ================== Profiler::start("g0w0_exx_ks_kgrid"); exx.build_KS_kgrid(); Profiler::stop("g0w0_exx_ks_kgrid"); @@ -256,7 +263,7 @@ void task_g0w0_band(std::map, ComplexMatrix> &sinvS) // imaginary freqencies for analytic continuation std::vector imagfreqs; - for (const auto &freq : chi0.tfg.get_freq_nodes()) + for (const auto &freq : tfg.get_freq_nodes()) { imagfreqs.push_back(cplxdb{0.0, freq}); } @@ -284,7 +291,7 @@ void task_g0w0_band(std::map, ComplexMatrix> &sinvS) const auto &exx_state = exx.Eexx[i_spin][i_kpoint][i_state]; const auto &vxc_state = vxc[i_spin](i_kpoint, i_state); std::vector sigc_state; - for (const auto &freq : chi0.tfg.get_freq_nodes()) + for (const auto &freq : tfg.get_freq_nodes()) { sigc_state.push_back(sigc_sk.at(freq)(i_state, i_state)); } @@ -342,7 +349,36 @@ void task_g0w0_band(std::map, ComplexMatrix> &sinvS) printf("\n"); } } - + if (Params::output_energy_qp) { + std::ofstream ofs("energy_qp"); + ofs << " state occ_num e_gs(Ha) e_qp(Ha)"<, ComplexMatrix> &sinvS) } Profiler::stop("g0w0_solve_qpe_kgrid"); - /* - * Compute the QP energies on band k-paths - */ + // ================== Compute the QP energies on k-path ================== // Reset k-space EXX and Sigmac matrices to avoid warning from internal reset exx.reset_kspace(); s_g0w0.reset_kspace(); @@ -472,7 +506,7 @@ void task_g0w0_band(std::map, ComplexMatrix> &sinvS) const auto &exx_state = exx.Eexx[i_spin][i_kpoint][i_state]; const auto &vxc_state = vxc_band[i_spin](i_kpoint, i_state); std::vector sigc_state; - for (const auto &freq : chi0.tfg.get_freq_nodes()) + for (const auto &freq : tfg.get_freq_nodes()) { sigc_state.push_back(sigc_sk.at(freq)(i_state, i_state)); } @@ -565,7 +599,7 @@ void task_g0w0_band(std::map, ComplexMatrix> &sinvS) << std::setprecision(7) << k.z; for (int i_state = 0; i_state < meanfield.get_n_bands(); i_state++) { - const auto &occ_state = mf.get_weight()[i_spin](i_kpoint, i_state); + const auto &occ_state = mf.get_weight()[i_spin](i_kpoint, i_state) * mf.get_n_kpoints(); const auto &eks_state = mf.get_eigenvals()[i_spin](i_kpoint, i_state) * HA2EV; const auto &exx_state = exx.Eexx[i_spin][i_kpoint][i_state] * HA2EV; const auto &vxc_state = vxc_band[i_spin](i_kpoint, i_state) * HA2EV; diff --git a/driver/task_screened_coulomb.cpp b/driver/task_screened_coulomb.cpp index 7ba66ba0..0a5c02ba 100644 --- a/driver/task_screened_coulomb.cpp +++ b/driver/task_screened_coulomb.cpp @@ -99,7 +99,7 @@ void task_screened_coulomb_real_freq(std::map, ComplexMatr if (iatom_row == iatom_col) continue; for (const auto &Wc_q : Wc.at(iatom_row).at(iatom_col)) { - Wc[iatom_col][iatom_row][Wc_q.first] = conj(Wc_q.second); + Wc[iatom_col][iatom_row][Wc_q.first] = Wc_q.second.get_transpose(true); } } } diff --git a/src/epsilon.cpp b/src/epsilon.cpp index bd03edb3..0d1c2b9a 100644 --- a/src/epsilon.cpp +++ b/src/epsilon.cpp @@ -2440,7 +2440,7 @@ FT_Wc_freq_q( const int ngrids = tfg.get_n_grids(); if (Params::debug) { - if (mpi_comm_global_h.is_root()) lib_printf("Converting Wc q,w -> R,t\n"); + if (mpi_comm_global_h.is_root()) lib_printf("Converting Wc q,w -> R,w\n"); mpi_comm_global_h.barrier(); } set> atpairs_unique; @@ -2673,8 +2673,240 @@ CT_FT_Wc_freq_q( return Wc_tau_R; } +map, matrix_m>>>::pair_t_old> +CT_FT_Wc_q2R_freq2time( + map, matrix_m>>>::pair_t_old> + &Wc_freq_q, + const TFGrids &tfg, const int &n_k_points, const vector> &Rlist) +{ + // major of Wc_freq_q input and Wc_tau_R output + const MAJOR major_Wc = MAJOR::ROW; + + map, matrix_m>>>::pair_t_old> + Wc_tau_R, Wc_freq_R; + if (!tfg.has_time_grids()) throw logic_error("TFGrids object does not have time grids"); + const int ngrids = tfg.get_n_grids(); + set> atpairs_unique; + + LIBRPA::utils::lib_printf_root("Converting Wc(q,w) -> Wc(R,w) -> W(R,t)\n"); + + Profiler::start("construct_Wc_lower_half", "Construct Lower Half of Wc(q,w)"); + // NOTE: only upper half of Wc is built now + // here we recover the other half before transform to R space using the Hermitian property + // and Hermitize the diagonal blocks (due to numerical noise) + for (int ifreq = 0; ifreq < ngrids; ifreq++) + { + const auto freq = tfg.get_freq_nodes()[ifreq]; + auto &Wc = Wc_freq_q.at(freq); + vector iatoms_row; + for (const auto &Mu_NuqWc : Wc) iatoms_row.push_back(Mu_NuqWc.first); + for (auto iatom_row: iatoms_row) + { + vector iatoms_col; + for (const auto &Nu_qWc : Wc.at(iatom_row)) + { + iatoms_col.push_back(Nu_qWc.first); + } + for (auto iatom_col : iatoms_col) + { + atpairs_unique.insert({iatom_row, iatom_col}); + atpairs_unique.insert({iatom_col, iatom_row}); + for (const auto &q_Wc : Wc.at(iatom_row).at(iatom_col)) + { + assert(q_Wc.second.major() == major_Wc); + if(iatom_row != iatom_col) + Wc[iatom_col][iatom_row][q_Wc.first] = q_Wc.second.get_transpose(true); + else // Hermitize the diagonal blocks + { + auto Wc_mat = q_Wc.second; + Wc_mat = (Wc_mat + Wc_mat.get_transpose(true)) * 0.5; + Wc[iatom_row][iatom_row][q_Wc.first] = Wc_mat; + } + } + } + } + } + mpi_comm_global_h.barrier(); + Profiler::stop("construct_Wc_lower_half"); + + Profiler::start("Wc(q,w) -> Wc(R,w)", "Convert Wc(q,w) -> Wc(R,w)"); + vector>, pair>> ifreqtau_R_atpair_all; + // allocate Wc(R,w) before hand + for (auto R : Rlist) + { + for (int ifreq = 0; ifreq != ngrids; ifreq++) + { + auto tau = tfg.get_time_nodes()[ifreq]; + auto freq = tfg.get_freq_nodes()[ifreq]; + for (auto atpair_unique : atpairs_unique) + { + const auto Mu = atpair_unique.first; + const int n_mu = atom_mu[Mu]; + const auto Nu = atpair_unique.second; + const int n_nu = atom_mu[Nu]; + Wc_freq_R[freq][Mu][Nu][R] = matrix_m>(n_mu, n_nu, major_Wc); + ifreqtau_R_atpair_all.push_back({{ifreq, R}, atpair_unique});// ifreq is same as itauR + } + } + } + + LIBRPA::utils::lib_printf_coll("Task %4d: distributing %d {I, J, R, freq} on %d threads\n", + LIBRPA::envs::myid_global, ifreqtau_R_atpair_all.size(), + omp_get_max_threads()); + +#pragma omp parallel for schedule(dynamic) + for (auto ifreqR_atpair : ifreqtau_R_atpair_all) + { + const auto ifreq = ifreqR_atpair.first.first; + const auto freq = tfg.get_freq_nodes()[ifreq]; + const auto R = ifreqR_atpair.first.second; + const auto Mu = ifreqR_atpair.second.first; + const auto Nu = ifreqR_atpair.second.second; + const int n_mu = atom_mu[Mu]; + const int n_nu = atom_mu[Nu]; + + // thread local temporary matrix + matrix_m> WfR_temp(n_mu, n_nu, major_Wc); + + if (Wc_freq_q.count(freq) == 0) continue; + if (Wc_freq_q.at(freq).count(Mu) == 0) continue; + if (Wc_freq_q.at(freq).at(Mu).count(Nu) == 0) continue; + + for (auto &Wc_q : Wc_freq_q.at(freq).at(Mu).at(Nu)) + { + const auto q = Wc_q.first; + const auto &Wc = Wc_q.second; + for (auto q_bz : map_irk_ks[q]) + { + const double ang = -q_bz * (R * latvec) * TWO_PI; + const complex weight = + complex(cos(ang), sin(ang)) / double(n_k_points); + if (q == q_bz) + WfR_temp += Wc * weight; + else + WfR_temp += conj(Wc) * weight; + } + } + // omp_set_lock(&lock_Wc); + Wc_freq_R[freq][Mu][Nu][R] += WfR_temp; + // omp_unset_lock(&lock_Wc); + } + mpi_comm_global_h.barrier(); + // HACK: Free up Wc_freq_q to save memory, especially for large Coulomb matrix case and many + // minimax grids + Wc_freq_q.clear(); + LIBRPA::utils::lib_printf_root("Done converting Wc(q,w) -> Wc(R,w)\n"); + + Profiler::stop("Wc(q,w) -> Wc(R,w)"); + + if (Params::output_Wc_Rf_mat > 0) + { + Profiler::start("write_Wc_freq_R", "Export Wc(R,w) to file"); + int write_freq = Params::output_Wc_Rf_mat==1 ? 1 : ngrids; + for (int ifreq = 0; ifreq != write_freq; ifreq++) + { + char fn[80]; + const auto freq = tfg.get_freq_nodes()[ifreq]; + auto& freq_MuNuRWc = Wc_freq_R.at(freq); + for (const auto &Mu_NuRWc : freq_MuNuRWc) + { + auto Mu = Mu_NuRWc.first; + // const int n_mu = atom_mu[Mu]; + for (const auto &Nu_RWc : Mu_NuRWc.second) + { + auto Nu = Nu_RWc.first; + // const int n_nu = atom_mu[Nu]; + for (const auto &R_Wc : Nu_RWc.second) + { + auto R = R_Wc.first; + auto Wc = R_Wc.second; + auto iteR = std::find(Rlist.cbegin(), Rlist.cend(), R); + auto iR = std::distance(Rlist.cbegin(), iteR); + sprintf(fn, "Wc_Mu_%zu_Nu_%zu_iR_%zu_ifreq_%d.mtx", Mu, Nu, iR, ifreq); + std::string info = "Wc at iR " + std::to_string(iR) + + " ( " + std::to_string(R.x) + " " + std::to_string(R.y) + " " + std::to_string(R.z) + + " ) and ifreq " + std::to_string(ifreq) + + " ( " + std::to_string(freq) + " a.u. )"; + print_matrix_mm_file(Wc, Params::output_dir + "/" + fn, info, 1e-10); + } + } + } + } + Profiler::stop("write_Wc_freq_R"); + } + + Profiler::start("Wc(R,w) -> Wc(R,t)", "Convert Wc(R,w) -> Wc(R,t)"); + if (mpi_comm_global_h.is_root()) + { + lib_printf("Start converting Wc(R,w) -> Wc(R,t)\n"); + } + // allocate Wc(R,t) before hand + for (auto R : Rlist) + { + for (int itau = 0; itau != ngrids; itau++) + { + auto tau = tfg.get_time_nodes()[itau]; + auto freq = tfg.get_freq_nodes()[itau]; + for (auto atpair_unique : atpairs_unique) + { + const auto Mu = atpair_unique.first; + const int n_mu = atom_mu[Mu]; + const auto Nu = atpair_unique.second; + const int n_nu = atom_mu[Nu]; + Wc_tau_R[tau][Mu][Nu][R] = matrix_m>(n_mu, n_nu, major_Wc); + } + } + } + + LIBRPA::utils::lib_printf_coll("Task %4d: distributing %d {I, J, R, tau} on %d threads\n", + LIBRPA::envs::myid_global, ifreqtau_R_atpair_all.size(), + omp_get_max_threads()); + +#pragma omp parallel for schedule(dynamic) + for (auto itauR_atpair : ifreqtau_R_atpair_all) + { + const auto itau = itauR_atpair.first.first; + const auto tau = tfg.get_time_nodes()[itau]; + const auto R = itauR_atpair.first.second; + const auto Mu = itauR_atpair.second.first; + const auto Nu = itauR_atpair.second.second; + const int n_mu = atom_mu[Mu]; + const int n_nu = atom_mu[Nu]; + + // thread local temporary matrix + matrix_m> WtR_temp(n_mu, n_nu, major_Wc); + + for (int ifreq = 0; ifreq < ngrids; ifreq++) + { + const auto freq = tfg.get_freq_nodes()[ifreq]; + const auto f2t = tfg.get_costrans_f2t()(itau, ifreq); + // ofs_myid << "f2t cos eff for freq " << freq << " -> tau " << tau << ": " << f2t << + // "\n"; + if (Wc_freq_R.count(freq) == 0) continue; + if (Wc_freq_R.at(freq).count(Mu) == 0) continue; + if (Wc_freq_R.at(freq).at(Mu).count(Nu) == 0) continue; + if (Wc_freq_R.at(freq).at(Mu).at(Nu).count(R) == 0) continue; + // cout << "freq: " << freq << "\n"; + + const auto &Wc = Wc_freq_R.at(freq).at(Mu).at(Nu).at(R); + WtR_temp += Wc * f2t; + } + // omp_set_lock(&lock_Wc); + Wc_tau_R[tau][Mu][Nu][R] += WtR_temp; + // omp_unset_lock(&lock_Wc); + } + // NOTE: Wc(R,w) will not be used any more, clean to free up memory + Wc_freq_R.clear(); + LIBRPA::utils::release_free_mem(); + mpi_comm_global_h.barrier(); + LIBRPA::utils::lib_printf_root("Done converting Wc(R,w) -> Wc(R,t)\n"); + Profiler::stop("Wc(R,w) -> Wc(R,t)"); + + return Wc_tau_R; +} + map, matrix_m>>>::pair_t_old> -CT_FT_Wc_freq2time_q( +CT_Wc_freq2time_q( const map, matrix_m>>>::pair_t_old> &Wc_freq_q, @@ -2776,25 +3008,22 @@ CT_FT_Wc_freq2time_q( return Wc_tau_q; } -atom_mapping, matrix_m>>>::pair_t_old CT_FT_Wc_tau_R2q( +/// @brief Wc(q,w) -> Wc(R,w) or Wc(q,t) -> W(R,t) +atom_mapping, matrix_m>>>::pair_t_old FT_Wc_q2R( const atom_mapping, matrix_m>>>::pair_t_old - &Wc_tau_q, - const TFGrids &tfg, const int &n_kpoints, const vector> &Rlist, - const int &itau) + &Wc_q, + const TFGrids &tfg, const int &n_kpoints, const vector> &Rlist, const bool is_freq) { - const int tau = tfg.get_time_nodes()[itau]; // major of Wc_freq_q input and Wc_tau_R output const MAJOR major_Wc = MAJOR::ROW; - atom_mapping, matrix_m>>>::pair_t_old Wc_tau_R; - if (!tfg.has_time_grids()) throw logic_error("TFGrids object does not have time grids"); - const int ngrids = tfg.get_n_grids(); + atom_mapping, matrix_m>>>::pair_t_old Wc_R; - LIBRPA::utils::lib_printf_root("Converting Wc(q,t) -> W(R,t)\n"); + LIBRPA::utils::lib_printf_root("Converting Wc(q) -> W(R)\n"); mpi_comm_global_h.barrier(); set> atpairs_unique; - for (const auto &MuNuqWc : Wc_tau_q) + for (const auto &MuNuqWc : Wc_q) { const auto Mu = MuNuqWc.first; for (const auto &Nu_qWc : MuNuqWc.second) @@ -2818,12 +3047,12 @@ atom_mapping, matrix_m>>>::pair_t_ol const int n_mu = atom_mu_l[Mu]; const auto Nu = atpair_unique.second; const int n_nu = atom_mu_l[Nu]; - Wc_tau_R[Mu][Nu][R] = matrix_m>(n_mu, n_nu, major_Wc); + Wc_R[Mu][Nu][R] = matrix_m>(n_mu, n_nu, major_Wc); iR_atpair_all.push_back({R, atpair_unique}); } } - LIBRPA::utils::lib_printf_coll("Task %4d: distributing %d {I, J, R, tau} on %d threads\n", + LIBRPA::utils::lib_printf_coll("Task %4d: distributing %d {I, J, R} on %d threads\n", LIBRPA::envs::myid_global, iR_atpair_all.size(), omp_get_max_threads()); @@ -2837,13 +3066,12 @@ atom_mapping, matrix_m>>>::pair_t_ol const int n_nu = atom_mu_l[Nu]; // thread local temporary matrix - matrix_m> WtR_temp(n_mu, n_nu, major_Wc); + matrix_m> WR_temp(n_mu, n_nu, major_Wc); - if (Wc_tau_q.count(Mu) == 0) continue; - if (Wc_tau_q.at(Mu).count(Nu) == 0) continue; - // cout << "freq: " << freq << "\n"; + if (Wc_q.count(Mu) == 0) continue; + if (Wc_q.at(Mu).count(Nu) == 0) continue; - const auto &Wc_q_all = Wc_tau_q.at(Mu).at(Nu); + const auto &Wc_q_all = Wc_q.at(Mu).at(Nu); for (auto &Wc_q : Wc_q_all) { const auto q = Wc_q.first; @@ -2854,19 +3082,51 @@ atom_mapping, matrix_m>>>::pair_t_ol const complex weight = complex(cos(ang), sin(ang)) / double(n_kpoints); if (q == q_bz) - WtR_temp += Wc * weight; + WR_temp += Wc * weight; else - WtR_temp += conj(Wc) * weight; + WR_temp += conj(Wc) * weight; } } // omp_set_lock(&lock_Wc); - Wc_tau_R[Mu][Nu][R] += WtR_temp; + Wc_R[Mu][Nu][R] += WR_temp; // omp_unset_lock(&lock_Wc); } - - LIBRPA::utils::lib_printf_root("Done converting Wc q,t -> R,t\n"); mpi_comm_global_h.barrier(); - return Wc_tau_R; + LIBRPA::utils::lib_printf_root("Done converting Wc q -> R\n"); + + if (is_freq && Params::output_Wc_Rf_mat==1) + { + Profiler::start("write_Wc_freq_R", "Export Wc(R,w) to file"); + int ifreq = 0; + const auto freq = tfg.get_freq_nodes()[ifreq]; + char fn[80]; + for (const auto &Mu_NuRWc : Wc_R) + { + auto Mu = Mu_NuRWc.first; + // const int n_mu = atom_mu[Mu]; + for (const auto &Nu_RWc : Mu_NuRWc.second) + { + auto Nu = Nu_RWc.first; + // const int n_nu = atom_mu[Nu]; + for (const auto &R_Wc : Nu_RWc.second) + { + auto R = R_Wc.first; + auto Wc = R_Wc.second; + auto iteR = std::find(Rlist.cbegin(), Rlist.cend(), R); + auto iR = std::distance(Rlist.cbegin(), iteR); + sprintf(fn, "Wc_Mu_%zu_Nu_%zu_iR_%zu_ifreq_%d.mtx", Mu, Nu, iR, ifreq); + std::string info = "Wc at iR " + std::to_string(iR) + + " ( " + std::to_string(R.x) + " " + std::to_string(R.y) + " " + std::to_string(R.z) + + " ) and ifreq " + std::to_string(ifreq) + + " ( " + std::to_string(freq) + " a.u. )"; + print_matrix_mm_file(Wc, Params::output_dir + "/" + fn, info, 1e-10); + } + } + } + Profiler::stop("write_Wc_freq_R"); + } + + return Wc_R; } void test_libcomm_for_system(const atpair_k_cplx_mat_t &coulmat) diff --git a/src/epsilon.h b/src/epsilon.h index 656d3926..11dd0f13 100644 --- a/src/epsilon.h +++ b/src/epsilon.h @@ -58,19 +58,29 @@ CT_FT_Wc_freq_q( &Wc_freq_q, const TFGrids &tfg, const int &n_kpoints, const vector> &Rlist); +/// @brief Fourier transform screened Coulomb Wc(q,w) -> Wc(R,w) -> W(R,t) +/// @details transform step by step to output Wc_freq_R, and return full atom-pair matrix +/// @attention CT_FT_Wc_freq_q only return upper atom-pair, but final Wc_libri should be same +map, matrix_m>>>::pair_t_old> +CT_FT_Wc_q2R_freq2time( + map, matrix_m>>>::pair_t_old> + &Wc_freq_q, // upper atom-pair input +const TFGrids &tfg, const int &n_kpoints, const vector> &Rlist); + +/// @brief Wc(q,w) -> Wc(q,t) map, matrix_m>>>::pair_t_old> -CT_FT_Wc_freq2time_q( +CT_Wc_freq2time_q( const map, matrix_m>>>::pair_t_old> &Wc_freq_q, const TFGrids &tfg, const int &n_kpoints, const vector> &Rlist, const vector> &qlist); -atom_mapping, matrix_m>>>::pair_t_old CT_FT_Wc_tau_R2q( +/// @brief Wc(q,w) -> Wc(R,w) or Wc(q,t) -> W(R,t) +atom_mapping, matrix_m>>>::pair_t_old FT_Wc_q2R( const atom_mapping, matrix_m>>>::pair_t_old &Wc_tau_q, - const TFGrids &tfg, const int &n_kpoints, const vector> &Rlist, - const int &itau); + const TFGrids &tfg, const int &n_kpoints, const vector> &Rlist, const bool is_freq); ComplexMatrix compute_Pi_freq_q_row_ri(const Vector3_Order &ik_vec, const atom_mapping::pair_t_old &chi0_freq_q, diff --git a/src/gw.cpp b/src/gw.cpp index 073a6f95..c432d8b2 100644 --- a/src/gw.cpp +++ b/src/gw.cpp @@ -324,8 +324,55 @@ void G0W0::build_spacetime( // Transform if (Params::use_shrink_abfs) { + if (Params::output_Wc_Rf_mat == 1) + { + Profiler::start("unfold_Wc_q", "Unfold Wc (q,w=0)"); + const double freq = tfg.get_freq_nodes()[0]; + auto Wc_q_f0 = Wc_freq_q.at(freq); + unfold_abfs_Wc(sinvS, Wc_q_f0, qlist, atom_mu_l, atom_mu_s); + Profiler::stop("unfold_Wc_q"); + Profiler::start("construct_Wc_lower_half", "Construct Lower Half of Wc(q,w)"); + // NOTE: only upper half of Wc is built now + // here we recover the other half before transform to R space using the Hermitian property + // and Hermitize the diagonal blocks (due to numerical noise) + vector iatoms_row; + for (const auto &Mu_NuqWc : Wc_q_f0) iatoms_row.push_back(Mu_NuqWc.first); + for (auto iatom_row: iatoms_row) + { + vector iatoms_col; + for (const auto &Nu_qWc : Wc_q_f0.at(iatom_row)) + { + iatoms_col.push_back(Nu_qWc.first); + } + for (auto iatom_col : iatoms_col) + { + for (const auto &q_Wc : Wc_q_f0.at(iatom_row).at(iatom_col)) + { + assert(q_Wc.second.major() == MAJOR::ROW); + if(iatom_row != iatom_col) + Wc_q_f0[iatom_col][iatom_row][q_Wc.first] = q_Wc.second.get_transpose(true); + else // Hermitize the diagonal blocks + { + auto Wc_mat = q_Wc.second; + Wc_mat = (Wc_mat + Wc_mat.get_transpose(true)) * 0.5; + Wc_q_f0[iatom_row][iatom_row][q_Wc.first] = Wc_mat; + } + } + } + } + + mpi_comm_global_h.barrier(); + Profiler::stop("construct_Wc_lower_half"); + FT_Wc_q2R(Wc_q_f0, tfg, meanfield.get_n_kpoints(), Rlist, true); + Wc_q_f0.clear(); + } + else if (Params::output_Wc_Rf_mat > 1) + { + throw std::logic_error("use_shrink_abfs doesn't support output_Wc_Rf_mat > 1"); + } Profiler::start("g0w0_build_spacetime_wt_ft_wc", "Tranform Wc (q,w) -> (q,t)"); - Wc_tau_q = CT_FT_Wc_freq2time_q(Wc_freq_q, tfg, meanfield.get_n_kpoints(), Rlist, qlist); + Wc_tau_q = CT_Wc_freq2time_q(Wc_freq_q, tfg, meanfield.get_n_kpoints(), Rlist, qlist); + // HACK: Free up Wc_freq_q to save memory, especially for large Coulomb matrix case and many // minimax grids Wc_freq_q.clear(); @@ -339,11 +386,18 @@ void G0W0::build_spacetime( else { Profiler::start("g0w0_build_spacetime_ct_ft_wc", "Tranform Wc (q,w) -> (R,t)"); - Wc_tau_R = CT_FT_Wc_freq_q(Wc_freq_q, tfg, meanfield.get_n_kpoints(), Rlist); - // HACK: Free up Wc_freq_q to save memory, especially for large Coulomb matrix case and many - // minimax grids - Wc_freq_q.clear(); - utils::release_free_mem(); + if (Params::output_Wc_Rf_mat > 0) + { + Wc_tau_R = CT_FT_Wc_q2R_freq2time(Wc_freq_q, tfg, meanfield.get_n_kpoints(), Rlist); + } + else + { + Wc_tau_R = CT_FT_Wc_freq_q(Wc_freq_q, tfg, meanfield.get_n_kpoints(), Rlist); + // HACK: Free up Wc_freq_q to save memory, especially for large Coulomb matrix case and many + // minimax grids + Wc_freq_q.clear(); + utils::release_free_mem(); + } Profiler::stop("g0w0_build_spacetime_ct_ft_wc"); LIBRPA::utils::lib_printf_root( "Time for Fourier transform of Wc in GW (seconds, Wall/CPU): %f %f\n", @@ -404,7 +458,7 @@ void G0W0::build_spacetime( Profiler::stop("unfold_Wc_abfs"); Profiler::start("g0w0_build_spacetime_Rq_ft_wc", "Tranform Wc (q,t) -> (R,t)"); Wc_tau_R[tau] = - CT_FT_Wc_tau_R2q(Wc_tau_q[tau], tfg, meanfield.get_n_kpoints(), Rlist, itau); + FT_Wc_q2R(Wc_tau_q[tau], tfg, meanfield.get_n_kpoints(), Rlist, false); Profiler::stop("g0w0_build_spacetime_Rq_ft_wc"); Wc_tau_q[tau].clear(); atom_mu = atom_mu_l; @@ -446,9 +500,11 @@ void G0W0::build_spacetime( else Wc_libri[static_cast(I)][{static_cast(J), {R.x, R.y, R.z}}] = RI::Tensor({nabf_I, nabf_J}, R_Wc.second.get_real().sptr()); - // cout << "I " << I << " J " << J << " R " << R << " tau " << tau << endl - // ; cout << Wc_libri[I][{J, {R.x, R.y, R.z}}] << endl; handle the - // block + // cout << "I " << I << " J " << J << " R " << R << " tau " << tau << endl; + // cout << Wc_libri[static_cast(I)][{static_cast(J), {R.x, R.y, R.z}}] << endl; + // std::cout << "R_Wc.second: " << R_Wc.second << std::endl; + // handle the block + if (Params::output_Wc_Rf_mat > 0 && !Params::use_shrink_abfs) continue; // full atom-pair has been constructed if (I == J) continue; auto minusR = (-R) % this->period_; if (J_RWc.second.count(minusR) == 0) continue; @@ -619,6 +675,7 @@ void G0W0::build_spacetime( std::shared_ptr> mat_ptr = std::make_shared>( gf_IJ_block.c, gf_IJ_block.size); + tau_gf_libri[t][static_cast(I)] [{static_cast(J), {R.x, R.y, R.z}}] = RI::Tensor({n_I, n_J}, mat_ptr); @@ -915,6 +972,60 @@ void G0W0::build_spacetime( } } +void G0W0::read_sigc(const std::string &input_dir, const vector> &Rlist) +{ + using LIBRPA::envs::mpi_comm_global_h; + + LIBRPA::utils::lib_printf_root("Reading real-space imaginary-frequency NAO sigma_c matrices\n"); + Profiler::start("g0w0_read_sigc(R,iw) in NAO"); + for (int ispin = 0; ispin != mf.get_n_spins(); ispin++) + { + for (int isoc1 = 0; isoc1 != mf.get_n_soc(); isoc1++) + { + for (int isoc2 = 0; isoc2 != mf.get_n_soc(); isoc2++) + { + for (auto iomega = 0; iomega != tfg.get_n_grids(); iomega++) + { + size_t n_IJR_myid = 0; + std::stringstream ss; + ss << input_dir << "SigcRF_ispin_" << std::setfill('0') << std::setw(2) << ispin + << "_s_" << std::setw(1) << isoc1 << std::setw(1) << isoc2 << "_iomega_" + << std::setfill('0') << std::setw(3) << iomega << "_myid_" + << std::setfill('0') << std::setw(5) << envs::myid_global << ".dat"; + std::ifstream ifs_sigmac_r(ss.str(), std::ios::in | std::ios::binary); + if (!ifs_sigmac_r.is_open()) + { + throw std::runtime_error("Cannot open file " + ss.str()); + } + ifs_sigmac_r.read((char *)&n_IJR_myid, sizeof(size_t)); + + const auto omega = tfg.get_freq_nodes()[iomega]; + for (size_t idx = 0; idx != n_IJR_myid; idx++) + { + size_t dims[5]; + ifs_sigmac_r.read((char *)dims, 5 * sizeof(size_t)); + const auto iR = dims[0]; + const auto I = dims[1]; + const auto J = dims[2]; + const auto n_I = dims[3]; + const auto n_J = dims[4]; + Matz sigc(n_I, n_J, MAJOR::ROW); + ifs_sigmac_r.read((char *)sigc.ptr(), n_I * n_J * 2 * sizeof(double)); + const auto R = Rlist[iR]; + sigc_is_f_R_IJ[ispin][isoc1][isoc2][omega][R][I][J] = std::move(sigc); + } + ifs_sigmac_r.close(); + } + } + } + } + + is_rspace_built_ = true; + mpi_comm_global_h.barrier(); + LIBRPA::utils::lib_printf_root("Finished reading real-space imaginary-frequency NAO sigma_c matrices.\n"); + Profiler::stop("g0w0_read_sigc(R,iw) in NAO"); +} + void G0W0::build_sigc_matrix_KS( const std::vector>> &wfc_target, const std::vector> &kfrac_target) diff --git a/src/gw.h b/src/gw.h index 20947dee..7355e39b 100644 --- a/src/gw.h +++ b/src/gw.h @@ -29,7 +29,7 @@ class G0W0 // std::map, // atom_mapping::pair_t_old>>> sigc_is_f_k_IJ; - //! frequency-domain reciprocal-space correlation self-energy, indices + //! frequency-domain real-space correlation self-energy, indices //! [ispin][isoc1][isoc2][freq][R][I][J](n_I, n_J) std::map< int, @@ -74,10 +74,12 @@ class G0W0 const vector> &Rlist, const vector> &qlist, std::map, ComplexMatrix> &sinvS); + //! read the real-space correlation self-energy matrix on imaginary frequencies from disk + void read_sigc(const std::string &input_dir, const vector> &Rlist); //! build the correlation self-energy matrix in Kohn-Sham basis at the SCF k-points void build_sigc_matrix_KS_kgrid(); void build_sigc_matrix_KS_kgrid0(); - //! build the correlation self-energy matrix in Kohn-Sham basis at the SCF k-points + //! build the correlation self-energy matrix in Kohn-Sham basis at the band k-points void build_sigc_matrix_KS_band( const std::vector>> &wfc_target, const std::vector> &kfrac_band); diff --git a/src/matrix_m.h b/src/matrix_m.h index 2042a464..186a3667 100644 --- a/src/matrix_m.h +++ b/src/matrix_m.h @@ -1204,6 +1204,18 @@ void print_matrix_mm_file(const matrix_m &mat, const std::string &fn, Treal t fs.close(); } +template ::type> +void print_matrix_mm_file(const matrix_m &mat, const std::string &fn, const std::string& header_comment, Treal threshold = 1e-15, bool row_first = true) +{ + ofstream fs; + fs.open(fn); + if (!header_comment.empty()){ + fs << header_comment << std::endl; + } + print_matrix_mm(mat, fs, threshold, row_first); + fs.close(); +} + //! print the matrix in ELSI CSC format template ::type> void print_matrix_elsi_csc(const matrix_m &mat, const std::string &fn, Treal threshold = 1e-15) diff --git a/src/params.cpp b/src/params.cpp index c01e0b08..05b22355 100644 --- a/src/params.cpp +++ b/src/params.cpp @@ -44,9 +44,13 @@ bool Params::use_shrink_chi = true; bool Params::use_shrink_abfs = false; bool Params::use_soc = false; +bool Params::band_continue = false; + /* ========================================================== * output options begin * ========================================================== */ +int Params::output_Wc_Rf_mat = 0; +bool Params::output_energy_qp = false; bool Params::output_gw_sigc_mat = false; bool Params::output_gw_sigc_mat_rt = false; bool Params::output_gw_sigc_mat_rf = false; @@ -87,6 +91,7 @@ void Params::print() {"nfreq", nfreq}, {"n_params_anacon", n_params_anacon}, {"option_dielect_func", option_dielect_func}, + {"output_Wc_Rf_mat", output_Wc_Rf_mat}, {"nbands_G", nbands_G}, }; @@ -100,9 +105,13 @@ void Params::print() const std::vector> bool_params{ {"debug", debug}, + {"band_continue", band_continue}, {"use_scalapack_ecrpa", use_scalapack_ecrpa}, {"use_scalapack_gw_wc", use_scalapack_gw_wc}, + {"output_energy_qp", output_energy_qp}, {"output_gw_sigc_mat", output_gw_sigc_mat}, + {"output_gw_sigc_mat_rt", output_gw_sigc_mat_rt}, + {"output_gw_sigc_mat_rf", output_gw_sigc_mat_rf}, {"replace_w_head", replace_w_head}, {"use_shrink_abfs", use_shrink_abfs}, {"use_shrink_chi", use_shrink_chi}, diff --git a/src/params.h b/src/params.h index 0d48aa6c..83a01ff8 100644 --- a/src/params.h +++ b/src/params.h @@ -104,6 +104,9 @@ struct Params //! switch of using spin-orbit coupling correction static bool use_soc; + //! in task "g0w0_band", continue from previous self-energy matrix in NAO (R, iw) + static bool band_continue; + //! option of computing dielectric function on imaginary axis /*! * Available values: @@ -116,6 +119,18 @@ struct Params /* ========================================================== * output options */ + /* + switch of outputting Wc matrix in Abs (real space, imaginary frequency domain) + Available values: + - 0: do not output + - 1: output lowerest frequency + - 2: output all frequencies + */ + static int output_Wc_Rf_mat; + + //! output energy_qp file for BSE calculation outside + static bool output_energy_qp; + //! output correlation self-energy matrix (reciprocal space, imaginary frequency domain) static bool output_gw_sigc_mat;