From 5477d6c4b16b22604fe33ff7902356fcf02c6cf4 Mon Sep 17 00:00:00 2001 From: memii-x Date: Thu, 15 Oct 2026 14:02:46 +0800 Subject: [PATCH] fix(rt): convert sHG density to ordinary NAO frame for RI exchange --- source/source_esolver/esolver_ks_lcao.cpp | 8 +- .../source_estate/module_pot/potential_new.h | 3 + source/source_hamilt/module_xc/exx_info.cpp | 5 + .../module_parameter/read_inp_exx_dftu.cpp | 3 +- .../source_lcao/module_ri/exx_lri_interface.h | 19 +- .../module_ri/exx_lri_interface.hpp | 127 ++++++--- .../module_ri/hybrid_gauge_phase.h | 195 +++++++++++++ .../module_ri/unittests/CMakeLists.txt | 12 + .../module_ri/unittests/test_exx_info.cpp | 61 ++++ .../unittests/test_hybrid_gauge_phase.cpp | 266 ++++++++++++++++++ 10 files changed, 655 insertions(+), 44 deletions(-) create mode 100644 source/source_lcao/module_ri/hybrid_gauge_phase.h create mode 100644 source/source_lcao/module_ri/unittests/test_exx_info.cpp create mode 100644 source/source_lcao/module_ri/unittests/test_hybrid_gauge_phase.cpp diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 76a9bf3b9f..ef85d6b686 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -431,13 +431,17 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const // calculate exact-exchange if (this->inp_->calculation != "nscf") { + const auto* td_info = this->get_td_info(); + const auto* field = td_info != nullptr ? &td_info->field() : nullptr; if (exx_info_.info_ri.real_number) { - this->exx_nao.exd->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, this->pv, iter); + this->exx_nao.exd->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, this->pv, iter, + field); } else { - this->exx_nao.exc->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, this->pv, iter); + this->exx_nao.exc->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, this->pv, iter, + field); } } #endif diff --git a/source/source_estate/module_pot/potential_new.h b/source/source_estate/module_pot/potential_new.h index e1af826981..dc1825af21 100644 --- a/source/source_estate/module_pot/potential_new.h +++ b/source/source_estate/module_pot/potential_new.h @@ -80,6 +80,9 @@ class Potential : public PotBase PotBase* get_pot_type(const std::string& pot_type); + /// @brief Read the shared RT-TDDFT field state. + const TDFieldManager* td_field_manager() const { return td_field_manager_.get(); } + /** * @brief Inject the shared RT-TDDFT field state before potential setup. * diff --git a/source/source_hamilt/module_xc/exx_info.cpp b/source/source_hamilt/module_xc/exx_info.cpp index 7098c1dd7d..53355a08a0 100644 --- a/source/source_hamilt/module_xc/exx_info.cpp +++ b/source/source_hamilt/module_xc/exx_info.cpp @@ -218,6 +218,11 @@ bool init_exx_info(Exx_Info& exx_info, const Input_para& inp) { exx_info.info_lip.lambda = std::stod(inp.exx_fock_lambda[0]); } exx_info.info_ri.real_number = std::stoi(inp.exx_real_number); + if (exx_info.info_global.cal_exx && inp.basis_type == "lcao" + && inp.esolver_type == "tddft" && inp.td_stype == 2 && exx_info.info_ri.real_number) + { + throw std::invalid_argument("td_stype=2 EXX requires exx_real_number=0"); + } exx_info.info_ri.pca_threshold = inp.exx_pca_threshold; exx_info.info_ri.C_threshold = inp.exx_c_threshold; exx_info.info_ri.V_threshold = inp.exx_v_threshold; diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp index bbecf5cace..f978fecc7a 100644 --- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp @@ -444,7 +444,8 @@ void ReadInput::item_exx() item.category = "Exact Exchange (LCAO)"; item.type = "String"; item.description = R"(* True: Enforce LibRI to use double data type. -* False: Enforce LibRI to use complex data type. Setting it to True can effectively improve the speed of self-consistent calculations with hybrid functionals.)"; +* False: Enforce LibRI to use complex data type. Setting it to True can effectively improve the speed of self-consistent calculations with hybrid functionals. +* LCAO RT-TDDFT exact exchange with td_stype=2 requires exx_real_number=0 to represent the input density phase.)"; item.default_value = "depends on the gamma_only option"; item.unit = ""; read_sync_string(input.exx_real_number); diff --git a/source/source_lcao/module_ri/exx_lri_interface.h b/source/source_lcao/module_ri/exx_lri_interface.h index 8b5bc61cf6..dd814bdc9b 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.h +++ b/source/source_lcao/module_ri/exx_lri_interface.h @@ -4,6 +4,7 @@ #include "exx_lri.h" #include "source_estate/module_charge/chg_mix.h" #include "source_lcao/module_ri/mix_dmk_2d.h" +#include "source_lcao/module_ri/hybrid_gauge_phase.h" #include "source_lcao/module_ri/module_exx_symmetry/symm_rotation.h" #include "source_estate/module_dm/density_matrix.h" // mohan add 2025-11-04 #include "source_hamilt/hamilt.h" @@ -15,6 +16,7 @@ class Charge_Mixing; namespace elecstate { class ElecState; + class TDFieldManager; } namespace module_dm { @@ -105,7 +107,8 @@ class Exx_LRI_Interface const module_dm::DensityMatrix& dm/**< double should be Tdata if complex-PBE-DM is supported*/, const K_Vectors& kv, const Parallel_Orbitals& pv, - const int& iter); + const int& iter, + const elecstate::TDFieldManager* field); /// @brief in hamilt2rho: calculate Hexx and Eexx void exx_hamilt2rho(elecstate::ElecState& elec, const Parallel_Orbitals& pv, const int iter); @@ -132,7 +135,8 @@ class Exx_LRI_Interface int& iter, const int& istep, const double& etot, - const double& scf_ene_thr); + const double& scf_ene_thr, + const elecstate::TDFieldManager* field); /// @brief the step of the outer loop. /// nullptr: no dependence on the number of two_level_step, contributeHk will do enerything normally. @@ -150,6 +154,17 @@ class Exx_LRI_Interface private: Mix_DMk_2D mix_DMk_2D; + Mix_DMk_2D> mix_common_density_; + RI_2D_Comm::HybridDensityLayout common_density_layout_; + + std::vector>>> prepare_exx_density( + const std::vector>& dmk, const UnitCell& ucell, + const K_Vectors& kv, const Parallel_Orbitals& pv, const int nspin, + const bool restart, const ModuleBase::Vector3* source_A); + void mix_common_density(RI_2D_Comm::HybridGaugeComplexDensity& density, + const UnitCell& ucell, const ModuleBase::Vector3& source_A, const bool restart); + void mix_common_density(std::vector>>>& density, + const UnitCell& ucell, const ModuleBase::Vector3& source_A, const bool restart); Exx_Info_Global info_global; size_t hybrid_step_ = 1; diff --git a/source/source_lcao/module_ri/exx_lri_interface.hpp b/source/source_lcao/module_ri/exx_lri_interface.hpp index b08514faad..c3a7efb1ac 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.hpp +++ b/source/source_lcao/module_ri/exx_lri_interface.hpp @@ -1,6 +1,7 @@ #ifndef EXX_LRI_INTERFACE_HPP #define EXX_LRI_INTERFACE_HPP #include "exx_lri_interface.h" +#include "source_estate/module_pot/td_field_manager.h" #include "source_cell/unitcell.h" #include "source_base/formatter.h" #include "source_base/parallel_common.h" @@ -79,6 +80,64 @@ void Exx_LRI_Interface::cal_exx_stress(const double& omega, const doub this->flag_finish.stress = true; } +template +void Exx_LRI_Interface::mix_common_density( + RI_2D_Comm::HybridGaugeComplexDensity& density, + const UnitCell& ucell, + const ModuleBase::Vector3& source_A, + const bool restart) +{ + RI_2D_Comm::mix_hybrid_density(ucell.latvec, ucell.lat0, source_A, restart, + this->mix_common_density_, this->common_density_layout_, density); +} + +template +void Exx_LRI_Interface::mix_common_density( + std::vector>>>& density, + const UnitCell& ucell, + const ModuleBase::Vector3& source_A, + const bool restart) +{ + throw std::invalid_argument("sHG EXX density requires complex RI tensors (exx_real_number 0)"); +} + +template +std::vector::TA, + std::map::TAC, RI::Tensor>>> +Exx_LRI_Interface::prepare_exx_density( + const std::vector>& dmk, + const UnitCell& ucell, + const K_Vectors& kv, + const Parallel_Orbitals& pv, + const int nspin, + const bool restart, + const ModuleBase::Vector3* source_A) +{ + std::vector> restored_storage; + const std::vector>* restored = &dmk; + if (this->exx_spacegroup_symmetry) + { + restored_storage = this->symrot_.restore_dm(kv, dmk, pv); + restored = &restored_storage; + } + if (source_A == nullptr) + { + this->mix_DMk_2D.mix(*restored, restart); + return RI_2D_Comm::split_m2D_ktoR( + ucell, *this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_out(), pv, nspin, + this->exx_spacegroup_symmetry); + } + std::vector*> raw; + for (const auto& block : *restored) + { + raw.push_back(&block); + } + auto density = RI_2D_Comm::split_m2D_ktoR( + ucell, *this->exx_ptr->p_kv, raw, pv, nspin, this->exx_spacegroup_symmetry); + this->mix_common_density(density, ucell, *source_A, restart); + return density; +} + template void Exx_LRI_Interface::exx_before_all_runners( const K_Vectors& kv, @@ -163,9 +222,11 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, { this->mix_DMk_2D.set_nks(kv.get_nks()); } if (this->info_global.separate_loop) - { this->mix_DMk_2D.set_mixing_plain(this->info_global.mixing_beta_for_loop1); } + { this->mix_DMk_2D.set_mixing_plain(this->info_global.mixing_beta_for_loop1); + this->mix_common_density_.set_mixing_plain(this->info_global.mixing_beta_for_loop1); } else - { this->mix_DMk_2D.set_mixing(chgmix.get_mixing()); } + { this->mix_DMk_2D.set_mixing(chgmix.get_mixing()); + this->mix_common_density_.set_mixing(chgmix.get_mixing()); } // remember chgmix so exx_eachiterinit can re-borrow its mixing pointer after // mixing_restart's init_mixing() has reallocated it (else use-after-free -> SIGSEGV) @@ -183,7 +244,8 @@ void Exx_LRI_Interface::exx_eachiterinit(const int istep, const module_dm::DensityMatrix& dm, const K_Vectors& kv, const Parallel_Orbitals& pv, - const int& iter) + const int& iter, + const elecstate::TDFieldManager* field) { ModuleBase::TITLE("Exx_LRI_Interface","exx_eachiterinit"); if (this->info_global.cal_exx) @@ -212,32 +274,23 @@ void Exx_LRI_Interface::exx_eachiterinit(const int istep, { flag_restart = true; } this->last_borrowed_mixing_ = cur_mixing; this->mix_DMk_2D.set_mixing(this->p_chgmix_->get_mixing()); + this->mix_common_density_.set_mixing(this->p_chgmix_->get_mixing()); } - auto cal = [this, &ucell,&kv, &pv, &flag_restart](const module_dm::DensityMatrix& dm_in) + const auto density_time = RI_2D_Comm::hybrid_density_time(istep, iter); + const bool use_left_density = density_time == RI_2D_Comm::HybridDensityTime::Left; + const auto& source_density = use_left_density ? *this->dm_last_step : dm; + const ModuleBase::Vector3* source_A = nullptr; + if (field != nullptr && field->gauge() == 2) { - if (this->exx_spacegroup_symmetry) - { this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm_in.get_dmk_vec(), pv), flag_restart); } - else - { this->mix_DMk_2D.mix(dm_in.get_dmk_vec(), flag_restart); } - const std::vector>>> Ds = - RI_2D_Comm::split_m2D_ktoR( - ucell, - *this->exx_ptr->p_kv, - this->mix_DMk_2D.get_DMk_out(), - pv, - PARAM.inp.nspin, - this->exx_spacegroup_symmetry); - if(this->exx_spacegroup_symmetry && this->exx_ptr->info.exx_symmetry_realspace) - { this->cal_exx_elec(Ds, ucell, pv, &this->symrot_); } - else - { this->cal_exx_elec(Ds, ucell, pv); } - }; - - if(istep > 0 && flag_restart) - { cal(*this->dm_last_step); } + source_A = use_left_density ? &field->A_left_ha() : &field->A_right_ha(); + } + const auto Ds = this->prepare_exx_density(source_density.get_dmk_vec(), ucell, kv, pv, + PARAM.inp.nspin, flag_restart, source_A); + if (this->exx_spacegroup_symmetry && this->exx_ptr->info.exx_symmetry_realspace) + { this->cal_exx_elec(Ds, ucell, pv, &this->symrot_); } else - { cal(dm); } + { this->cal_exx_elec(Ds, ucell, pv); } } } } @@ -341,7 +394,7 @@ void Exx_LRI_Interface::exx_iter_finish(const K_Vectors& kv, iter, istep, elec.f_en.etot, - scf_ene_thr); + scf_ene_thr, elec.pot->td_field_manager()); // On the 0->1 transition, exx_after_converge switches the XC functional (PBE->hybrid), // but v_eff was already set for this (converged) iteration under the OLD functional. @@ -371,7 +424,8 @@ bool Exx_LRI_Interface::exx_after_converge( int& iter, const int& istep, const double& etot, - const double& scf_ene_thr) + const double& scf_ene_thr, + const elecstate::TDFieldManager* field) { // only called if (this->info_global.cal_exx) ModuleBase::TITLE("Exx_LRI_Interface","exx_after_converge"); auto restart_reset = [this]() @@ -432,18 +486,13 @@ bool Exx_LRI_Interface::exx_after_converge( // if init_wfc == "file", DM is calculated in the 1st iter of the 1st two-level step, so we mix it here const bool flag_restart = (this->two_level_step == 0 && PARAM.inp.init_wfc != "file") ? true : false; - if(this->exx_spacegroup_symmetry) - { this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm.get_dmk_vec(), pv), flag_restart); } - else - { this->mix_DMk_2D.mix(dm.get_dmk_vec(), flag_restart); } - const std::vector>>> Ds = - RI_2D_Comm::split_m2D_ktoR( - ucell, - *this->exx_ptr->p_kv, - this->mix_DMk_2D.get_DMk_out(), - pv, - nspin, - this->exx_spacegroup_symmetry); + const ModuleBase::Vector3* source_A = nullptr; + if (field != nullptr && field->gauge() == 2) + { + source_A = &field->A_right_ha(); + } + const auto Ds = this->prepare_exx_density(dm.get_dmk_vec(), ucell, kv, pv, + nspin, flag_restart, source_A); if(this->exx_spacegroup_symmetry && this->exx_ptr->info.exx_symmetry_realspace) { this->cal_exx_elec(Ds, ucell, pv, &this->symrot_); } else diff --git a/source/source_lcao/module_ri/hybrid_gauge_phase.h b/source/source_lcao/module_ri/hybrid_gauge_phase.h new file mode 100644 index 0000000000..a911488168 --- /dev/null +++ b/source/source_lcao/module_ri/hybrid_gauge_phase.h @@ -0,0 +1,195 @@ +#ifndef HYBRID_GAUGE_PHASE_H +#define HYBRID_GAUGE_PHASE_H + +#include "source_base/vector3.h" +#include "source_base/matrix3.h" +#include "source_lcao/module_ri/mix_dmk_2d.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace RI_2D_Comm +{ +// The RI Fourier convention is exp(-ik.R). Its sHG input therefore needs +// exp(-i A_source.R) to enter the ordinary NAO frame. A_source must identify +// the density's time layer, not merely the current Hamiltonian's stage. +inline std::complex hybrid_input_phase( + const ModuleBase::Vector3& vector_potential, + const ModuleBase::Vector3& lattice_translation) +{ + const double argument = -(vector_potential * lattice_translation); + return std::exp(std::complex(0.0, argument)); +} + +enum class HybridDensityTime { Left, Right }; + +// Restarting SCF mixing does not change the time layer of the current density. +inline HybridDensityTime hybrid_density_time(const int step, const int iteration) +{ + return step > 0 && iteration == 1 ? HybridDensityTime::Left : HybridDensityTime::Right; +} + +using HybridGaugeComplexDensity = std::vector< + std::map>, RI::Tensor>>>>; + +using HybridDensityBlock = std::tuple, std::vector>; +using HybridDensityLayout = std::vector>; + +inline HybridDensityLayout hybrid_density_layout(const HybridGaugeComplexDensity& density) +{ + HybridDensityLayout layout(density.size()); + for (std::size_t is = 0; is < density.size(); ++is) + { + for (const auto& row_atom : density[is]) + { + for (const auto& block : row_atom.second) + { + if (block.second.empty()) + { + throw std::invalid_argument("Empty RI density block in common-frame mixing"); + } + const std::vector shape(block.second.shape.begin(), block.second.shape.end()); + layout[is].emplace_back(row_atom.first, block.first.first, block.first.second, shape); + } + } + } + return layout; +} + +inline std::vector>> flatten_hybrid_density( + const HybridGaugeComplexDensity& density) +{ + std::vector>> flat(density.size()); + for (std::size_t is = 0; is < density.size(); ++is) + { + for (const auto& row_atom : density[is]) + { + for (const auto& block : row_atom.second) + { + const auto& tensor = block.second; + if (tensor.empty()) + { + throw std::invalid_argument("Empty RI density block in common-frame mixing"); + } + const std::size_t size = tensor.shape.get_shape_all(); + flat[is].insert(flat[is].end(), tensor.ptr(), tensor.ptr() + size); + } + } + } + return flat; +} + +inline void assign_flat_hybrid_density( + const std::vector>*>& flat, + const HybridDensityLayout& expected_layout, + HybridGaugeComplexDensity& density) +{ + const auto actual_layout = hybrid_density_layout(density); + if (flat.size() != density.size() || actual_layout != expected_layout) + { + throw std::invalid_argument("RI density topology changed during common-frame mixing"); + } + // Validate the whole transaction before writing any block. + for (std::size_t is = 0; is < density.size(); ++is) + { + std::size_t expected_size = 0; + for (const auto& row_atom : density[is]) + { + for (const auto& block : row_atom.second) + { + expected_size += block.second.shape.get_shape_all(); + } + } + if (flat[is] == nullptr || flat[is]->size() != expected_size) + { + throw std::invalid_argument("RI density size changed during common-frame mixing"); + } + } + for (std::size_t is = 0; is < density.size(); ++is) + { + std::size_t offset = 0; + for (auto& row_atom : density[is]) + { + for (auto& block : row_atom.second) + { + const std::size_t size = block.second.shape.get_shape_all(); + // Avoid writing through a shallow copy into accepted-state storage. + block.second = block.second.copy(); + std::copy(flat[is]->begin() + offset, flat[is]->begin() + offset + size, block.second.ptr()); + offset += size; + } + } + } +} + +// Tensor copies share storage in LibRI. Allocate a new phased tensor rather +// than changing entries in place, so cached source density remains untouched. +inline void apply_hybrid_input_phase(const ModuleBase::Matrix3& lattice, + const double lat0, + const ModuleBase::Vector3& source_A, + HybridGaugeComplexDensity& density) +{ + for (int idir = 0; idir < 3; ++idir) + { + if (!std::isfinite(source_A[idir])) + { + throw std::invalid_argument("Nonfinite RI density source phase"); + } + } + for (auto& spin_density : density) + { + for (auto& row_atom : spin_density) + { + for (auto& block : row_atom.second) + { + const auto& cell = block.first.second; + const ModuleBase::Vector3 cell_index(cell[0], cell[1], cell[2]); + const ModuleBase::Vector3 translation = (cell_index * lattice) * lat0; + const std::complex phase = hybrid_input_phase(source_A, translation); + block.second = block.second * phase; + } + } + } +} + +// The production transaction: transform first, validate topology before mixing, +// then copy the mixed values back without changing shared source storage. +inline void mix_hybrid_density(const ModuleBase::Matrix3& lattice, + const double lat0, + const ModuleBase::Vector3& source_A, + const bool restart, + Mix_DMk_2D>& mixer, + HybridDensityLayout& history_layout, + HybridGaugeComplexDensity& density) +{ + if (density.empty()) + { + throw std::invalid_argument("Missing spin blocks in RI density mixing"); + } + const auto layout = hybrid_density_layout(density); + if (!restart && layout != history_layout) + { + throw std::invalid_argument("RI density topology changed without a mixing restart"); + } + apply_hybrid_input_phase(lattice, lat0, source_A, density); + const auto flat = flatten_hybrid_density(density); + if (restart) + { + mixer.set_nks(static_cast(density.size())); + history_layout = layout; + } + mixer.mix(flat, restart); + assign_flat_hybrid_density(mixer.get_DMk_out(), history_layout, density); +} +} // namespace RI_2D_Comm + +#endif diff --git a/source/source_lcao/module_ri/unittests/CMakeLists.txt b/source/source_lcao/module_ri/unittests/CMakeLists.txt index 5635d5d68a..8db92f2449 100644 --- a/source/source_lcao/module_ri/unittests/CMakeLists.txt +++ b/source/source_lcao/module_ri/unittests/CMakeLists.txt @@ -28,3 +28,15 @@ AddTest( ) file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/support" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}") + +AddTest( + TARGET MODULE_RI_hybrid_gauge_phase_test + LIBS parameter base device + SOURCES test_hybrid_gauge_phase.cpp ../mix_dmk_2d.cpp +) + +AddTest( + TARGET MODULE_RI_exx_info_test + LIBS parameter base device + SOURCES test_exx_info.cpp ../../../source_hamilt/module_xc/exx_info.cpp +) diff --git a/source/source_lcao/module_ri/unittests/test_exx_info.cpp b/source/source_lcao/module_ri/unittests/test_exx_info.cpp new file mode 100644 index 0000000000..69ca823616 --- /dev/null +++ b/source/source_lcao/module_ri/unittests/test_exx_info.cpp @@ -0,0 +1,61 @@ +#include "source_hamilt/module_xc/exx_info.h" +#include "source_io/module_parameter/input_parameter.h" + +#include +#include + +namespace +{ +Input_para hse_input() +{ + Input_para input; + input.basis_type = "lcao"; + input.esolver_type = "tddft"; + input.dft_functional = "hse"; + input.exx_fock_alpha.clear(); + input.exx_fock_lambda.clear(); + input.exx_erfc_alpha = {"0.25"}; + input.exx_erfc_omega = {"0.11"}; + input.exx_ccp_rmesh_times = "1.5"; + input.exx_real_number = "1"; + input.td_stype = 2; + return input; +} +} + +TEST(ExxInfo, RejectsRealHybridGaugeBeforeRiSetup) +{ + const auto input = hse_input(); + Exx_Info info; + try + { + init_exx_info(info, input); + FAIL() << "Real RI must be rejected for hybrid-gauge EXX"; + } + catch (const std::invalid_argument& error) + { + EXPECT_STREQ(error.what(), "td_stype=2 EXX requires exx_real_number=0"); + } +} + +TEST(ExxInfo, AcceptsComplexHybridGauge) +{ + auto input = hse_input(); + input.exx_real_number = "0"; + Exx_Info info; + EXPECT_NO_THROW(init_exx_info(info, input)); + EXPECT_FALSE(info.info_ri.real_number); +} + +TEST(ExxInfo, PreservesRealRiForOtherGaugesAndGroundState) +{ + auto input = hse_input(); + input.td_stype = 1; + Exx_Info info; + EXPECT_NO_THROW(init_exx_info(info, input)); + EXPECT_TRUE(info.info_ri.real_number); + input.td_stype = 2; + input.esolver_type = "ksdft"; + EXPECT_NO_THROW(init_exx_info(info, input)); + EXPECT_TRUE(info.info_ri.real_number); +} diff --git a/source/source_lcao/module_ri/unittests/test_hybrid_gauge_phase.cpp b/source/source_lcao/module_ri/unittests/test_hybrid_gauge_phase.cpp new file mode 100644 index 0000000000..1bdb8c0759 --- /dev/null +++ b/source/source_lcao/module_ri/unittests/test_hybrid_gauge_phase.cpp @@ -0,0 +1,266 @@ +#include "source_lcao/module_ri/hybrid_gauge_phase.h" +#include "source_lcao/module_ri/mix_dmk_2d.h" +#include "source_base/module_mixing/pulay_mixing.h" +#include "source_base/module_mixing/broyden_mixing.h" + +#include "gtest/gtest.h" + +#include +#include + +namespace +{ +void check_borrowed_charge_mixer(Base_Mixing::Mixing& charge_mixer) +{ + const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0); + const auto key = std::make_pair(1, std::array{{1, 0, 0}}); + const std::complex physical(0.8, -0.1); + RI_2D_Comm::HybridGaugeComplexDensity density(1); + density[0][0][key] = RI::Tensor>({1, 1}); + Mix_DMk_2D> mixer; + mixer.set_mixing(&charge_mixer); + RI_2D_Comm::HybridDensityLayout layout; + Base_Mixing::Mixing_Data charge_history; + charge_mixer.init_mixing_data(charge_history, 2, sizeof(double)); + const auto inner_product = [](double* first, double* second) + { + return first[0] * second[0] + first[1] * second[1]; + }; + for (int iteration = 1; iteration <= 6; ++iteration) + { + // Also exercise a reset inside the same time step. Source time remains + // the endpoint even though both charge and EXX histories are restarted. + const bool restart = iteration == 1 || iteration == 4; + if (restart) + { + charge_mixer.reset(); + charge_mixer.init_mixing_data(charge_history, 2, sizeof(double)); + } + const auto time = RI_2D_Comm::hybrid_density_time(2, iteration); + const double left_value = 0.1; + const double right_value = 0.2 * iteration; + const double source_value = time == RI_2D_Comm::HybridDensityTime::Left ? left_value : right_value; + const ModuleBase::Vector3 source_A(source_value, 0.0, 0.0); + const double raw_source_value = iteration == 1 ? left_value : right_value; + density[0][0][key](0, 0) = physical * std::exp(std::complex(0.0, raw_source_value)); + RI_2D_Comm::mix_hybrid_density(lattice, 1.0, source_A, restart, mixer, layout, density); + EXPECT_NEAR(std::abs(density[0][0][key](0, 0) - physical), 0.0, 1e-12); + // Charge mixing supplies the coefficients consumed by EXX next iteration. + const double input[2] = {0.2, 0.4}; + const double output[2] = {0.3 + 0.1 * iteration, 0.8 - 0.01 * iteration * iteration}; + charge_mixer.push_data(charge_history, input, output, nullptr, true); + charge_mixer.cal_coef(charge_history, inner_product); + } +} +} // namespace + +TEST(HybridGaugePhase, BorrowsPulayCoefficientsAcrossRestart) +{ + Base_Mixing::Pulay_Mixing charge_mixer(2, 0.4); + check_borrowed_charge_mixer(charge_mixer); +} + +TEST(HybridGaugePhase, BorrowsBroydenCoefficientsAcrossRestart) +{ + Base_Mixing::Broyden_Mixing charge_mixer(2, 0.4); + check_borrowed_charge_mixer(charge_mixer); +} + +TEST(HybridGaugePhase, NegativeRIConventionAndRoundTrip) +{ + const ModuleBase::Vector3 source_A(0.3, -0.1, 0.2); + const ModuleBase::Vector3 translation(2.0, 3.0, -1.0); + const std::complex phase = RI_2D_Comm::hybrid_input_phase(source_A, translation); + EXPECT_NEAR(phase.real(), std::cos(0.1), 1e-15); + EXPECT_NEAR(phase.imag(), -std::sin(0.1), 1e-15); + const ModuleBase::Vector3 reverse_A(-0.3, 0.1, -0.2); + const std::complex reverse = RI_2D_Comm::hybrid_input_phase(reverse_A, translation); + EXPECT_NEAR(std::abs(phase * reverse - std::complex(1.0, 0.0)), 0.0, 1e-15); +} + +TEST(HybridGaugePhase, SelectsExplicitDensityTimeLayer) +{ + using RI_2D_Comm::HybridDensityTime; + EXPECT_EQ(RI_2D_Comm::hybrid_density_time(0, 1), HybridDensityTime::Right); + EXPECT_EQ(RI_2D_Comm::hybrid_density_time(1, 1), HybridDensityTime::Left); + EXPECT_EQ(RI_2D_Comm::hybrid_density_time(1, 2), HybridDensityTime::Right); + EXPECT_EQ(RI_2D_Comm::hybrid_density_time(2, 5), HybridDensityTime::Right); +} + +TEST(HybridGaugePhase, ExplicitSourceFrameAndUnchangedCache) +{ + const std::array cell = {{1, -1, 0}}; + const auto key = std::make_pair(2, cell); + RI_2D_Comm::HybridGaugeComplexDensity original(2); + original[1][0][key] = RI::Tensor>({1, 2}); + original[1][0][key](0, 0) = std::complex(0.7, 0.4); + original[1][0][key](0, 1) = std::complex(-0.2, 0.1); + auto transformed = original; + const ModuleBase::Matrix3 lattice(2.0, 0.5, 0.0, + 0.0, 3.0, 0.2, + 0.0, 0.0, 4.0); + const ModuleBase::Vector3 source_A(0.2, 0.1, -0.3); + RI_2D_Comm::apply_hybrid_input_phase(lattice, 2.0, source_A, transformed); + // R=(a1-a2)*lat0=(4,-5,-0.4), A.R=0.42. + const std::complex expected_phase(std::cos(0.42), -std::sin(0.42)); + const std::complex source(0.7, 0.4); + EXPECT_NEAR(std::abs(transformed[1][0][key](0, 0) - source * expected_phase), 0.0, 1e-15); + EXPECT_EQ(original[1][0][key](0, 0), source); + EXPECT_NE(original[1][0][key].ptr(), transformed[1][0][key].ptr()); + EXPECT_TRUE(transformed[0].empty()); +} + +TEST(HybridGaugePhase, EqualPhysicalDensityAcrossTwoSourceFrames) +{ + const ModuleBase::Vector3 translation(1.1, -0.5, 0.3); + const ModuleBase::Vector3 left_A(0.1, 0.2, -0.1); + const ModuleBase::Vector3 right_A(0.4, -0.1, 0.2); + const std::complex physical(0.7, -0.2); + const auto left_phase = RI_2D_Comm::hybrid_input_phase(left_A, translation); + const auto right_phase = RI_2D_Comm::hybrid_input_phase(right_A, translation); + const auto left_raw = physical / left_phase; + const auto right_raw = physical / right_phase; + const auto common_mixed = 0.4 * left_raw * left_phase + 0.6 * right_raw * right_phase; + EXPECT_NEAR(std::abs(common_mixed - physical), 0.0, 1e-15); + const auto wrong_stage_mixed = (0.4 * left_raw + 0.6 * right_raw) * right_phase; + EXPECT_GT(std::abs(wrong_stage_mixed - physical), 0.01); +} + +TEST(HybridGaugePhase, ImageRelabelCommutesWithOrdinaryFrameConversion) +{ + // Algebraic oracle for the inspected conj-first DMK / RI transpose + // conventions. It does not enable transport of production snapshots. + const ModuleBase::Matrix3 lattice(2.0, 0.5, 0.0, + 0.0, 3.0, 0.2, + 0.1, 0.0, 4.0); + const double lat0 = 1.7; + const ModuleBase::Vector3 a(0.2, -0.1, 0.07); + const ModuleBase::Vector3 k_phase(0.31, 0.19, -0.23); + const ModuleBase::Vector3 bi(-1, 0, 1); + const ModuleBase::Vector3 bj(0, 1, -1); + const ModuleBase::Vector3 r(1, -2, 0); + const auto delta = bi - bj; + const auto rn = r + delta; + const auto cartesian = [&](const ModuleBase::Vector3& image) + { + return ModuleBase::Vector3(image.x, image.y, image.z) * lattice * lat0; + }; + const auto coefficient_phase = [&](const ModuleBase::Vector3& image) + { + const double argument = k_phase.x * image.x + k_phase.y * image.y + + k_phase.z * image.z + a * cartesian(image); + return std::exp(std::complex(0.0, argument)); + }; + const std::complex ci(0.7, -0.2); + const std::complex cj(-0.1, 0.6); + const std::complex ui = coefficient_phase(bi); + const std::complex uj = coefficient_phase(bj); + // dm_from_psi stores conj(C_i)*C_j; RI's row-major tensor view + // transposes its column-major buffer before Fourier conversion. + const auto stored_old = std::conj(ci) * cj; + const auto stored_new = std::conj(ui * ci) * (uj * cj); + EXPECT_NEAR(std::abs(stored_new - std::conj(ui) * stored_old * uj), 0.0, 1e-14); + EXPECT_GT(std::abs(stored_new - ui * stored_old * std::conj(uj)), 1e-3); + + const double k_delta = k_phase.x * delta.x + k_phase.y * delta.y + k_phase.z * delta.z; + const auto transposed_new_ratio = ui * std::conj(uj); + const auto fourier_relabel = std::exp(std::complex(0.0, -k_delta)); + const auto raw_ri_ratio = transposed_new_ratio * fourier_relabel; + EXPECT_NEAR(std::abs(raw_ri_ratio - std::exp(std::complex(0.0, a * cartesian(delta)))), + 0.0, 1e-14); + + const auto old_key = std::make_pair(1, std::array{{r.x, r.y, r.z}}); + const auto new_key = std::make_pair(1, std::array{{rn.x, rn.y, rn.z}}); + RI_2D_Comm::HybridGaugeComplexDensity old_density(1); + RI_2D_Comm::HybridGaugeComplexDensity new_density(1); + old_density[0][0][old_key] = RI::Tensor>({1, 1}); + new_density[0][0][new_key] = RI::Tensor>({1, 1}); + const std::complex raw(0.8, -0.3); + old_density[0][0][old_key](0, 0) = raw; + new_density[0][0][new_key](0, 0) = raw * raw_ri_ratio; + const auto old_source = old_density; + const auto new_source = new_density; + RI_2D_Comm::apply_hybrid_input_phase(lattice, lat0, a, old_density); + RI_2D_Comm::apply_hybrid_input_phase(lattice, lat0, a, new_density); + EXPECT_NEAR(std::abs(old_density[0][0][old_key](0, 0) - new_density[0][0][new_key](0, 0)), + 0.0, 1e-14); + EXPECT_EQ(old_source[0].at(0).at(old_key)(0, 0), raw); + EXPECT_EQ(new_source[0].at(0).at(new_key)(0, 0), raw * raw_ri_ratio); +} + +TEST(HybridGaugePhase, ProductionPlainMixerAndTransactionalTopologyGuard) +{ + const std::array cell = {{1, 0, 0}}; + const auto key = std::make_pair(1, cell); + RI_2D_Comm::HybridGaugeComplexDensity density(1); + density[0][0][key] = RI::Tensor>({1, 2}); + density[0][0][key](0, 0) = std::complex(0.5, 0.2); + density[0][0][key](0, 1) = std::complex(0.3, -0.1); + const auto accepted = density; + const auto layout = RI_2D_Comm::hybrid_density_layout(density); + auto first = RI_2D_Comm::flatten_hybrid_density(density); + auto second = first; + second[0][0] += std::complex(0.4, -0.2); + Mix_DMk_2D> mixer; + mixer.set_nks(1); + mixer.set_mixing_plain(0.25); + mixer.mix(first, true); + mixer.mix(second, false); + const auto output = mixer.get_DMk_out(); + RI_2D_Comm::assign_flat_hybrid_density(output, layout, density); + EXPECT_NEAR(std::abs(density[0][0][key](0, 0) - std::complex(0.6, 0.15)), 0.0, 1e-15); + EXPECT_EQ(accepted[0].at(0).at(key)(0, 0), first[0][0]); + + auto changed = density; + const std::array other_cell = {{2, 0, 0}}; + const auto other_key = std::make_pair(1, other_cell); + changed[0][0][other_key] = changed[0][0][key]; + changed[0][0].erase(key); + EXPECT_THROW(RI_2D_Comm::assign_flat_hybrid_density(output, layout, changed), std::invalid_argument); + EXPECT_EQ(changed[0][0][other_key](0, 0), density[0][0][key](0, 0)); + changed = density; + changed[0][0][key] = RI::Tensor>({2, 1}); + EXPECT_THROW(RI_2D_Comm::assign_flat_hybrid_density(output, layout, changed), std::invalid_argument); + const std::vector>*> invalid(1, nullptr); + EXPECT_THROW(RI_2D_Comm::assign_flat_hybrid_density(invalid, layout, density), std::invalid_argument); +} + +TEST(HybridGaugePhase, ProductionMixerReceivesCommonFrameNotRawSourceFrames) +{ + const ModuleBase::Matrix3 lattice(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0); + const ModuleBase::Vector3 left_A(0.2, 0.0, 0.0); + const ModuleBase::Vector3 right_A(0.5, 0.0, 0.0); + const std::array cell = {{1, 0, 0}}; + const auto key = std::make_pair(1, cell); + const std::complex physical(0.8, -0.1); + RI_2D_Comm::HybridGaugeComplexDensity left(1); + RI_2D_Comm::HybridGaugeComplexDensity right(1); + left[0][0][key] = RI::Tensor>({1, 1}); + right[0][0][key] = RI::Tensor>({1, 1}); + left[0][0][key](0, 0) = physical * std::exp(std::complex(0.0, 0.2)); + right[0][0][key](0, 0) = physical * std::exp(std::complex(0.0, 0.5)); + Mix_DMk_2D> mixer; + mixer.set_mixing_plain(0.6); + RI_2D_Comm::HybridDensityLayout layout; + RI_2D_Comm::mix_hybrid_density(lattice, 1.0, left_A, true, mixer, layout, left); + RI_2D_Comm::mix_hybrid_density(lattice, 1.0, right_A, false, mixer, layout, right); + EXPECT_NEAR(std::abs(right[0][0][key](0, 0) - physical), 0.0, 1e-15); + + const auto before = *mixer.get_DMk_out()[0]; + auto changed = right; + const auto other_key = std::make_pair(1, std::array{{2, 0, 0}}); + changed[0][0][other_key] = changed[0][0][key]; + changed[0][0].erase(key); + EXPECT_THROW(RI_2D_Comm::mix_hybrid_density(lattice, 1.0, right_A, false, + mixer, layout, changed), std::invalid_argument); + EXPECT_EQ(*mixer.get_DMk_out()[0], before); + EXPECT_EQ(changed[0][0][other_key](0, 0), right[0][0][key](0, 0)); + const ModuleBase::Vector3 invalid_A(std::numeric_limits::infinity(), 0.0, 0.0); + EXPECT_THROW(RI_2D_Comm::mix_hybrid_density(lattice, 1.0, invalid_A, true, + mixer, layout, changed), std::invalid_argument); + EXPECT_EQ(*mixer.get_DMk_out()[0], before); + RI_2D_Comm::mix_hybrid_density(lattice, 1.0, right_A, true, mixer, layout, changed); + EXPECT_EQ(layout, RI_2D_Comm::hybrid_density_layout(changed)); +}