diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 31b025652b..acbe4e2bb8 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -343,18 +343,28 @@ OBJS_HAMILT_OF=kedf_tf.o\ OBJS_HAMILT_LCAO=hamilt_lcao.o\ operator_lcao.o\ ekinetic.o\ + ekinetic_fs.o\ + ekinetic_dh.o\ nonlocal.o\ + nonlocal_dh.o\ + nonlocal_fs.o\ overlap.o\ + overlap_fs.o\ td_ekinetic_lcao.o\ td_nonlocal_lcao.o\ td_pot_hybrid.o\ + td_pot_hybrid_fs.o\ veff_lcao.o\ + veff_dh.o\ meta_lcao.o\ op_dftu_lcao.o\ deepks_lcao.o\ op_exx_lcao.o\ dspin_lcao.o\ + dspin_fs.o\ dftu_lcao.o\ + dftu_fs.o\ + operator_fs_utils.o\ OBJS_HCONTAINER=base_matrix.o\ atom_pair.o\ @@ -721,7 +731,7 @@ OBJS_SRCPW=H_Ewald_pw.o\ setup_pwrho.o\ setup_pwwfc.o\ update_cell_pw.o\ - dftu_pw.o\ + setup_dftu_pw.o\ deltaspin_pw.o\ forces.o\ forces_us.o\ @@ -793,7 +803,7 @@ OBJS_DFTU=dftu.o\ dftu_tools.o\ dftu_occup.o\ dftu_hamilt.o\ - dftu_pw.o + setup_dftu_pw.o OBJS_DELTASPIN=basic_funcs.o\ cal_mw_from_lambda.o\ diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 666ea1278d..ef08d799bf 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -5,7 +5,7 @@ #include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_lcao/module_deltaspin/deltaspin_lcao.h" #include "source_lcao/dftu_lcao.h" -#include "source_lcao/hs_matrix_k.hpp" // there may be multiple definitions if using hpp +#include "source_hamilt/hs_matrix_k.h" #include "source_estate/module_charge/symmetry_rho.h" #include "source_lcao/LCAO_domain.h" // need DeePKS_init #include "source_lcao/FORCE_STRESS.h" diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 387137e349..2a7920ade0 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -1,21 +1,14 @@ #include "esolver_ks_pw.h" -#include "source_cell/cal_ux.h" #include "source_estate/elecstate_pw.h" #include "source_estate/module_charge/symmetry_rho.h" -#include "source_hamilt/module_xc/xc_functional.h" // use XC_Functional #include "source_hsolver/diago_iter_assist.h" #include "source_hsolver/diago_params.h" #include "source_hsolver/hsolver_pw.h" -#include "source_hsolver/kernels/hegvd_op.h" #include "source_io/module_parameter/parameter.h" -#include "source_lcao/module_deltaspin/spin_constrain.h" -#include "source_lcao/module_dftu/dftu.h" #include "source_pw/module_pwdft/forces.h" #include "source_pw/module_pwdft/hamilt_pw.h" -#include "source_pw/module_pwdft/onsite_proj.h" #include "source_pw/module_pwdft/stress_pw.h" -#include "source_pw/module_pwdft/vsep_pw.h" #ifdef __DSP #include "source_base/kernels/dsp/dsp_connector.h" @@ -23,13 +16,12 @@ #include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006 #include "source_estate/setup_estate_pw.h" // mohan add 20251005 -#include "source_estate/update_pot.h" // mohan add 20251016 #include "source_hamilt/module_xc/exx_info.h" // use GlobalC::exx_info #include "source_io/module_ctrl/ctrl_output_pw.h" // mohan add 20250927 #include "source_pw/module_pwdft/deltaspin_pw.h" // mohan add 20250309 -#include "source_pw/module_pwdft/dftu_pw.h" // mohan add 20250309 #include "source_pw/module_pwdft/setup_pot.h" // mohan add 20250929 #include "source_pw/module_pwdft/update_cell_pw.h" // mohan add 20250309 +#include "source_pw/module_pwdft/setup_dftu_pw.h" // mohan add 20250309 namespace ModuleESolver { @@ -150,7 +142,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) if (ucell.cell_parameter_updated) { - this->stp.p_psi_init->prepare_init(PARAM.inp.pw_seed); + this->stp.p_psi_init->prepare_init(PARAM.inp.pw_seed, istep); } //! Init Hamiltonian (cell changed) diff --git a/source/source_hamilt/CMakeLists.txt b/source/source_hamilt/CMakeLists.txt index 278389af0c..eaad08cf3f 100644 --- a/source/source_hamilt/CMakeLists.txt +++ b/source/source_hamilt/CMakeLists.txt @@ -9,6 +9,7 @@ endif() list(APPEND objects operator.cpp + hs_matrix_k.cpp module_ewald/H_Ewald_pw.cpp module_ewald/dnrm2.cpp ) diff --git a/source/source_hamilt/hs_matrix_k.cpp b/source/source_hamilt/hs_matrix_k.cpp new file mode 100644 index 0000000000..348701b547 --- /dev/null +++ b/source/source_hamilt/hs_matrix_k.cpp @@ -0,0 +1,7 @@ +#include "source_hamilt/hs_matrix_k.h" + +namespace hamilt +{ + template class HS_Matrix_K; + template class HS_Matrix_K>; +} \ No newline at end of file diff --git a/source/source_lcao/hs_matrix_k.hpp b/source/source_hamilt/hs_matrix_k.h similarity index 92% rename from source/source_lcao/hs_matrix_k.hpp rename to source/source_hamilt/hs_matrix_k.h index a347bd2e31..2df0ecf1f7 100644 --- a/source/source_lcao/hs_matrix_k.hpp +++ b/source/source_hamilt/hs_matrix_k.h @@ -1,8 +1,9 @@ -#ifndef HS_MATRIX_K_HPP -#define HS_MATRIX_K_HPP +#ifndef HS_MATRIX_K_H +#define HS_MATRIX_K_H #include "source_basis/module_ao/parallel_orbitals.h" +#include #include namespace hamilt { @@ -42,4 +43,4 @@ namespace hamilt }; } -#endif +#endif \ No newline at end of file diff --git a/source/source_io/module_dhs/write_dH_terms.cpp b/source/source_io/module_dhs/write_dH_terms.cpp index 9f725c6102..43130abc4f 100644 --- a/source/source_io/module_dhs/write_dH_terms.cpp +++ b/source/source_io/module_dhs/write_dH_terms.cpp @@ -6,7 +6,7 @@ #include "source_hamilt/module_hcontainer/output_hcontainer.h" #include "source_lcao/module_operator_lcao/ekinetic.h" #include "source_lcao/module_operator_lcao/nonlocal.h" -#include "source_lcao/module_operator_lcao/operator_force_stress_utils.h" +#include "source_lcao/module_operator_lcao/operator_fs_utils.h" #include "source_lcao/module_operator_lcao/veff_lcao.h" #include "source_hamilt/module_gint/gint_interface.h" #include "source_lcao/module_lr/utils/lr_util_xc.hpp" diff --git a/source/source_io/module_energy/write_eband_terms.hpp b/source/source_io/module_energy/write_eband_terms.hpp index 45a51bf125..327fdfea66 100644 --- a/source/source_io/module_energy/write_eband_terms.hpp +++ b/source/source_io/module_energy/write_eband_terms.hpp @@ -29,8 +29,8 @@ void write_eband_terms(const int nspin, const TwoCenterBundle& two_center_bundle #ifdef __EXX , - std::vector>>>* Hexxd = nullptr, - std::vector>>>>* Hexxc = nullptr + std::vector>>>* Hexxd = nullptr, + std::vector>>>>* Hexxc = nullptr #endif ) { diff --git a/source/source_io/module_hs/write_H_terms.cpp b/source/source_io/module_hs/write_H_terms.cpp index 481f49cbed..20115c1d2d 100644 --- a/source/source_io/module_hs/write_H_terms.cpp +++ b/source/source_io/module_hs/write_H_terms.cpp @@ -14,7 +14,7 @@ #include "source_hamilt/module_hcontainer/output_hcontainer.h" #include "source_lcao/module_operator_lcao/ekinetic.h" #include "source_lcao/module_operator_lcao/nonlocal.h" -#include "source_lcao/module_operator_lcao/operator_force_stress_utils.h" +#include "source_lcao/module_operator_lcao/operator_fs_utils.h" #ifdef __EXX #include "source_lcao/module_operator_lcao/op_exx_lcao.h" #include "source_lcao/module_ri/Exx_LRI_interface.h" diff --git a/source/source_io/module_hs/write_vxc.hpp b/source/source_io/module_hs/write_vxc.hpp index d6757f3cf0..7f8d2c8134 100644 --- a/source/source_io/module_hs/write_vxc.hpp +++ b/source/source_io/module_hs/write_vxc.hpp @@ -156,8 +156,8 @@ void write_Vxc(const int nspin, bool cal_exx #ifdef __EXX , - std::vector>>>* Hexxd = nullptr, - std::vector>>>>* Hexxc = nullptr + std::vector>>>* Hexxd = nullptr, + std::vector>>>>* Hexxc = nullptr #endif ) { diff --git a/source/source_io/module_hs/write_vxc_r.hpp b/source/source_io/module_hs/write_vxc_r.hpp index 81d4b38a10..4aa3ad0e52 100644 --- a/source/source_io/module_hs/write_vxc_r.hpp +++ b/source/source_io/module_hs/write_vxc_r.hpp @@ -41,8 +41,8 @@ void write_Vxc_R(const int nspin, bool real_number #ifdef __EXX , - const std::vector>>>* const Hexxd, - const std::vector>>>>* const Hexxc + const std::vector>>>* const Hexxd, + const std::vector>>>>* const Hexxc #endif , const double sparse_thr = 1e-10) diff --git a/source/source_lcao/CMakeLists.txt b/source/source_lcao/CMakeLists.txt index 1658c95737..b9fd77a947 100644 --- a/source/source_lcao/CMakeLists.txt +++ b/source/source_lcao/CMakeLists.txt @@ -9,19 +9,28 @@ if(ENABLE_LCAO) hamilt_lcao.cpp module_operator_lcao/operator_lcao.cpp module_operator_lcao/veff_lcao.cpp + module_operator_lcao/veff_dh.cpp module_operator_lcao/meta_lcao.cpp module_operator_lcao/op_dftu_lcao.cpp module_operator_lcao/deepks_lcao.cpp module_operator_lcao/op_exx_lcao.cpp module_operator_lcao/overlap.cpp + module_operator_lcao/overlap_fs.cpp module_operator_lcao/ekinetic.cpp + module_operator_lcao/ekinetic_fs.cpp + module_operator_lcao/ekinetic_dh.cpp module_operator_lcao/nonlocal.cpp + module_operator_lcao/nonlocal_dh.cpp + module_operator_lcao/nonlocal_fs.cpp module_operator_lcao/td_ekinetic_lcao.cpp module_operator_lcao/td_nonlocal_lcao.cpp module_operator_lcao/td_pot_hybrid.cpp + module_operator_lcao/td_pot_hybrid_fs.cpp module_operator_lcao/dspin_lcao.cpp + module_operator_lcao/dspin_fs.cpp module_operator_lcao/dftu_lcao.cpp - module_operator_lcao/operator_force_stress_utils.cpp + module_operator_lcao/dftu_fs.cpp + module_operator_lcao/operator_fs_utils.cpp dftu_lcao.cpp pulay_fs_center2.cpp FORCE_STRESS.cpp diff --git a/source/source_lcao/hamilt_lcao.h b/source/source_lcao/hamilt_lcao.h index 71392d069d..0c7a35f2af 100644 --- a/source/source_lcao/hamilt_lcao.h +++ b/source/source_lcao/hamilt_lcao.h @@ -5,7 +5,7 @@ #include "source_cell/klist.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_hamilt/hamilt.h" -#include "source_lcao/hs_matrix_k.hpp" +#include "source_hamilt/hs_matrix_k.h" #include "source_hamilt/module_hcontainer/hcontainer.h" #include diff --git a/source/source_lcao/module_deepks/test/deepks_test_e_deltabands.cpp b/source/source_lcao/module_deepks/test/deepks_test_e_deltabands.cpp index bf6906b91b..92c58867ac 100644 --- a/source/source_lcao/module_deepks/test/deepks_test_e_deltabands.cpp +++ b/source/source_lcao/module_deepks/test/deepks_test_e_deltabands.cpp @@ -1,6 +1,6 @@ #include "deepks_test_runner.h" -#include "source_lcao/hs_matrix_k.hpp" +#include "source_hamilt/hs_matrix_k.h" #include "source_lcao/module_operator_lcao/deepks_lcao.h" #include diff --git a/source/source_lcao/module_operator_lcao/CMakeLists.txt b/source/source_lcao/module_operator_lcao/CMakeLists.txt index e9dbbe381b..7e09e099a0 100644 --- a/source/source_lcao/module_operator_lcao/CMakeLists.txt +++ b/source/source_lcao/module_operator_lcao/CMakeLists.txt @@ -5,16 +5,25 @@ add_library( op_dftu_lcao.cpp meta_lcao.cpp veff_lcao.cpp + veff_dh.cpp deepks_lcao.cpp overlap.cpp + overlap_fs.cpp ekinetic.cpp + ekinetic_fs.cpp + ekinetic_dh.cpp nonlocal.cpp + nonlocal_dh.cpp + nonlocal_fs.cpp td_ekinetic_lcao.cpp td_nonlocal_lcao.cpp td_pot_hybrid.cpp + td_pot_hybrid_fs.cpp dspin_lcao.cpp + dspin_fs.cpp dftu_lcao.cpp - operator_force_stress_utils.cpp + dftu_fs.cpp + operator_fs_utils.cpp ) if(ENABLE_COVERAGE) diff --git a/source/source_lcao/module_operator_lcao/dftu_force_stress.hpp b/source/source_lcao/module_operator_lcao/dftu_fs.cpp similarity index 84% rename from source/source_lcao/module_operator_lcao/dftu_force_stress.hpp rename to source/source_lcao/module_operator_lcao/dftu_fs.cpp index e9a546c2b9..54014250bc 100644 --- a/source/source_lcao/module_operator_lcao/dftu_force_stress.hpp +++ b/source/source_lcao/module_operator_lcao/dftu_fs.cpp @@ -1,4 +1,3 @@ -#pragma once #include "dftu_lcao.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" @@ -426,4 +425,77 @@ void DFTU>::cal_stress_IJR(const int& iat1, } } +// explicit member function instantiations +template void DFTU>::cal_force_stress( + const bool cal_force, const bool cal_stress, + ModuleBase::matrix& force, ModuleBase::matrix& stress); +template void DFTU, double>>::cal_force_stress( + const bool cal_force, const bool cal_stress, + ModuleBase::matrix& force, ModuleBase::matrix& stress); +template void DFTU, std::complex>>::cal_force_stress( + const bool cal_force, const bool cal_stress, + ModuleBase::matrix& force, ModuleBase::matrix& stress); + +template void DFTU>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& vu_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + double* force1, double* force2); +template void DFTU, double>>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& vu_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + double* force1, double* force2); +template void DFTU, std::complex>>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& vu_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + double* force1, double* force2); + +template void DFTU>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& vu_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); +template void DFTU, double>>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& vu_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); +template void DFTU, std::complex>>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const std::vector& vu_in, + const hamilt::BaseMatrix** dmR_pointer, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); + } // namespace hamilt diff --git a/source/source_lcao/module_operator_lcao/dftu_lcao.cpp b/source/source_lcao/module_operator_lcao/dftu_lcao.cpp index f7b142e2f1..590d6596dd 100644 --- a/source/source_lcao/module_operator_lcao/dftu_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/dftu_lcao.cpp @@ -687,8 +687,6 @@ void hamilt::DFTU>::cal_v_of_u(const std::vector>; template class hamilt::DFTU, double>>; template class hamilt::DFTU, std::complex>>; diff --git a/source/source_lcao/module_operator_lcao/dspin_force_stress.hpp b/source/source_lcao/module_operator_lcao/dspin_fs.cpp similarity index 82% rename from source/source_lcao/module_operator_lcao/dspin_force_stress.hpp rename to source/source_lcao/module_operator_lcao/dspin_fs.cpp index 1fe8812e9b..cb0cbd96db 100644 --- a/source/source_lcao/module_operator_lcao/dspin_force_stress.hpp +++ b/source/source_lcao/module_operator_lcao/dspin_fs.cpp @@ -1,7 +1,8 @@ -#pragma once #include "dspin_lcao.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" +#include "source_lcao/module_deltaspin/spin_constrain.h" +#include "source_io/module_parameter/parameter.h" namespace hamilt { @@ -393,4 +394,80 @@ void DeltaSpin>::cal_stress_IJR(const int& iat1, } } +// explicit member function instantiations +template void DeltaSpin>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); +template void DeltaSpin, double>>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); +template void DeltaSpin, std::complex>>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); + +template void DeltaSpin>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& lambda, + const int nspin, + double* force1, double* force2); +template void DeltaSpin, double>>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& lambda, + const int nspin, + double* force1, double* force2); +template void DeltaSpin, std::complex>>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& lambda, + const int nspin, + double* force1, double* force2); + +template void DeltaSpin>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& lambda, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); +template void DeltaSpin, double>>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& lambda, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); +template void DeltaSpin, std::complex>>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& lambda, + const int nspin, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); + } // namespace hamilt diff --git a/source/source_lcao/module_operator_lcao/dspin_lcao.cpp b/source/source_lcao/module_operator_lcao/dspin_lcao.cpp index d1377a6f0c..2d994ac37f 100644 --- a/source/source_lcao/module_operator_lcao/dspin_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/dspin_lcao.cpp @@ -628,8 +628,6 @@ void hamilt::DeltaSpin>::cal_PI_sub( } } -#include "dspin_force_stress.hpp" - template class hamilt::DeltaSpin>; template class hamilt::DeltaSpin, double>>; template class hamilt::DeltaSpin, std::complex>>; \ No newline at end of file diff --git a/source/source_lcao/module_operator_lcao/dspin_lcao.h b/source/source_lcao/module_operator_lcao/dspin_lcao.h index a07d3676da..c0916d5b80 100644 --- a/source/source_lcao/module_operator_lcao/dspin_lcao.h +++ b/source/source_lcao/module_operator_lcao/dspin_lcao.h @@ -75,7 +75,7 @@ class DeltaSpin> : public OperatorLCAO * spin switch. sc_hr_done must be reset here so each spin's HR is computed * independently. */ - void set_current_spin(const int current_spin_in) + void set_current_spin(const int current_spin_in) override { if (this->current_spin != current_spin_in) { diff --git a/source/source_lcao/module_operator_lcao/ekinetic.cpp b/source/source_lcao/module_operator_lcao/ekinetic.cpp index be7425225c..cd7ae91083 100644 --- a/source/source_lcao/module_operator_lcao/ekinetic.cpp +++ b/source/source_lcao/module_operator_lcao/ekinetic.cpp @@ -257,10 +257,6 @@ void hamilt::EKinetic>::contributeHR() return; } -// Include force/stress implementation -#include "ekinetic_force_stress.hpp" -#include "ekinetic_dh.hpp" - template class hamilt::EKinetic>; template class hamilt::EKinetic, double>>; template class hamilt::EKinetic, std::complex>>; diff --git a/source/source_lcao/module_operator_lcao/ekinetic_dh.hpp b/source/source_lcao/module_operator_lcao/ekinetic_dh.cpp similarity index 93% rename from source/source_lcao/module_operator_lcao/ekinetic_dh.hpp rename to source/source_lcao/module_operator_lcao/ekinetic_dh.cpp index f06643e485..6f1f66fc78 100644 --- a/source/source_lcao/module_operator_lcao/ekinetic_dh.hpp +++ b/source/source_lcao/module_operator_lcao/ekinetic_dh.cpp @@ -1,6 +1,4 @@ -#pragma once #include "ekinetic.h" -#include "operator_force_stress_utils.hpp" #include "source_base/timer.h" namespace hamilt @@ -162,4 +160,12 @@ void EKinetic>::cal_dH(std::array>::cal_dH( + std::array*>, 3>& dhR); +template void EKinetic, double>>::cal_dH( + std::array*>, 3>& dhR); +template void EKinetic, std::complex>>::cal_dH( + std::array*>, 3>& dhR); + } // namespace hamilt diff --git a/source/source_lcao/module_operator_lcao/ekinetic_force_stress.hpp b/source/source_lcao/module_operator_lcao/ekinetic_force_stress.hpp deleted file mode 100644 index 5204c2c952..0000000000 --- a/source/source_lcao/module_operator_lcao/ekinetic_force_stress.hpp +++ /dev/null @@ -1,67 +0,0 @@ -#pragma once -#include "ekinetic.h" -#include "operator_force_stress_utils.hpp" -#include "source_base/timer.h" - -namespace hamilt -{ - -template -void EKinetic>::cal_force_stress(const bool cal_force, - const bool cal_stress, - const HContainer* dmR, - ModuleBase::matrix& force, - ModuleBase::matrix& stress) -{ - ModuleBase::TITLE("EKinetic", "cal_force_stress"); - ModuleBase::timer::start("EKinetic", "cal_force_stress"); - - // Lambda function to calculate kinetic integral and its gradient - auto integral_calc = [this](int T1, int L1, int N1, int M1, - int T2, int L2, int N2, int M2, - const ModuleBase::Vector3& dtau, - double* olm) { - this->intor_->calculate(T1, L1, N1, M1, T2, L2, N2, M2, - dtau * this->ucell->lat0, &olm[0], &olm[1]); - }; - - // Use unified template with ForceSign=+1, StressSign=-1 for kinetic operator - OperatorForceStress::cal_force_stress_2center( - cal_force, cal_stress, dmR, this->ucell, this->gridD, - this->orb_cutoff_, dmR->get_paraV(), integral_calc, force, stress); - - ModuleBase::timer::end("EKinetic", "cal_force_stress"); -} - -// Dummy implementations for cal_force_IJR and cal_stress_IJR -// These are not used in the simplified approach above -template -void EKinetic>::cal_force_IJR( - const int& iat1, - const int& iat2, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const hamilt::BaseMatrix* dmR_pointer, - double* force1, - double* force2) -{ - // Not used in current implementation -} - -template -void EKinetic>::cal_stress_IJR( - const int& iat1, - const int& iat2, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const hamilt::BaseMatrix* dmR_pointer, - const ModuleBase::Vector3& dis1, - const ModuleBase::Vector3& dis2, - double* stress) -{ - // Not used in current implementation -} - -} // namespace hamilt diff --git a/source/source_lcao/module_operator_lcao/ekinetic_fs.cpp b/source/source_lcao/module_operator_lcao/ekinetic_fs.cpp new file mode 100644 index 0000000000..dabaaa6791 --- /dev/null +++ b/source/source_lcao/module_operator_lcao/ekinetic_fs.cpp @@ -0,0 +1,130 @@ +#include "ekinetic.h" +#include "operator_fs_utils.hpp" +#include "source_base/timer.h" + +namespace hamilt +{ + +template +void EKinetic>::cal_force_stress(const bool cal_force, + const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, + ModuleBase::matrix& stress) +{ + ModuleBase::TITLE("EKinetic", "cal_force_stress"); + ModuleBase::timer::start("EKinetic", "cal_force_stress"); + + // Lambda function to calculate kinetic integral and its gradient + auto integral_calc = [this](int T1, int L1, int N1, int M1, + int T2, int L2, int N2, int M2, + const ModuleBase::Vector3& dtau, + double* olm) { + this->intor_->calculate(T1, L1, N1, M1, T2, L2, N2, M2, + dtau * this->ucell->lat0, &olm[0], &olm[1]); + }; + + // Use unified template with ForceSign=+1, StressSign=-1 for kinetic operator + OperatorForceStress::cal_force_stress_2center( + cal_force, cal_stress, dmR, this->ucell, this->gridD, + this->orb_cutoff_, dmR->get_paraV(), integral_calc, force, stress); + + ModuleBase::timer::end("EKinetic", "cal_force_stress"); +} + +// Dummy implementations for cal_force_IJR and cal_stress_IJR +// These are not used in the simplified approach above +template +void EKinetic>::cal_force_IJR( + const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + double* force1, + double* force2) +{ + // Not used in current implementation +} + +template +void EKinetic>::cal_stress_IJR( + const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress) +{ + // Not used in current implementation +} + +// explicit member function instantiations +template void EKinetic>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); +template void EKinetic, double>>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); +template void EKinetic, std::complex>>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); + +template void EKinetic>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + double* force1, double* force2); +template void EKinetic, double>>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + double* force1, double* force2); +template void EKinetic, std::complex>>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix>* dmR_pointer, + double* force1, double* force2); + +template void EKinetic>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); +template void EKinetic, double>>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); +template void EKinetic, std::complex>>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix>* dmR_pointer, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); + +} // namespace hamilt diff --git a/source/source_lcao/module_operator_lcao/nonlocal.cpp b/source/source_lcao/module_operator_lcao/nonlocal.cpp index c90f5c17e2..b73cc022d0 100644 --- a/source/source_lcao/module_operator_lcao/nonlocal.cpp +++ b/source/source_lcao/module_operator_lcao/nonlocal.cpp @@ -1,5 +1,5 @@ #include "nonlocal.h" - +#include "operator_fs_utils.h" #include "source_base/timer.h" #include "source_base/tool_title.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" @@ -323,9 +323,6 @@ void hamilt::Nonlocal>::contributeHR() return; } -#include "nonlocal_force_stress.hpp" -#include "nonlocal_dh.hpp" - template class hamilt::Nonlocal>; template class hamilt::Nonlocal, double>>; template class hamilt::Nonlocal, std::complex>>; diff --git a/source/source_lcao/module_operator_lcao/nonlocal.h b/source/source_lcao/module_operator_lcao/nonlocal.h index 3c76fe397c..4dae60781c 100644 --- a/source/source_lcao/module_operator_lcao/nonlocal.h +++ b/source/source_lcao/module_operator_lcao/nonlocal.h @@ -61,7 +61,71 @@ class Nonlocal> : public OperatorLCAO ModuleBase::matrix& force, ModuleBase::matrix& stress); - // per-atom-I derivative d/dtau_I; one HContainer per atom I (size nat each) + /** + * @brief Calculate the derivative of the non-local pseudopotential Hamiltonian + * with respect to atomic displacements: dH^{NL}/dτ_I + * + * For each atom I and each Cartesian direction d ∈ {0,1,2} (x,y,z), + * fills an HContainer matrix dhR[d][I] that stores + * dH^{NL}_{IJ}/dτ_I for all atom-pair blocks (I,J,R). + * + * ### Mathematical Background + * + * Kleinman-Bylander non-local pseudopotential: + * \f[ + * V^{NL}(r, r') = \sum_{I} \sum_{m,m'} \beta_{m}^{I}(r)\, D_{mm'}^{I}\, \beta_{m'}^{I}(r') + * \f] + * + * In LCAO basis: + * \f[ + * H_{ij}^{NL} = \sum_{m,m'} \langle \phi_i | \beta_{m}^{I} \rangle + * D_{mm'}^{I} + * \langle \beta_{m'}^{I} | \phi_j \rangle + * \f] + * + * ### Derivative + * \f[ + * \frac{d H_{ij}^{NL}}{d \tau_I} + * = \sum_{m,m'} \Bigg[ + * \langle \frac{d \phi_i}{d \tau_I} | \beta_{m}^{I} \rangle D_{mm'}^{I} + * \langle \beta_{m'}^{I} | \phi_j \rangle + * + \langle \phi_i | \frac{d \beta_{m}^{I}}{d \tau_I} \rangle D_{mm'}^{I} + * \langle \beta_{m'}^{I} | \phi_j \rangle + * + \langle \phi_i | \beta_{m}^{I} \rangle D_{mm'}^{I} + * \langle \beta_{m'}^{I} | \frac{d \phi_j}{d \tau_I} \rangle + * \Bigg] + * \f] + * + * ### Hellmann-Feynman relation + * \f[ + * \langle \phi | \frac{d \beta}{d \tau} \rangle + * = - \langle \phi | \nabla \beta \rangle + * = \langle \nabla \phi | \beta \rangle + * \f] + * + * ### Computational strategy + * + * Intermediate vectors per direction d: + * \f[ + * tU[d] = \sum_{p_1,p_2} + * \langle \nabla \phi_{I1} | \beta \rangle_{p_1, d} + * \, D_{p_1, p_2} \, + * \langle \beta | \phi_{I2} \rangle_{p_2} + * \f] + * \f[ + * tV[d] = \sum_{p_1,p_2} + * \langle \phi_{I1} | \beta \rangle_{p_1} + * \, D_{p_1, p_2} \, + * \langle \beta | \nabla \phi_{I2} \rangle_{p_2, d} + * \f] + * + * Distribution: + * - dhR[d][I1] -= tU[d] (orbital 1 moves) + * - dhR[d][I2] -= tV[d] (orbital 2 moves) + * - dhR[d][I0] += tU[d]+tV[d] (Hellmann-Feynman: nucleus moves) + * + * @param[out] dhR Array of 3 vectors (x,y,z), each nat HContainer pointers. + */ void cal_dH(std::array*>, 3>& dhR); virtual void set_HR_fixed(void*) override; diff --git a/source/source_lcao/module_operator_lcao/nonlocal_dh.hpp b/source/source_lcao/module_operator_lcao/nonlocal_dh.cpp similarity index 91% rename from source/source_lcao/module_operator_lcao/nonlocal_dh.hpp rename to source/source_lcao/module_operator_lcao/nonlocal_dh.cpp index 432f08ad04..0085554be2 100644 --- a/source/source_lcao/module_operator_lcao/nonlocal_dh.hpp +++ b/source/source_lcao/module_operator_lcao/nonlocal_dh.cpp @@ -1,11 +1,13 @@ -#pragma once #include "nonlocal.h" -#include "operator_force_stress_utils.h" +#include "operator_fs_utils.h" #include "source_base/timer.h" namespace hamilt { +/** + * @see Nonlocal>::cal_dH in nonlocal.h for full documentation. + */ template void Nonlocal>::cal_dH(std::array*>, 3>& dhR) { @@ -43,7 +45,9 @@ void Nonlocal>::cal_dH(std::arrayucell->itia2iat(T1, I1); @@ -52,7 +56,9 @@ void Nonlocal>::cal_dH(std::arrayucell->itia2iat(T2, I2); @@ -69,7 +75,9 @@ void Nonlocal>::cal_dH(std::arrayinsert_pair(ap); + } } } } @@ -78,7 +86,9 @@ void Nonlocal>::cal_dH(std::arrayallocate(nullptr, true); + } } #pragma omp parallel @@ -112,7 +122,9 @@ void Nonlocal>::cal_dH(std::array>::cal_dH(std::array>::cal_dH(std::arrayucell->itia2iat(T1, I1); @@ -160,7 +176,9 @@ void Nonlocal>::cal_dH(std::arrayucell->itia2iat(T2, I2); @@ -183,7 +201,9 @@ void Nonlocal>::cal_dH(std::arrayget_pointer(), m1[1]->get_pointer(), m1[2]->get_pointer()}; double* p2[3] = {m2[0]->get_pointer(), m2[1]->get_pointer(), m2[2]->get_pointer()}; @@ -200,7 +220,9 @@ void Nonlocal>::cal_dH(std::array& nlm1 = it1->second; const size_t length = nlm1.size() / 4; const int iw1_row = static_cast(iw1l); @@ -209,7 +231,9 @@ void Nonlocal>::cal_dH(std::array& nlm2 = it2->second; const int iw2_col = static_cast(iw2l); @@ -254,4 +278,12 @@ void Nonlocal>::cal_dH(std::array>::cal_dH( + std::array*>, 3>& dhR); +template void Nonlocal, double>>::cal_dH( + std::array*>, 3>& dhR); +template void Nonlocal, std::complex>>::cal_dH( + std::array*>, 3>& dhR); + +} // namespace hamilt \ No newline at end of file diff --git a/source/source_lcao/module_operator_lcao/nonlocal_force_stress.hpp b/source/source_lcao/module_operator_lcao/nonlocal_fs.cpp similarity index 88% rename from source/source_lcao/module_operator_lcao/nonlocal_force_stress.hpp rename to source/source_lcao/module_operator_lcao/nonlocal_fs.cpp index f3eeb8aa00..83abeaf021 100644 --- a/source/source_lcao/module_operator_lcao/nonlocal_force_stress.hpp +++ b/source/source_lcao/module_operator_lcao/nonlocal_fs.cpp @@ -1,6 +1,5 @@ -#pragma once #include "nonlocal.h" -#include "operator_force_stress_utils.h" +#include "operator_fs_utils.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" @@ -219,7 +218,7 @@ void Nonlocal>::cal_force_stress(const bool cal_force, } template <> -void Nonlocal, std::complex>>::cal_force_IJR(const int& iat1, +inline void Nonlocal, std::complex>>::cal_force_IJR(const int& iat1, const int& iat2, const int& T0, const Parallel_Orbitals* paraV, @@ -284,7 +283,7 @@ void Nonlocal, std::complex>>::cal_for } template <> -void Nonlocal, std::complex>>::cal_stress_IJR(const int& iat1, +inline void Nonlocal, std::complex>>::cal_stress_IJR(const int& iat1, const int& iat2, const int& T0, const Parallel_Orbitals* paraV, @@ -454,4 +453,54 @@ void Nonlocal>::cal_stress_IJR(const int& iat1, } } +// explicit member function instantiations for cal_force_stress +template void Nonlocal>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); +template void Nonlocal, double>>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); +template void Nonlocal, std::complex>>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer>* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); + +// explicit member function instantiations for cal_force_IJR (generic template) +template void Nonlocal>::cal_force_IJR( + const int& iat1, const int& iat2, const int& T0, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + double* force1, double* force2); +template void Nonlocal, double>>::cal_force_IJR( + const int& iat1, const int& iat2, const int& T0, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + double* force1, double* force2); + +// explicit member function instantiations for cal_stress_IJR (generic template) +template void Nonlocal>::cal_stress_IJR( + const int& iat1, const int& iat2, const int& T0, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); +template void Nonlocal, double>>::cal_stress_IJR( + const int& iat1, const int& iat2, const int& T0, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); + } // namespace hamilt diff --git a/source/source_lcao/module_operator_lcao/op_exx_lcao.cpp b/source/source_lcao/module_operator_lcao/op_exx_lcao.cpp index 62eb7fe256..1a101973d4 100644 --- a/source/source_lcao/module_operator_lcao/op_exx_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/op_exx_lcao.cpp @@ -1,6 +1,15 @@ #ifdef __EXX #include "op_exx_lcao.h" #include "source_base/module_external/blacs_connector.h" +#include "source_base/parallel_reduce.h" +#include "source_hamilt/module_xc/xc_functional.h" +#include "source_io/module_parameter/parameter.h" +#include "source_io/module_restart/restart.h" +#include "source_io/module_restart/restart_exx_csr.h" +#include "source_hamilt/module_hcontainer/read_hcontainer.h" +#include "source_lcao/module_ri/Exx_LRI_interface.h" +#include "source_lcao/module_ri/RI_2D_Comm.h" +#include "source_lcao/module_rt/td_info.h" namespace hamilt { @@ -38,4 +47,696 @@ void OperatorEXX, std::complex>>::add_ } } // namespace hamilt + +// Begin content migrated from op_exx_lcao.hpp +namespace hamilt +{ +using TAC = std::pair>; + +// allocate according to the read-in HexxR, used in nscf +template +void reallocate_hcontainer(const std::vector>>>& Hexxs, + HContainer* hR, + const RI::Cell_Nearest* const cell_nearest) +{ + auto* pv = hR->get_paraV(); + bool need_allocate = false; + for (auto& Htmp1: Hexxs[0]) + { + const int& iat0 = Htmp1.first; + for (auto& Htmp2: Htmp1.second) + { + const int& iat1 = Htmp2.first.first; + if (pv->get_nrow_atom(iat0) > 0 && pv->get_ncol_atom(iat1) > 0) + { + const Abfs::Vector3_Order& R = RI_Util::array3_to_Vector3( + (cell_nearest ? cell_nearest->get_cell_nearest_discrete(iat0, iat1, Htmp2.first.second) + : Htmp2.first.second)); + BaseMatrix* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z); + if (HlocR == nullptr) + { // add R to HContainer + need_allocate = true; + AtomPair tmp(iat0, iat1, R.x, R.y, R.z, pv); + hR->insert_pair(tmp); + } + } + } + } + if (need_allocate) + { + hR->allocate(nullptr, true); + } +} + +/// allocate according to BvK cells, used in scf +template +void reallocate_hcontainer(const int nat, + HContainer* hR, + const std::array& Rs_period, + const RI::Cell_Nearest* const cell_nearest) +{ + auto* pv = hR->get_paraV(); + auto Rs = RI_Util::get_Born_von_Karmen_cells(Rs_period); + bool need_allocate = false; + for (int iat0 = 0; iat0 < nat; ++iat0) + { + for (int iat1 = 0; iat1 < nat; ++iat1) + { + // complete the atom pairs that has orbitals in this processor but not in hR due to the adj_list + // but adj_list is not enought for EXX, which is more nonlocal than Nonlocal + if (pv->get_nrow_atom(iat0) > 0 && pv->get_ncol_atom(iat1) > 0) + { + for (auto& cell: Rs) + { + const Abfs::Vector3_Order& R = RI_Util::array3_to_Vector3( + (cell_nearest ? cell_nearest->get_cell_nearest_discrete(iat0, iat1, cell) : cell)); + BaseMatrix* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z); + + if (HlocR == nullptr) + { // add R to HContainer + need_allocate = true; + AtomPair tmp(iat0, iat1, R.x, R.y, R.z, pv); + hR->insert_pair(tmp); + } + } + } + } + } + if (need_allocate) + { + hR->allocate(nullptr, true); + } +} + +template +OperatorEXX>::OperatorEXX( + HS_Matrix_K* hsk_in, + hamilt::HContainer* hR_in, + const UnitCell& ucell, + const K_Vectors& kv_in, + std::vector>>>* Hexxd_in, + std::vector>>>>* Hexxc_in, + Add_Hexx_Type add_hexx_type_in) + : OperatorLCAO(hsk_in, kv_in.kvec_d, hR_in), ucell(ucell), kv(kv_in), Hexxd(Hexxd_in), Hexxc(Hexxc_in), + add_hexx_type(add_hexx_type_in) +{ + this->cal_type = calculation_type::lcao_exx; + // This one-shot constructor never builds cell_nearest, so cal_dH() must not use it: + // the (d)Hexxs from LibRI are in native cells, and the dH output mirrors the H-term + // writer (write_h_exx_impl), which also passes a nullptr cell_nearest. + this->use_cell_nearest = false; +} + +template +OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, + HContainer* hR_in, + const UnitCell& ucell_in, + const K_Vectors& kv_in, + Exx_LRI_Interface* exd_in, + Exx_LRI_Interface>* exc_in, + Add_Hexx_Type add_hexx_type_in, + const int istep_in, + const bool restart_in) + : OperatorEXX>(hsk_in, + hR_in, + ucell_in, + kv_in, + exd_in ? &exd_in->get_Hexxs() : nullptr, + exc_in ? &exc_in->get_Hexxs() : nullptr, + add_hexx_type_in) +{ + this->exd = exd_in; + this->exc = exc_in; + const_cast(this->istep) = istep_in; + this->restart = restart_in; + ModuleBase::TITLE("OperatorEXX", "OperatorEXX"); + const Parallel_Orbitals* const pv = hR_in->get_paraV(); + + if (PARAM.inp.calculation == "nscf" && GlobalC::exx_info.info_global.cal_exx) + { // for nscf, calculate HexxR from the read-in DM, or read HexxR in + auto file_name_list_csr = []() -> std::vector { + std::vector file_name_list; + for (int irank = 0; irank < PARAM.globalv.nproc; ++irank) + { + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + file_name_list.push_back(PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(irank) + "_" + + std::to_string(is) + ".csr"); + } + } + return file_name_list; + }; + auto file_name_list_cereal = []() -> std::vector { + std::vector file_name_list; + for (int irank = 0; irank < PARAM.globalv.nproc; ++irank) + { + file_name_list.push_back("HexxR_" + std::to_string(irank)); + } + return file_name_list; + }; + auto check_exist = [](const std::vector& file_name_list) -> bool { + for (const std::string& file_name: file_name_list) + { + std::ifstream ifs(file_name); + if (!ifs.is_open()) + { + return false; + } + } + return true; + }; + + if (PARAM.inp.init_chg == "dm" || PARAM.inp.init_chg == "dm_no_renormalize") + { + // 1. cal Cs, Vs + if (GlobalC::exx_info.info_ri.real_number) + { + this->exd->cal_exx_ions(ucell, PARAM.inp.out_ri_cv); + } + else + { + this->exc->cal_exx_ions(ucell, PARAM.inp.out_ri_cv); + } + + // 2. read DM + const int nspin_dm = (PARAM.inp.nspin == 2) ? 2 : 1; + std::vector*> dmR_vec(nspin_dm); + for (int is = 0; is < nspin_dm; ++is) + { + const std::string dmfile + = PARAM.globalv.global_readin_dir + "/dmrs" + std::to_string(is + 1) + "_nao.csr"; + dmR_vec[is] = new hamilt::HContainer(const_cast(pv)); + hamilt::Read_HContainer reader_dm(dmR_vec[is], dmfile, PARAM.globalv.nlocal, &ucell); + reader_dm.read(); + } + + // 3. DM->Ds->Hexx (do not use symmetry for nscf) + XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); + if (GlobalC::exx_info.info_ri.real_number) + { + const auto& Ds = RI_2D_Comm::dm_container_to_Ds(dmR_vec, ucell, *pv, PARAM.inp.nspin); + this->exd->cal_exx_elec(Ds, ucell, *pv); + } + else + { + const auto& Ds = RI_2D_Comm::dm_container_to_Ds>(dmR_vec, + ucell, + *pv, + PARAM.inp.nspin); + this->exc->cal_exx_elec(Ds, ucell, *pv); + } + } + else // need to read HexxR + { + std::cout << " Attention: The number of MPI processes must be strictly identical between SCF and NSCF when " + "computing exact-exchange." + << std::endl; + if (check_exist(file_name_list_csr())) + { + // read HexxR first and reallocate hR according to the read-in HexxR + const std::string file_name_exx_csr + = PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(PARAM.globalv.myrank); + // Read HexxR in CSR format + if (GlobalC::exx_info.info_ri.real_number) + { + ModuleIO::read_Hexxs_csr(file_name_exx_csr, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd); + } + else + { + ModuleIO::read_Hexxs_csr(file_name_exx_csr, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc); + } + } + else if (check_exist(file_name_list_cereal())) + { + // Read HexxR in binary format (old version) + const std::string file_name_exx_cereal + = PARAM.globalv.global_readin_dir + "HexxR_" + std::to_string(PARAM.globalv.myrank); + std::ifstream ifs(file_name_exx_cereal, std::ios::binary); + if (!ifs) + { + ModuleBase::WARNING_QUIT("OperatorEXX", "Can't open EXX file < " + file_name_exx_cereal + " >."); + } + if (GlobalC::exx_info.info_ri.real_number) + { + ModuleIO::read_Hexxs_cereal(file_name_exx_cereal, *Hexxd); + } + else + { + ModuleIO::read_Hexxs_cereal(file_name_exx_cereal, *Hexxc); + } + } + else + { + ModuleBase::WARNING_QUIT("OperatorEXX", "Can't open EXX file in " + PARAM.globalv.global_readin_dir); + } + } + // reallocate hR according to Hexx(R) + if (this->add_hexx_type == Add_Hexx_Type::R) + { + if (GlobalC::exx_info.info_ri.real_number) + { + reallocate_hcontainer(*this->Hexxd, this->hR); + } + else + { + reallocate_hcontainer(*this->Hexxc, this->hR); + } + } + this->use_cell_nearest = false; + } + else + { // if scf and Add_Hexx_Type::R, init cell_nearest and reallocate hR according to BvK cells + if (this->add_hexx_type == Add_Hexx_Type::R) + { + // if k points has no shift, use cell_nearest to reduce the memory cost + this->use_cell_nearest = (ModuleBase::Vector3(std::fmod(this->kv.get_koffset(0), 1.0), + std::fmod(this->kv.get_koffset(1), 1.0), + std::fmod(this->kv.get_koffset(2), 1.0)) + .norm() + < 1e-10); + + const std::array Rs_period = {this->kv.nmp[0], this->kv.nmp[1], this->kv.nmp[2]}; + if (this->use_cell_nearest) + { + this->cell_nearest = init_cell_nearest(ucell, Rs_period); + reallocate_hcontainer(ucell.nat, this->hR, Rs_period, &this->cell_nearest); + } + else + { + reallocate_hcontainer(ucell.nat, this->hR, Rs_period); + } + } + + if (this->restart) + { /// Now only Hexx depends on DM, so we can directly read Hexx to reduce the computational cost. + /// If other operators depends on DM, we can also read DM and then calculate the operators to save the + /// memory to store operator terms. + + if (this->add_hexx_type == Add_Hexx_Type::k) + { + /// read in Hexx(k) + if (std::is_same::value) + { + this->Hexxd_k_load.resize(this->kv.get_nks()); + for (int ik = 0; ik < this->kv.get_nks(); ik++) + { + this->Hexxd_k_load[ik].resize(pv->get_local_size(), 0.0); + this->restart = GlobalC::restart.load_disk("Hexx", + ik, + pv->get_local_size(), + this->Hexxd_k_load[ik].data(), + false); + if (!this->restart) + { + break; + } + } + } + else + { + this->Hexxc_k_load.resize(this->kv.get_nks()); + for (int ik = 0; ik < this->kv.get_nks(); ik++) + { + this->Hexxc_k_load[ik].resize(pv->get_local_size(), 0.0); + this->restart = GlobalC::restart.load_disk("Hexx", + ik, + pv->get_local_size(), + this->Hexxc_k_load[ik].data(), + false); + if (!this->restart) + { + break; + } + } + } + } + else if (this->add_hexx_type == Add_Hexx_Type::R) + { + // read in Hexx(R) + const std::string restart_HR_path + = GlobalC::restart.folder + "HexxR" + std::to_string(PARAM.globalv.myrank); + int all_exist = 1; + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + std::ifstream ifs(restart_HR_path + "_" + std::to_string(is) + ".csr"); + if (!ifs) + { + all_exist = 0; + break; + } + } + // Add MPI communication to synchronize all_exist across processes +#ifdef __MPI + Parallel_Reduce::reduce_min(all_exist); +#endif + if (all_exist) + { + // Read HexxR in CSR format + if (GlobalC::exx_info.info_ri.real_number) + { + ModuleIO::read_Hexxs_csr(restart_HR_path, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd); + } + else + { + ModuleIO::read_Hexxs_csr(restart_HR_path, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc); + } + } + else + { + // Read HexxR in binary format (old version) + const std::string restart_HR_path_cereal + = GlobalC::restart.folder + "HexxR_" + std::to_string(PARAM.globalv.myrank); + std::ifstream ifs(restart_HR_path_cereal, std::ios::binary); + int all_exist_cereal = ifs ? 1 : 0; +#ifdef __MPI + Parallel_Reduce::reduce_min(all_exist_cereal); +#endif + if (!all_exist_cereal) + { + // no HexxR file in CSR or binary format + this->restart = false; + } + else + { + if (GlobalC::exx_info.info_ri.real_number) + { + ModuleIO::read_Hexxs_cereal(restart_HR_path_cereal, *Hexxd); + } + else + { + ModuleIO::read_Hexxs_cereal(restart_HR_path_cereal, *Hexxc); + } + } + } + } + + if (!this->restart) + { + std::cout << "WARNING: Hexx not found, restart from the non-exx loop." << std::endl + << "If the loaded charge density is EXX-solved, this may lead to poor convergence." + << std::endl; + } + GlobalC::restart.info_load.load_H_finish = this->restart; + } + } +} +template +void OperatorEXX>::contributeHR() +{ + ModuleBase::TITLE("OperatorEXX", "contributeHR"); + // Peize Lin add 2016-12-03 + + // 1. For NSCF + if (PARAM.inp.calculation == "nscf") + { + // Do nothing here, allow the code to proceed and calculate EXX. + } + // 2. For the first ionic step of SCF, relaxation, or MD: + else if (this->istep == 0) + { + const int two_level_step + = GlobalC::exx_info.info_ri.real_number ? this->exd->get_two_level_step() : this->exc->get_two_level_step(); + + // Check if we are in the pre-convergence stage of the two-level SCF (i.e., the pure GGA loop) + bool in_gga_pre_loop = (two_level_step == 0); + + // Check if a high-quality initial guess is missing (neither reading wavefunctions from a file nor restarting) + bool lacks_good_guess = (PARAM.inp.init_wfc != "file" && !this->restart); + + // If in the pre-convergence loop and lacking a good initial guess, skip adding the EXX contribution + if (in_gga_pre_loop && lacks_good_guess) + { + return; // In the non-EXX loop, skip adding EXX contribution + } + } + // 3. For subsequent ionic steps (istep > 0), add EXX normally + + if (this->add_hexx_type == Add_Hexx_Type::k) + { + return; + } + + if (XC_Functional::get_func_type() == 4 || XC_Functional::get_func_type() == 5) + { + // add H(R) normally + if (GlobalC::exx_info.info_ri.real_number) + { + RI_2D_Comm::add_HexxR(this->current_spin, + GlobalC::exx_info.info_global.hybrid_alpha, + *this->Hexxd, + *this->hR->get_paraV(), + PARAM.globalv.npol, + *this->hR, + this->use_cell_nearest ? &this->cell_nearest : nullptr); + } + else + { + RI_2D_Comm::add_HexxR(this->current_spin, + GlobalC::exx_info.info_global.hybrid_alpha, + *this->Hexxc, + *this->hR->get_paraV(), + PARAM.globalv.npol, + *this->hR, + this->use_cell_nearest ? &this->cell_nearest : nullptr); + } + } + if (PARAM.inp.nspin == 2) + { + this->current_spin = 1 - this->current_spin; + } +} + +template +void OperatorEXX>::contributeHk(int ik) +{ + ModuleBase::TITLE("OperatorEXX", "constributeHk"); + const bool has_workflow = GlobalC::exx_info.info_ri.real_number ? (this->exd != nullptr) : (this->exc != nullptr); + int two_level_step = 0; + if (has_workflow) + { + two_level_step + = GlobalC::exx_info.info_ri.real_number ? this->exd->get_two_level_step() : this->exc->get_two_level_step(); + } + + // Peize Lin add 2016-12-03 + + // Taoni Bao add 2026-05-15 + // In RT-TDDFT, contributeHk is used, but two_level_step is reset to 0 at each ionic step. + // In order to add EXX correctly in for istep > 0, this->istep == 0 is needed to avoid skipping EXX calculation. + // 1. For NSCF + if (PARAM.inp.calculation == "nscf" || !has_workflow) + { + // Do nothing here, allow the code to proceed and calculate EXX. + } + // 2. For the first ionic step: + else if (this->istep == 0) + { + // If EXX is once turned on (two_level_step > 0), let OperatorEXX remember this + if (two_level_step > 0) + { + this->initial_gga_done = true; + } + + // Check if we are in the pre-convergence stage of the two-level SCF (i.e., the pure GGA loop) + bool in_gga_pre_loop = (two_level_step == 0); + + // Check if a high-quality initial guess is missing + bool lacks_good_guess = (!this->restart); + + // If in the pre-convergence loop and lacking a good initial guess, skip adding the EXX contribution + // Taoni Bao add 2026-05-18, only skip EXX if initial GGA loop is not done + // Fix RT-TDDFT EXX missing problem in the evolution + if (in_gga_pre_loop && lacks_good_guess && !this->initial_gga_done) + { + return; // In the non-EXX loop, skip adding EXX contribution + } + } + // 3. For subsequent ionic steps (istep > 0), add EXX normally + + if (this->add_hexx_type == Add_Hexx_Type::R) + { + OperatorLCAO::contributeHk(ik); + } + + if (XC_Functional::get_func_type() == 4 || XC_Functional::get_func_type() == 5) + { + if (this->restart) + { + if (two_level_step == 0) + { + this->add_loaded_Hexx(ik); + return; + } + else // clear loaded Hexx and release memory + { + if (this->Hexxd_k_load.size() > 0) + { + this->Hexxd_k_load.clear(); + this->Hexxd_k_load.shrink_to_fit(); + } + else if (this->Hexxc_k_load.size() > 0) + { + this->Hexxc_k_load.clear(); + this->Hexxc_k_load.shrink_to_fit(); + } + } + } + // cal H(k) from H(R) normally + if (PARAM.inp.esolver_type == "tddft" && PARAM.inp.td_stype == 2) + { + RI_2D_Comm::add_Hexx_td(ucell, + this->kv, + ik, + GlobalC::exx_info.info_global.hybrid_alpha, + *this->Hexxc, + *this->hR->get_paraV(), + TD_info::td_vel_op->cart_At, + TD_info::td_vel_op->get_phase_hybrid(), + this->hsk->get_hk()); + } + else + { + if (GlobalC::exx_info.info_ri.real_number) + { + RI_2D_Comm::add_Hexx(ucell, + this->kv, + ik, + GlobalC::exx_info.info_global.hybrid_alpha, + *this->Hexxd, + *this->hR->get_paraV(), + this->hsk->get_hk()); + } + else + { + RI_2D_Comm::add_Hexx(ucell, + this->kv, + ik, + GlobalC::exx_info.info_global.hybrid_alpha, + *this->Hexxc, + *this->hR->get_paraV(), + this->hsk->get_hk()); + } + } + } +} + +template +template +void OperatorEXX>::cal_dH( + const int ispin, + std::array*>, 3>& dhR, + const std::array>>>>, 3>& dHexxs) +{ + // dhR is the set of per-atom-I HContainers to fill (not this->hR, which may be a dummy here). + const Parallel_Orbitals* const paraV = dhR[0][0]->get_paraV(); + const RI::Cell_Nearest* const cell_nearest + = this->use_cell_nearest ? &this->cell_nearest : nullptr; + for (int idir = 0; idir < 3; ++idir) + { + for (int iat = 0; iat < ucell.nat; ++iat) + { + // add_HexxR only fills existing matrices, so first allocate the atom-pair + // structure of this per-I container from the exx-form data (same cell mapping). + reallocate_hcontainer(dHexxs[idir][iat], dhR[idir][iat], cell_nearest); + RI_2D_Comm::add_HexxR(ispin, + GlobalC::exx_info.info_global.hybrid_alpha, + dHexxs[idir][iat], + *paraV, + PARAM.globalv.npol, + *dhR[idir][iat], + cell_nearest); + } + } +} + +// explicit member function instantiations for constructors +template OperatorEXX>::OperatorEXX( + HS_Matrix_K*, HContainer*, const UnitCell&, const K_Vectors&, + std::vector>>>*, + std::vector>>>>*, + Add_Hexx_Type); +template OperatorEXX, double>>::OperatorEXX( + HS_Matrix_K>*, HContainer*, const UnitCell&, const K_Vectors&, + std::vector>>>*, + std::vector>>>>*, + Add_Hexx_Type); +template OperatorEXX, std::complex>>::OperatorEXX( + HS_Matrix_K>*, HContainer>*, const UnitCell&, const K_Vectors&, + std::vector>>>*, + std::vector>>>>*, + Add_Hexx_Type); + +// explicit member function instantiations for second constructor +template OperatorEXX>::OperatorEXX( + HS_Matrix_K*, HContainer*, const UnitCell&, const K_Vectors&, + Exx_LRI_Interface*, Exx_LRI_Interface>*, + Add_Hexx_Type, const int, const bool); +template OperatorEXX, double>>::OperatorEXX( + HS_Matrix_K>*, HContainer*, const UnitCell&, const K_Vectors&, + Exx_LRI_Interface, double>*, Exx_LRI_Interface, std::complex>*, + Add_Hexx_Type, const int, const bool); +template OperatorEXX, std::complex>>::OperatorEXX( + HS_Matrix_K>*, HContainer>*, const UnitCell&, const K_Vectors&, + Exx_LRI_Interface, double>*, Exx_LRI_Interface, std::complex>*, + Add_Hexx_Type, const int, const bool); + +// explicit member function instantiations for contributeHR +template void OperatorEXX>::contributeHR(); +template void OperatorEXX, double>>::contributeHR(); +template void OperatorEXX, std::complex>>::contributeHR(); + +// explicit member function instantiations for contributeHk +template void OperatorEXX>::contributeHk(int); +template void OperatorEXX, double>>::contributeHk(int); +template void OperatorEXX, std::complex>>::contributeHk(int); + +// explicit member function instantiations for cal_dH (template member function) +template void OperatorEXX>::cal_dH( + const int, std::array*>, 3>&, + const std::array>>>>, 3>&); +template void OperatorEXX>::cal_dH>( + const int, std::array*>, 3>&, + const std::array>>>>>, 3>&); +template void OperatorEXX, double>>::cal_dH( + const int, std::array*>, 3>&, + const std::array>>>>, 3>&); +template void OperatorEXX, double>>::cal_dH>( + const int, std::array*>, 3>&, + const std::array>>>>>, 3>&); +template void OperatorEXX, std::complex>>::cal_dH( + const int, std::array*>, 3>&, + const std::array>>>>, 3>&); +template void OperatorEXX, std::complex>>::cal_dH>( + const int, std::array*>, 3>&, + const std::array>>>>>, 3>&); + +// explicit instantiations for reallocate_hcontainer (first overload) +template void reallocate_hcontainer( + const std::vector>>>&, + HContainer*, + const RI::Cell_Nearest* const); +template void reallocate_hcontainer, double>( + const std::vector>>>>&, + HContainer*, + const RI::Cell_Nearest* const); +template void reallocate_hcontainer>( + const std::vector>>>&, + HContainer>*, + const RI::Cell_Nearest* const); +template void reallocate_hcontainer, std::complex>( + const std::vector>>>>&, + HContainer>*, + const RI::Cell_Nearest* const); + +// explicit instantiations for reallocate_hcontainer (second overload) +template void reallocate_hcontainer( + const int, HContainer*, const std::array&, + const RI::Cell_Nearest* const); +template void reallocate_hcontainer>( + const int, HContainer>*, const std::array&, + const RI::Cell_Nearest* const); + +} // namespace hamilt + +// End content migrated from op_exx_lcao.hpp #endif \ No newline at end of file diff --git a/source/source_lcao/module_operator_lcao/op_exx_lcao.h b/source/source_lcao/module_operator_lcao/op_exx_lcao.h index 0fa20ae449..1cb25bf054 100644 --- a/source/source_lcao/module_operator_lcao/op_exx_lcao.h +++ b/source/source_lcao/module_operator_lcao/op_exx_lcao.h @@ -122,5 +122,4 @@ void reallocate_hcontainer(const int nat, } // namespace hamilt #endif // __EXX -#include "op_exx_lcao.hpp" #endif // OPEXXLCAO_H \ No newline at end of file diff --git a/source/source_lcao/module_operator_lcao/op_exx_lcao.hpp b/source/source_lcao/module_operator_lcao/op_exx_lcao.hpp deleted file mode 100644 index 73130524d0..0000000000 --- a/source/source_lcao/module_operator_lcao/op_exx_lcao.hpp +++ /dev/null @@ -1,618 +0,0 @@ -#ifndef OPEXXLCAO_HPP -#define OPEXXLCAO_HPP -#ifdef __EXX - -#include "op_exx_lcao.h" -#include "source_base/parallel_reduce.h" -#include "source_hamilt/module_xc/xc_functional.h" -#include "source_io/module_parameter/parameter.h" -#include "source_io/module_restart/restart.h" -#include "source_io/module_restart/restart_exx_csr.h" -#include "source_hamilt/module_hcontainer/read_hcontainer.h" -#include "source_lcao/module_ri/Exx_LRI_interface.h" -#include "source_lcao/module_ri/RI_2D_Comm.h" -#include "source_lcao/module_rt/td_info.h" - -namespace hamilt -{ -using TAC = std::pair>; - -// allocate according to the read-in HexxR, used in nscf -template -void reallocate_hcontainer(const std::vector>>>& Hexxs, - HContainer* hR, - const RI::Cell_Nearest* const cell_nearest) -{ - auto* pv = hR->get_paraV(); - bool need_allocate = false; - for (auto& Htmp1: Hexxs[0]) - { - const int& iat0 = Htmp1.first; - for (auto& Htmp2: Htmp1.second) - { - const int& iat1 = Htmp2.first.first; - if (pv->get_nrow_atom(iat0) > 0 && pv->get_ncol_atom(iat1) > 0) - { - const Abfs::Vector3_Order& R = RI_Util::array3_to_Vector3( - (cell_nearest ? cell_nearest->get_cell_nearest_discrete(iat0, iat1, Htmp2.first.second) - : Htmp2.first.second)); - BaseMatrix* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z); - if (HlocR == nullptr) - { // add R to HContainer - need_allocate = true; - AtomPair tmp(iat0, iat1, R.x, R.y, R.z, pv); - hR->insert_pair(tmp); - } - } - } - } - if (need_allocate) - { - hR->allocate(nullptr, true); - } -} - -/// allocate according to BvK cells, used in scf -template -void reallocate_hcontainer(const int nat, - HContainer* hR, - const std::array& Rs_period, - const RI::Cell_Nearest* const cell_nearest) -{ - auto* pv = hR->get_paraV(); - auto Rs = RI_Util::get_Born_von_Karmen_cells(Rs_period); - bool need_allocate = false; - for (int iat0 = 0; iat0 < nat; ++iat0) - { - for (int iat1 = 0; iat1 < nat; ++iat1) - { - // complete the atom pairs that has orbitals in this processor but not in hR due to the adj_list - // but adj_list is not enought for EXX, which is more nonlocal than Nonlocal - if (pv->get_nrow_atom(iat0) > 0 && pv->get_ncol_atom(iat1) > 0) - { - for (auto& cell: Rs) - { - const Abfs::Vector3_Order& R = RI_Util::array3_to_Vector3( - (cell_nearest ? cell_nearest->get_cell_nearest_discrete(iat0, iat1, cell) : cell)); - BaseMatrix* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z); - - if (HlocR == nullptr) - { // add R to HContainer - need_allocate = true; - AtomPair tmp(iat0, iat1, R.x, R.y, R.z, pv); - hR->insert_pair(tmp); - } - } - } - } - } - if (need_allocate) - { - hR->allocate(nullptr, true); - } -} - -template -OperatorEXX>::OperatorEXX( - HS_Matrix_K* hsk_in, - hamilt::HContainer* hR_in, - const UnitCell& ucell, - const K_Vectors& kv_in, - std::vector>>>* Hexxd_in, - std::vector>>>>* Hexxc_in, - Add_Hexx_Type add_hexx_type_in) - : OperatorLCAO(hsk_in, kv_in.kvec_d, hR_in), ucell(ucell), kv(kv_in), Hexxd(Hexxd_in), Hexxc(Hexxc_in), - add_hexx_type(add_hexx_type_in) -{ - this->cal_type = calculation_type::lcao_exx; - // This one-shot constructor never builds cell_nearest, so cal_dH() must not use it: - // the (d)Hexxs from LibRI are in native cells, and the dH output mirrors the H-term - // writer (write_h_exx_impl), which also passes a nullptr cell_nearest. - this->use_cell_nearest = false; -} - -template -OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, - HContainer* hR_in, - const UnitCell& ucell_in, - const K_Vectors& kv_in, - Exx_LRI_Interface* exd_in, - Exx_LRI_Interface>* exc_in, - Add_Hexx_Type add_hexx_type_in, - const int istep_in, - const bool restart_in) - : OperatorEXX>(hsk_in, - hR_in, - ucell_in, - kv_in, - exd_in ? &exd_in->get_Hexxs() : nullptr, - exc_in ? &exc_in->get_Hexxs() : nullptr, - add_hexx_type_in) -{ - this->exd = exd_in; - this->exc = exc_in; - const_cast(this->istep) = istep_in; - this->restart = restart_in; - ModuleBase::TITLE("OperatorEXX", "OperatorEXX"); - const Parallel_Orbitals* const pv = hR_in->get_paraV(); - - if (PARAM.inp.calculation == "nscf" && GlobalC::exx_info.info_global.cal_exx) - { // for nscf, calculate HexxR from the read-in DM, or read HexxR in - auto file_name_list_csr = []() -> std::vector { - std::vector file_name_list; - for (int irank = 0; irank < PARAM.globalv.nproc; ++irank) - { - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - file_name_list.push_back(PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(irank) + "_" - + std::to_string(is) + ".csr"); - } - } - return file_name_list; - }; - auto file_name_list_cereal = []() -> std::vector { - std::vector file_name_list; - for (int irank = 0; irank < PARAM.globalv.nproc; ++irank) - { - file_name_list.push_back("HexxR_" + std::to_string(irank)); - } - return file_name_list; - }; - auto check_exist = [](const std::vector& file_name_list) -> bool { - for (const std::string& file_name: file_name_list) - { - std::ifstream ifs(file_name); - if (!ifs.is_open()) - { - return false; - } - } - return true; - }; - - if (PARAM.inp.init_chg == "dm" || PARAM.inp.init_chg == "dm_no_renormalize") - { - // 1. cal Cs, Vs - if (GlobalC::exx_info.info_ri.real_number) - { - this->exd->cal_exx_ions(ucell, PARAM.inp.out_ri_cv); - } - else - { - this->exc->cal_exx_ions(ucell, PARAM.inp.out_ri_cv); - } - - // 2. read DM - const int nspin_dm = (PARAM.inp.nspin == 2) ? 2 : 1; - std::vector*> dmR_vec(nspin_dm); - for (int is = 0; is < nspin_dm; ++is) - { - const std::string dmfile - = PARAM.globalv.global_readin_dir + "/dmrs" + std::to_string(is + 1) + "_nao.csr"; - dmR_vec[is] = new hamilt::HContainer(const_cast(pv)); - hamilt::Read_HContainer reader_dm(dmR_vec[is], dmfile, PARAM.globalv.nlocal, &ucell); - reader_dm.read(); - } - - // 3. DM->Ds->Hexx (do not use symmetry for nscf) - XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); - if (GlobalC::exx_info.info_ri.real_number) - { - const auto& Ds = RI_2D_Comm::dm_container_to_Ds(dmR_vec, ucell, *pv, PARAM.inp.nspin); - this->exd->cal_exx_elec(Ds, ucell, *pv); - } - else - { - const auto& Ds = RI_2D_Comm::dm_container_to_Ds>(dmR_vec, - ucell, - *pv, - PARAM.inp.nspin); - this->exc->cal_exx_elec(Ds, ucell, *pv); - } - } - else // need to read HexxR - { - std::cout << " Attention: The number of MPI processes must be strictly identical between SCF and NSCF when " - "computing exact-exchange." - << std::endl; - if (check_exist(file_name_list_csr())) - { - // read HexxR first and reallocate hR according to the read-in HexxR - const std::string file_name_exx_csr - = PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(PARAM.globalv.myrank); - // Read HexxR in CSR format - if (GlobalC::exx_info.info_ri.real_number) - { - ModuleIO::read_Hexxs_csr(file_name_exx_csr, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd); - } - else - { - ModuleIO::read_Hexxs_csr(file_name_exx_csr, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc); - } - } - else if (check_exist(file_name_list_cereal())) - { - // Read HexxR in binary format (old version) - const std::string file_name_exx_cereal - = PARAM.globalv.global_readin_dir + "HexxR_" + std::to_string(PARAM.globalv.myrank); - std::ifstream ifs(file_name_exx_cereal, std::ios::binary); - if (!ifs) - { - ModuleBase::WARNING_QUIT("OperatorEXX", "Can't open EXX file < " + file_name_exx_cereal + " >."); - } - if (GlobalC::exx_info.info_ri.real_number) - { - ModuleIO::read_Hexxs_cereal(file_name_exx_cereal, *Hexxd); - } - else - { - ModuleIO::read_Hexxs_cereal(file_name_exx_cereal, *Hexxc); - } - } - else - { - ModuleBase::WARNING_QUIT("OperatorEXX", "Can't open EXX file in " + PARAM.globalv.global_readin_dir); - } - } - // reallocate hR according to Hexx(R) - if (this->add_hexx_type == Add_Hexx_Type::R) - { - if (GlobalC::exx_info.info_ri.real_number) - { - reallocate_hcontainer(*this->Hexxd, this->hR); - } - else - { - reallocate_hcontainer(*this->Hexxc, this->hR); - } - } - this->use_cell_nearest = false; - } - else - { // if scf and Add_Hexx_Type::R, init cell_nearest and reallocate hR according to BvK cells - if (this->add_hexx_type == Add_Hexx_Type::R) - { - // if k points has no shift, use cell_nearest to reduce the memory cost - this->use_cell_nearest = (ModuleBase::Vector3(std::fmod(this->kv.get_koffset(0), 1.0), - std::fmod(this->kv.get_koffset(1), 1.0), - std::fmod(this->kv.get_koffset(2), 1.0)) - .norm() - < 1e-10); - - const std::array Rs_period = {this->kv.nmp[0], this->kv.nmp[1], this->kv.nmp[2]}; - if (this->use_cell_nearest) - { - this->cell_nearest = init_cell_nearest(ucell, Rs_period); - reallocate_hcontainer(ucell.nat, this->hR, Rs_period, &this->cell_nearest); - } - else - { - reallocate_hcontainer(ucell.nat, this->hR, Rs_period); - } - } - - if (this->restart) - { /// Now only Hexx depends on DM, so we can directly read Hexx to reduce the computational cost. - /// If other operators depends on DM, we can also read DM and then calculate the operators to save the - /// memory to store operator terms. - - if (this->add_hexx_type == Add_Hexx_Type::k) - { - /// read in Hexx(k) - if (std::is_same::value) - { - this->Hexxd_k_load.resize(this->kv.get_nks()); - for (int ik = 0; ik < this->kv.get_nks(); ik++) - { - this->Hexxd_k_load[ik].resize(pv->get_local_size(), 0.0); - this->restart = GlobalC::restart.load_disk("Hexx", - ik, - pv->get_local_size(), - this->Hexxd_k_load[ik].data(), - false); - if (!this->restart) - { - break; - } - } - } - else - { - this->Hexxc_k_load.resize(this->kv.get_nks()); - for (int ik = 0; ik < this->kv.get_nks(); ik++) - { - this->Hexxc_k_load[ik].resize(pv->get_local_size(), 0.0); - this->restart = GlobalC::restart.load_disk("Hexx", - ik, - pv->get_local_size(), - this->Hexxc_k_load[ik].data(), - false); - if (!this->restart) - { - break; - } - } - } - } - else if (this->add_hexx_type == Add_Hexx_Type::R) - { - // read in Hexx(R) - const std::string restart_HR_path - = GlobalC::restart.folder + "HexxR" + std::to_string(PARAM.globalv.myrank); - int all_exist = 1; - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - std::ifstream ifs(restart_HR_path + "_" + std::to_string(is) + ".csr"); - if (!ifs) - { - all_exist = 0; - break; - } - } - // Add MPI communication to synchronize all_exist across processes -#ifdef __MPI - Parallel_Reduce::reduce_min(all_exist); -#endif - if (all_exist) - { - // Read HexxR in CSR format - if (GlobalC::exx_info.info_ri.real_number) - { - ModuleIO::read_Hexxs_csr(restart_HR_path, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd); - } - else - { - ModuleIO::read_Hexxs_csr(restart_HR_path, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc); - } - } - else - { - // Read HexxR in binary format (old version) - const std::string restart_HR_path_cereal - = GlobalC::restart.folder + "HexxR_" + std::to_string(PARAM.globalv.myrank); - std::ifstream ifs(restart_HR_path_cereal, std::ios::binary); - int all_exist_cereal = ifs ? 1 : 0; -#ifdef __MPI - Parallel_Reduce::reduce_min(all_exist_cereal); -#endif - if (!all_exist_cereal) - { - // no HexxR file in CSR or binary format - this->restart = false; - } - else - { - if (GlobalC::exx_info.info_ri.real_number) - { - ModuleIO::read_Hexxs_cereal(restart_HR_path_cereal, *Hexxd); - } - else - { - ModuleIO::read_Hexxs_cereal(restart_HR_path_cereal, *Hexxc); - } - } - } - } - - if (!this->restart) - { - std::cout << "WARNING: Hexx not found, restart from the non-exx loop." << std::endl - << "If the loaded charge density is EXX-solved, this may lead to poor convergence." - << std::endl; - } - GlobalC::restart.info_load.load_H_finish = this->restart; - } - } -} -template -void OperatorEXX>::contributeHR() -{ - ModuleBase::TITLE("OperatorEXX", "contributeHR"); - // Peize Lin add 2016-12-03 - - // 1. For NSCF - if (PARAM.inp.calculation == "nscf") - { - // Do nothing here, allow the code to proceed and calculate EXX. - } - // 2. For the first ionic step of SCF, relaxation, or MD: - else if (this->istep == 0) - { - const int two_level_step - = GlobalC::exx_info.info_ri.real_number ? this->exd->get_two_level_step() : this->exc->get_two_level_step(); - - // Check if we are in the pre-convergence stage of the two-level SCF (i.e., the pure GGA loop) - bool in_gga_pre_loop = (two_level_step == 0); - - // Check if a high-quality initial guess is missing (neither reading wavefunctions from a file nor restarting) - bool lacks_good_guess = (PARAM.inp.init_wfc != "file" && !this->restart); - - // If in the pre-convergence loop and lacking a good initial guess, skip adding the EXX contribution - if (in_gga_pre_loop && lacks_good_guess) - { - return; // In the non-EXX loop, skip adding EXX contribution - } - } - // 3. For subsequent ionic steps (istep > 0), add EXX normally - - if (this->add_hexx_type == Add_Hexx_Type::k) - { - return; - } - - if (XC_Functional::get_func_type() == 4 || XC_Functional::get_func_type() == 5) - { - // add H(R) normally - if (GlobalC::exx_info.info_ri.real_number) - { - RI_2D_Comm::add_HexxR(this->current_spin, - GlobalC::exx_info.info_global.hybrid_alpha, - *this->Hexxd, - *this->hR->get_paraV(), - PARAM.globalv.npol, - *this->hR, - this->use_cell_nearest ? &this->cell_nearest : nullptr); - } - else - { - RI_2D_Comm::add_HexxR(this->current_spin, - GlobalC::exx_info.info_global.hybrid_alpha, - *this->Hexxc, - *this->hR->get_paraV(), - PARAM.globalv.npol, - *this->hR, - this->use_cell_nearest ? &this->cell_nearest : nullptr); - } - } - if (PARAM.inp.nspin == 2) - { - this->current_spin = 1 - this->current_spin; - } -} - -template -void OperatorEXX>::contributeHk(int ik) -{ - ModuleBase::TITLE("OperatorEXX", "constributeHk"); - const bool has_workflow = GlobalC::exx_info.info_ri.real_number ? (this->exd != nullptr) : (this->exc != nullptr); - int two_level_step = 0; - if (has_workflow) - { - two_level_step - = GlobalC::exx_info.info_ri.real_number ? this->exd->get_two_level_step() : this->exc->get_two_level_step(); - } - - // Peize Lin add 2016-12-03 - - // Taoni Bao add 2026-05-15 - // In RT-TDDFT, contributeHk is used, but two_level_step is reset to 0 at each ionic step. - // In order to add EXX correctly in for istep > 0, this->istep == 0 is needed to avoid skipping EXX calculation. - // 1. For NSCF - if (PARAM.inp.calculation == "nscf" || !has_workflow) - { - // Do nothing here, allow the code to proceed and calculate EXX. - } - // 2. For the first ionic step: - else if (this->istep == 0) - { - // If EXX is once turned on (two_level_step > 0), let OperatorEXX remember this - if (two_level_step > 0) - { - this->initial_gga_done = true; - } - - // Check if we are in the pre-convergence stage of the two-level SCF (i.e., the pure GGA loop) - bool in_gga_pre_loop = (two_level_step == 0); - - // Check if a high-quality initial guess is missing - bool lacks_good_guess = (!this->restart); - - // If in the pre-convergence loop and lacking a good initial guess, skip adding the EXX contribution - // Taoni Bao add 2026-05-18, only skip EXX if initial GGA loop is not done - // Fix RT-TDDFT EXX missing problem in the evolution - if (in_gga_pre_loop && lacks_good_guess && !this->initial_gga_done) - { - return; // In the non-EXX loop, skip adding EXX contribution - } - } - // 3. For subsequent ionic steps (istep > 0), add EXX normally - - if (this->add_hexx_type == Add_Hexx_Type::R) - { - OperatorLCAO::contributeHk(ik); - } - - if (XC_Functional::get_func_type() == 4 || XC_Functional::get_func_type() == 5) - { - if (this->restart) - { - if (two_level_step == 0) - { - this->add_loaded_Hexx(ik); - return; - } - else // clear loaded Hexx and release memory - { - if (this->Hexxd_k_load.size() > 0) - { - this->Hexxd_k_load.clear(); - this->Hexxd_k_load.shrink_to_fit(); - } - else if (this->Hexxc_k_load.size() > 0) - { - this->Hexxc_k_load.clear(); - this->Hexxc_k_load.shrink_to_fit(); - } - } - } - // cal H(k) from H(R) normally - if (PARAM.inp.esolver_type == "tddft" && PARAM.inp.td_stype == 2) - { - RI_2D_Comm::add_Hexx_td(ucell, - this->kv, - ik, - GlobalC::exx_info.info_global.hybrid_alpha, - *this->Hexxc, - *this->hR->get_paraV(), - TD_info::td_vel_op->cart_At, - TD_info::td_vel_op->get_phase_hybrid(), - this->hsk->get_hk()); - } - else - { - if (GlobalC::exx_info.info_ri.real_number) - { - RI_2D_Comm::add_Hexx(ucell, - this->kv, - ik, - GlobalC::exx_info.info_global.hybrid_alpha, - *this->Hexxd, - *this->hR->get_paraV(), - this->hsk->get_hk()); - } - else - { - RI_2D_Comm::add_Hexx(ucell, - this->kv, - ik, - GlobalC::exx_info.info_global.hybrid_alpha, - *this->Hexxc, - *this->hR->get_paraV(), - this->hsk->get_hk()); - } - } - } -} - -template -template -void OperatorEXX>::cal_dH( - const int ispin, - std::array*>, 3>& dhR, - const std::array>>>>, 3>& dHexxs) -{ - // dhR is the set of per-atom-I HContainers to fill (not this->hR, which may be a dummy here). - const Parallel_Orbitals* const paraV = dhR[0][0]->get_paraV(); - const RI::Cell_Nearest* const cell_nearest - = this->use_cell_nearest ? &this->cell_nearest : nullptr; - for (int idir = 0; idir < 3; ++idir) - { - for (int iat = 0; iat < ucell.nat; ++iat) - { - // add_HexxR only fills existing matrices, so first allocate the atom-pair - // structure of this per-I container from the exx-form data (same cell mapping). - reallocate_hcontainer(dHexxs[idir][iat], dhR[idir][iat], cell_nearest); - RI_2D_Comm::add_HexxR(ispin, - GlobalC::exx_info.info_global.hybrid_alpha, - dHexxs[idir][iat], - *paraV, - PARAM.globalv.npol, - *dhR[idir][iat], - cell_nearest); - } - } -} - -} // namespace hamilt -#endif // __EXX -#endif // OPEXXLCAO_HPP diff --git a/source/source_lcao/module_operator_lcao/operator_force_stress_utils.cpp b/source/source_lcao/module_operator_lcao/operator_fs_utils.cpp similarity index 96% rename from source/source_lcao/module_operator_lcao/operator_force_stress_utils.cpp rename to source/source_lcao/module_operator_lcao/operator_fs_utils.cpp index c485727868..34ccc7c4ee 100644 --- a/source/source_lcao/module_operator_lcao/operator_force_stress_utils.cpp +++ b/source/source_lcao/module_operator_lcao/operator_fs_utils.cpp @@ -1,4 +1,4 @@ -#include "operator_force_stress_utils.h" +#include "operator_fs_utils.h" #include "source_base/parallel_reduce.h" namespace OperatorForceStress { diff --git a/source/source_lcao/module_operator_lcao/operator_force_stress_utils.h b/source/source_lcao/module_operator_lcao/operator_fs_utils.h similarity index 96% rename from source/source_lcao/module_operator_lcao/operator_force_stress_utils.h rename to source/source_lcao/module_operator_lcao/operator_fs_utils.h index 26ec0245fb..e880aac43b 100644 --- a/source/source_lcao/module_operator_lcao/operator_force_stress_utils.h +++ b/source/source_lcao/module_operator_lcao/operator_fs_utils.h @@ -1,5 +1,5 @@ -#ifndef OPERATOR_FORCE_STRESS_UTILS_H -#define OPERATOR_FORCE_STRESS_UTILS_H +#ifndef OPERATOR_FS_UTILS_H +#define OPERATOR_FS_UTILS_H #include "source_base/matrix.h" #include "source_cell/unitcell.h" @@ -121,4 +121,4 @@ void finalize_force_stress( } // namespace OperatorForceStress -#endif // OPERATOR_FORCE_STRESS_UTILS_H +#endif // OPERATOR_FS_UTILS_H diff --git a/source/source_lcao/module_operator_lcao/operator_force_stress_utils.hpp b/source/source_lcao/module_operator_lcao/operator_fs_utils.hpp similarity index 97% rename from source/source_lcao/module_operator_lcao/operator_force_stress_utils.hpp rename to source/source_lcao/module_operator_lcao/operator_fs_utils.hpp index 7a97a57296..7720903fa0 100644 --- a/source/source_lcao/module_operator_lcao/operator_force_stress_utils.hpp +++ b/source/source_lcao/module_operator_lcao/operator_fs_utils.hpp @@ -1,7 +1,7 @@ -#ifndef OPERATOR_FORCE_STRESS_UTILS_HPP -#define OPERATOR_FORCE_STRESS_UTILS_HPP +#ifndef OPERATOR_FS_UTILS_HPP +#define OPERATOR_FS_UTILS_HPP -#include "operator_force_stress_utils.h" +#include "operator_fs_utils.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" @@ -193,4 +193,4 @@ void cal_force_stress_2center( } // namespace OperatorForceStress -#endif // OPERATOR_FORCE_STRESS_UTILS_HPP +#endif // OPERATOR_FS_UTILS_HPP diff --git a/source/source_lcao/module_operator_lcao/operator_lcao.cpp b/source/source_lcao/module_operator_lcao/operator_lcao.cpp index 78463b3fa1..f4cea3504c 100644 --- a/source/source_lcao/module_operator_lcao/operator_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/operator_lcao.cpp @@ -3,7 +3,6 @@ #include "source_base/timer.h" #include "source_base/tool_title.h" #include "source_hamilt/module_hcontainer/hcontainer_funcs.h" -#include "source_hsolver/hsolver_lcao.h" #include "source_io/module_parameter/parameter.h" diff --git a/source/source_lcao/module_operator_lcao/operator_lcao.h b/source/source_lcao/module_operator_lcao/operator_lcao.h index 5c9d3616c3..e35d2f2f4f 100644 --- a/source/source_lcao/module_operator_lcao/operator_lcao.h +++ b/source/source_lcao/module_operator_lcao/operator_lcao.h @@ -4,7 +4,7 @@ #include "source_hamilt/matrixblock.h" #include "source_hamilt/operator.h" #include "source_hamilt/module_hcontainer/hcontainer.h" -#include "source_lcao/hs_matrix_k.hpp" +#include "source_hamilt/hs_matrix_k.h" namespace hamilt { diff --git a/source/source_lcao/module_operator_lcao/overlap.cpp b/source/source_lcao/module_operator_lcao/overlap.cpp index 86f1144e24..b9e4d8d1b0 100644 --- a/source/source_lcao/module_operator_lcao/overlap.cpp +++ b/source/source_lcao/module_operator_lcao/overlap.cpp @@ -457,9 +457,6 @@ void hamilt::Overlap>::output_SR_async_csr(const in ModuleBase::timer::end("OverlapNew", "output_SR_async_csr"); } -// Include force/stress implementation -#include "overlap_force_stress.hpp" - template class hamilt::Overlap>; template class hamilt::Overlap, double>>; template class hamilt::Overlap, std::complex>>; diff --git a/source/source_lcao/module_operator_lcao/overlap.h b/source/source_lcao/module_operator_lcao/overlap.h index 88b65483d8..c9df271c2b 100644 --- a/source/source_lcao/module_operator_lcao/overlap.h +++ b/source/source_lcao/module_operator_lcao/overlap.h @@ -1,5 +1,5 @@ -#ifndef W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_HAMILT_LCAO_HAMILT_LCAODFT_OPERATOR_LCAO_OVERLAP_H -#define W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_HAMILT_LCAO_HAMILT_LCAODFT_OPERATOR_LCAO_OVERLAP_H +#ifndef OVERLAP_H +#define OVERLAP_H #include "source_basis/module_ao/parallel_orbitals.h" #include "source_basis/module_nao/two_center_integrator.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" diff --git a/source/source_lcao/module_operator_lcao/overlap_force_stress.hpp b/source/source_lcao/module_operator_lcao/overlap_force_stress.hpp deleted file mode 100644 index 642ea7cccf..0000000000 --- a/source/source_lcao/module_operator_lcao/overlap_force_stress.hpp +++ /dev/null @@ -1,67 +0,0 @@ -#pragma once -#include "overlap.h" -#include "operator_force_stress_utils.hpp" -#include "source_base/timer.h" - -namespace hamilt -{ - -template -void Overlap>::cal_force_stress(const bool cal_force, - const bool cal_stress, - const HContainer* dmR, - ModuleBase::matrix& force, - ModuleBase::matrix& stress) -{ - ModuleBase::TITLE("Overlap", "cal_force_stress"); - ModuleBase::timer::start("Overlap", "cal_force_stress"); - - // Lambda function to calculate overlap integral and its gradient - auto integral_calc = [this](int T1, int L1, int N1, int M1, - int T2, int L2, int N2, int M2, - const ModuleBase::Vector3& dtau, - double* olm) { - this->intor_->calculate(T1, L1, N1, M1, T2, L2, N2, M2, - dtau * this->ucell->lat0, &olm[0], &olm[1]); - }; - - // Use unified template with ForceSign=-1, StressSign=+1 for overlap operator - OperatorForceStress::cal_force_stress_2center( - cal_force, cal_stress, dmR, this->ucell, this->gridD, - this->orb_cutoff_, dmR->get_paraV(), integral_calc, force, stress); - - ModuleBase::timer::end("Overlap", "cal_force_stress"); -} - -// Dummy implementations for cal_force_IJR and cal_stress_IJR -// These are not used in the simplified approach above -template -void Overlap>::cal_force_IJR( - const int& iat1, - const int& iat2, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const hamilt::BaseMatrix* dmR_pointer, - double* force1, - double* force2) -{ - // Not used in current implementation -} - -template -void Overlap>::cal_stress_IJR( - const int& iat1, - const int& iat2, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const hamilt::BaseMatrix* dmR_pointer, - const ModuleBase::Vector3& dis1, - const ModuleBase::Vector3& dis2, - double* stress) -{ - // Not used in current implementation -} - -} // namespace hamilt diff --git a/source/source_lcao/module_operator_lcao/overlap_fs.cpp b/source/source_lcao/module_operator_lcao/overlap_fs.cpp new file mode 100644 index 0000000000..76472f6fd8 --- /dev/null +++ b/source/source_lcao/module_operator_lcao/overlap_fs.cpp @@ -0,0 +1,130 @@ +#include "overlap.h" +#include "operator_fs_utils.hpp" +#include "source_base/timer.h" + +namespace hamilt +{ + +template +void Overlap>::cal_force_stress(const bool cal_force, + const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, + ModuleBase::matrix& stress) +{ + ModuleBase::TITLE("Overlap", "cal_force_stress"); + ModuleBase::timer::start("Overlap", "cal_force_stress"); + + // Lambda function to calculate overlap integral and its gradient + auto integral_calc = [this](int T1, int L1, int N1, int M1, + int T2, int L2, int N2, int M2, + const ModuleBase::Vector3& dtau, + double* olm) { + this->intor_->calculate(T1, L1, N1, M1, T2, L2, N2, M2, + dtau * this->ucell->lat0, &olm[0], &olm[1]); + }; + + // Use unified template with ForceSign=-1, StressSign=+1 for overlap operator + OperatorForceStress::cal_force_stress_2center( + cal_force, cal_stress, dmR, this->ucell, this->gridD, + this->orb_cutoff_, dmR->get_paraV(), integral_calc, force, stress); + + ModuleBase::timer::end("Overlap", "cal_force_stress"); +} + +// Dummy implementations for cal_force_IJR and cal_stress_IJR +// These are not used in the simplified approach above +template +void Overlap>::cal_force_IJR( + const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + double* force1, + double* force2) +{ + // Not used in current implementation +} + +template +void Overlap>::cal_stress_IJR( + const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress) +{ + // Not used in current implementation +} + +// explicit member function instantiations +template void Overlap>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); +template void Overlap, double>>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); +template void Overlap, std::complex>>::cal_force_stress( + const bool cal_force, const bool cal_stress, + const HContainer* dmR, + ModuleBase::matrix& force, ModuleBase::matrix& stress); + +template void Overlap>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + double* force1, double* force2); +template void Overlap, double>>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + double* force1, double* force2); +template void Overlap, std::complex>>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix>* dmR_pointer, + double* force1, double* force2); + +template void Overlap>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); +template void Overlap, double>>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix* dmR_pointer, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); +template void Overlap, std::complex>>::cal_stress_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const hamilt::BaseMatrix>* dmR_pointer, + const ModuleBase::Vector3& dis1, + const ModuleBase::Vector3& dis2, + double* stress); + +} // namespace hamilt diff --git a/source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp b/source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp index 6a9fa10c2e..7f0703f675 100644 --- a/source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp +++ b/source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp @@ -293,7 +293,6 @@ template void hamilt::TD_pot_hybrid>::contributeHk(int ik) { return; } -#include "td_pot_hybrid_force.hpp" template class hamilt::TD_pot_hybrid>; template class hamilt::TD_pot_hybrid, double>>; template class hamilt::TD_pot_hybrid, std::complex>>; diff --git a/source/source_lcao/module_operator_lcao/td_pot_hybrid_force.hpp b/source/source_lcao/module_operator_lcao/td_pot_hybrid_fs.cpp similarity index 84% rename from source/source_lcao/module_operator_lcao/td_pot_hybrid_force.hpp rename to source/source_lcao/module_operator_lcao/td_pot_hybrid_fs.cpp index 751dcb0e20..83423b5741 100644 --- a/source/source_lcao/module_operator_lcao/td_pot_hybrid_force.hpp +++ b/source/source_lcao/module_operator_lcao/td_pot_hybrid_fs.cpp @@ -1,8 +1,8 @@ -#pragma once #include "td_pot_hybrid.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/libm/libm.h" +#include "source_estate/module_pot/H_TDDFT_pw.h" namespace hamilt { @@ -176,4 +176,29 @@ void TD_pot_hybrid>::cal_force_IJR(const int& iat1, } } } +// explicit member function instantiations for cal_force_stress +template void TD_pot_hybrid>::cal_force_stress( + const bool cal_force, const HContainer* dmR, ModuleBase::matrix& force); +template void TD_pot_hybrid, double>>::cal_force_stress( + const bool cal_force, const HContainer* dmR, ModuleBase::matrix& force); +template void TD_pot_hybrid, std::complex>>::cal_force_stress( + const bool cal_force, const HContainer>* dmR, ModuleBase::matrix& force); + +// explicit member function instantiations for cal_force_IJR (generic template, TR=double cases) +// Note: ,complex> has a template<> specialization above +template void TD_pot_hybrid>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const ModuleBase::Vector3& dtau, + const ModuleBase::Vector3& dR, + double* dmR_pointer, + double* force1, double* force2); +template void TD_pot_hybrid, double>>::cal_force_IJR( + const int& iat1, const int& iat2, + const Parallel_Orbitals* paraV, + const ModuleBase::Vector3& dtau, + const ModuleBase::Vector3& dR, + double* dmR_pointer, + double* force1, double* force2); + }// namespace hamilt \ No newline at end of file diff --git a/source/source_lcao/module_operator_lcao/test/CMakeLists.txt b/source/source_lcao/module_operator_lcao/test/CMakeLists.txt index 9318815ad8..412d43b8b6 100644 --- a/source/source_lcao/module_operator_lcao/test/CMakeLists.txt +++ b/source/source_lcao/module_operator_lcao/test/CMakeLists.txt @@ -4,7 +4,7 @@ abacus_disable_feature_definitions(__FFT_TWO_CENTER) AddTest( TARGET MODULE_LCAO_operator_overlap_test LIBS parameter psi base device container - SOURCES test_overlap.cpp ../overlap.cpp ../operator_force_stress_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp + SOURCES test_overlap.cpp ../overlap.cpp ../overlap_fs.cpp ../operator_fs_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp ../../../source_hamilt/module_hcontainer/base_matrix.cpp ../../../source_hamilt/module_hcontainer/hcontainer.cpp ../../../source_hamilt/module_hcontainer/atom_pair.cpp ../../../source_hamilt/module_hcontainer/func_transfer.cpp ../../../source_hamilt/module_hcontainer/output_hcontainer.cpp ../../../source_hamilt/module_hcontainer/transfer.cpp ../../../source_io/module_output/sparse_matrix.cpp @@ -17,7 +17,7 @@ AddTest( AddTest( TARGET MODULE_LCAO_operator_overlap_serial_test LIBS parameter psi base device container - SOURCES test_overlap_serial.cpp ../overlap.cpp ../operator_force_stress_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp + SOURCES test_overlap_serial.cpp ../overlap.cpp ../overlap_fs.cpp ../operator_fs_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp ../../../source_hamilt/module_hcontainer/base_matrix.cpp ../../../source_hamilt/module_hcontainer/hcontainer.cpp ../../../source_hamilt/module_hcontainer/atom_pair.cpp ../../../source_hamilt/module_hcontainer/func_transfer.cpp ../../../source_hamilt/module_hcontainer/output_hcontainer.cpp ../../../source_hamilt/module_hcontainer/transfer.cpp ../../../source_io/module_output/sparse_matrix.cpp @@ -30,7 +30,7 @@ AddTest( AddTest( TARGET MODULE_LCAO_operator_overlap_cd_test LIBS parameter psi base device container - SOURCES test_overlap_cd.cpp ../overlap.cpp ../operator_force_stress_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp + SOURCES test_overlap_cd.cpp ../overlap.cpp ../overlap_fs.cpp ../operator_fs_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp ../../../source_hamilt/module_hcontainer/base_matrix.cpp ../../../source_hamilt/module_hcontainer/hcontainer.cpp ../../../source_hamilt/module_hcontainer/atom_pair.cpp ../../../source_hamilt/module_hcontainer/func_transfer.cpp ../../../source_hamilt/module_hcontainer/output_hcontainer.cpp ../../../source_hamilt/module_hcontainer/transfer.cpp ../../../source_io/module_output/sparse_matrix.cpp @@ -43,7 +43,7 @@ AddTest( AddTest( TARGET MODULE_LCAO_operator_ekinetic_test LIBS parameter psi base device container - SOURCES test_ekinetic.cpp ../ekinetic.cpp ../operator_force_stress_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp + SOURCES test_ekinetic.cpp ../ekinetic.cpp ../ekinetic_fs.cpp ../ekinetic_dh.cpp ../operator_fs_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp ../../../source_hamilt/module_hcontainer/base_matrix.cpp ../../../source_hamilt/module_hcontainer/hcontainer.cpp ../../../source_hamilt/module_hcontainer/atom_pair.cpp ../../../source_basis/module_ao/parallel_orbitals.cpp ../../../source_basis/module_ao/ORB_atomic_lm.cpp @@ -53,7 +53,7 @@ AddTest( AddTest( TARGET MODULE_LCAO_operator_ekinetic_serial_test LIBS parameter psi base device container - SOURCES test_ekinetic_serial.cpp ../ekinetic.cpp ../operator_force_stress_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp + SOURCES test_ekinetic_serial.cpp ../ekinetic.cpp ../ekinetic_fs.cpp ../ekinetic_dh.cpp ../operator_fs_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp ../../../source_hamilt/module_hcontainer/base_matrix.cpp ../../../source_hamilt/module_hcontainer/hcontainer.cpp ../../../source_hamilt/module_hcontainer/atom_pair.cpp ../../../source_basis/module_ao/parallel_orbitals.cpp ../../../source_basis/module_ao/ORB_atomic_lm.cpp @@ -63,7 +63,7 @@ AddTest( AddTest( TARGET MODULE_LCAO_operator_nonlocal_test LIBS parameter psi base device container - SOURCES test_nonlocal.cpp ../nonlocal.cpp ../operator_force_stress_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp + SOURCES test_nonlocal.cpp ../nonlocal.cpp ../nonlocal_fs.cpp ../nonlocal_dh.cpp ../operator_fs_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp ../../../source_hamilt/module_hcontainer/base_matrix.cpp ../../../source_hamilt/module_hcontainer/hcontainer.cpp ../../../source_hamilt/module_hcontainer/atom_pair.cpp ../../../source_basis/module_ao/parallel_orbitals.cpp ../../../source_basis/module_ao/ORB_atomic_lm.cpp @@ -73,7 +73,7 @@ AddTest( AddTest( TARGET MODULE_LCAO_operator_T_NL_cd_test LIBS parameter psi base device container - SOURCES test_T_NL_cd.cpp ../nonlocal.cpp ../ekinetic.cpp ../operator_force_stress_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp + SOURCES test_T_NL_cd.cpp ../nonlocal.cpp ../nonlocal_fs.cpp ../nonlocal_dh.cpp ../ekinetic.cpp ../ekinetic_fs.cpp ../ekinetic_dh.cpp ../operator_fs_utils.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp ../../../source_hamilt/module_hcontainer/base_matrix.cpp ../../../source_hamilt/module_hcontainer/hcontainer.cpp ../../../source_hamilt/module_hcontainer/atom_pair.cpp ../../../source_basis/module_ao/parallel_orbitals.cpp ../../../source_basis/module_ao/ORB_atomic_lm.cpp @@ -83,7 +83,7 @@ AddTest( AddTest( TARGET MODULE_LCAO_operator_dftu_test LIBS parameter psi base device container - SOURCES test_dftu.cpp ../dftu_lcao.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp + SOURCES test_dftu.cpp ../dftu_lcao.cpp ../dftu_fs.cpp ../../../source_hamilt/module_hcontainer/func_folding.cpp ../../../source_hamilt/module_hcontainer/base_matrix.cpp ../../../source_hamilt/module_hcontainer/hcontainer.cpp ../../../source_hamilt/module_hcontainer/atom_pair.cpp ../../../source_basis/module_ao/parallel_orbitals.cpp ../../../source_basis/module_ao/ORB_atomic_lm.cpp diff --git a/source/source_lcao/module_operator_lcao/veff_dh.hpp b/source/source_lcao/module_operator_lcao/veff_dh.cpp similarity index 96% rename from source/source_lcao/module_operator_lcao/veff_dh.hpp rename to source/source_lcao/module_operator_lcao/veff_dh.cpp index b3afd26fda..046395570a 100644 --- a/source/source_lcao/module_operator_lcao/veff_dh.hpp +++ b/source/source_lcao/module_operator_lcao/veff_dh.cpp @@ -1,4 +1,3 @@ -#pragma once #include "source_base/timer.h" #include "source_estate/module_charge/charge.h" #include "source_estate/module_pot/H_Hartree_pw.h" @@ -524,4 +523,24 @@ void Veff>::cal_dH(std::array>::cal_dH( + std::array*>, 3>& dhR, + const std::string& hellmann_feynman_type, + const std::vector*>& dmR, + const Charge* chg, + const int ispin); +template void Veff, double>>::cal_dH( + std::array*>, 3>& dhR, + const std::string& hellmann_feynman_type, + const std::vector*>& dmR, + const Charge* chg, + const int ispin); +template void Veff, std::complex>>::cal_dH( + std::array*>, 3>& dhR, + const std::string& hellmann_feynman_type, + const std::vector*>& dmR, + const Charge* chg, + const int ispin); + } // namespace hamilt diff --git a/source/source_lcao/module_operator_lcao/veff_lcao.cpp b/source/source_lcao/module_operator_lcao/veff_lcao.cpp index 5e2f68e43d..3ee99660bc 100644 --- a/source/source_lcao/module_operator_lcao/veff_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/veff_lcao.cpp @@ -1,5 +1,4 @@ #include "veff_lcao.h" -#include "veff_dh.hpp" #include "source_base/timer.h" #include "source_io/module_parameter/parameter.h" #include "source_base/tool_title.h" diff --git a/source/source_lcao/module_rdmft/rdmft.h b/source/source_lcao/module_rdmft/rdmft.h index 38f328d655..e084568720 100644 --- a/source/source_lcao/module_rdmft/rdmft.h +++ b/source/source_lcao/module_rdmft/rdmft.h @@ -16,7 +16,7 @@ #include "source_lcao/module_operator_lcao/operator_lcao.h" #include "source_hamilt/module_hcontainer/hcontainer.h" -#include "source_lcao/hs_matrix_k.hpp" +#include "source_hamilt/hs_matrix_k.h" #ifdef __EXX // Exx_LRI forward declaration, full definition in Exx_LRI.h (moved to .cpp) diff --git a/source/source_lcao/module_rdmft/rdmft_tools.h b/source/source_lcao/module_rdmft/rdmft_tools.h index 4d8331b4d0..52b9817043 100644 --- a/source/source_lcao/module_rdmft/rdmft_tools.h +++ b/source/source_lcao/module_rdmft/rdmft_tools.h @@ -19,7 +19,7 @@ #include "source_estate/module_dm/density_matrix.h" #include "source_hamilt/module_hcontainer/hcontainer.h" -#include "source_lcao/hs_matrix_k.hpp" +#include "source_hamilt/hs_matrix_k.h" #include "source_lcao/module_operator_lcao/operator_lcao.h" diff --git a/source/source_lcao/module_ri/Exx_LRI_interface.hpp b/source/source_lcao/module_ri/Exx_LRI_interface.hpp index b25cbc34ed..a99e2c4fc9 100644 --- a/source/source_lcao/module_ri/Exx_LRI_interface.hpp +++ b/source/source_lcao/module_ri/Exx_LRI_interface.hpp @@ -9,6 +9,7 @@ #include "source_io/module_output/csr_reader.h" #include "source_io/module_parameter/parameter.h" #include "source_io/module_restart/restart.h" +#include "source_io/module_restart/restart_exx_csr.h" #include "source_lcao/module_operator_lcao/op_exx_lcao.h" #include "source_lcao/module_ri/exx_abfs-jle.h" diff --git a/source/source_psi/psi_prepare.cpp b/source/source_psi/psi_prepare.cpp index adc10eeb3f..eada5f13d1 100644 --- a/source/source_psi/psi_prepare.cpp +++ b/source/source_psi/psi_prepare.cpp @@ -35,7 +35,7 @@ PSIPrepare::PSIPrepare(const std::string& init_wfc_in, } template -void PSIPrepare::prepare_init(const int& random_seed) +void PSIPrepare::prepare_init(const int& random_seed, const int istep) { // under restriction of C++11, std::unique_ptr can not be allocate via std::make_unique @@ -54,17 +54,23 @@ void PSIPrepare::prepare_init(const int& random_seed) } else if ((this->init_wfc.substr(0, 6) == "atomic") && (this->ucell.natomwfc == 0)) { - std::cout << " WARNING: init_wfc = " + this->init_wfc + - " requires atomic pseudo wavefunctions(PP_PSWFC),\n but none available." - " Automatically switch to random initialization." << std::endl; + // The switch to random initialization still happens every ion step, + // but the warning is printed only on the first step to avoid + // spamming relax/cell-relax output with the same message. + if (istep == 0) + { + std::cout << " WARNING: init_wfc = " + this->init_wfc + + " requires atomic pseudo wavefunctions(PP_PSWFC),\n but none available." + " Automatically switch to random initialization." << std::endl; + GlobalV::ofs_running << "\n WARNING:\n init_wfc = " + this->init_wfc + " requires atomic pseudo wavefunctions(PP_PSWFC), but none available. \n" + " Automatically switch to random initialization.\n" + " Note: Random starting wavefunctions may slow down convergence.\n" + " For faster convergence, consider using:\n" + " 1) A pseudopotential file that includes atomic wavefunctions (with PP_PSWFC), or\n" + " 2) Numerical atomic orbitals with 'init_wfc = nao' or 'nao+random' if available.\n" + << std::endl; + } GlobalV::ofs_running << "\n Using RANDOM starting wave functions for all " << PARAM.inp.nbands << " bands\n"; - GlobalV::ofs_running << "\n WARNING:\n init_wfc = " + this->init_wfc + " requires atomic pseudo wavefunctions(PP_PSWFC), but none available. \n" - " Automatically switch to random initialization.\n" - " Note: Random starting wavefunctions may slow down convergence.\n" - " For faster convergence, consider using:\n" - " 1) A pseudopotential file that includes atomic wavefunctions (with PP_PSWFC), or\n" - " 2) Numerical atomic orbitals with 'init_wfc = nao' or 'nao+random' if available.\n" - << std::endl; this->psi_initer = std::unique_ptr>(new psi_init_random()); } else if (this->init_wfc == "atomic" diff --git a/source/source_psi/psi_prepare.h b/source/source_psi/psi_prepare.h index 4b35f54521..20ec2eeb1b 100644 --- a/source/source_psi/psi_prepare.h +++ b/source/source_psi/psi_prepare.h @@ -24,7 +24,10 @@ class PSIPrepare : public PSIPrepareBase ~PSIPrepare(){}; ///@brief prepare the wavefunction initialization - void prepare_init(const int& random_seed); + ///@param random_seed seed for random initialization + ///@param istep current ion/relax step; informational warnings are only + /// printed on the first step to avoid spamming relax output + void prepare_init(const int& random_seed, const int istep); //------------------------ only for psi_initializer -------------------- /** diff --git a/source/source_psi/psi_prepare_base.h b/source/source_psi/psi_prepare_base.h index c7a10718bf..1feb0d8d68 100644 --- a/source/source_psi/psi_prepare_base.h +++ b/source/source_psi/psi_prepare_base.h @@ -16,7 +16,11 @@ class PSIPrepareBase public: PSIPrepareBase() = default; virtual ~PSIPrepareBase() = default; - virtual void prepare_init(const int& random_seed) = 0; + ///@brief prepare the wavefunction initialization + ///@param random_seed seed for random initialization + ///@param istep current ion/relax step; used to suppress repeated + /// informational warnings (only printed on the first step) + virtual void prepare_init(const int& random_seed, const int istep) = 0; }; } // namespace psi diff --git a/source/source_psi/setup_psi_pw.cpp b/source/source_psi/setup_psi_pw.cpp index f5bc240292..816842d444 100644 --- a/source/source_psi/setup_psi_pw.cpp +++ b/source/source_psi/setup_psi_pw.cpp @@ -21,7 +21,10 @@ void Setup_Psi_pw::before_runner_impl( allocate_psi(this->psi_cpu, kv.get_nks(), kv.ngk, PARAM.globalv.nbands_l, pw_wfc.npwk_max); auto* p_psi_init = static_cast*>(this->p_psi_init); - p_psi_init->prepare_init(inp.pw_seed); + // before_runner_impl is invoked only once before the ion dynamics starts, + // so istep == 0 here ensures any one-time informational warnings are + // printed from this initial setup call. + p_psi_init->prepare_init(inp.pw_seed, 0); if (std::is_same::value) { precision_type_ = PrecisionType::Float; diff --git a/source/source_pw/module_pwdft/CMakeLists.txt b/source/source_pw/module_pwdft/CMakeLists.txt index dec6cd86d8..a45bdb70e2 100644 --- a/source/source_pw/module_pwdft/CMakeLists.txt +++ b/source/source_pw/module_pwdft/CMakeLists.txt @@ -16,7 +16,7 @@ list(APPEND objects setup_pwrho.cpp setup_pwwfc.cpp update_cell_pw.cpp - dftu_pw.cpp + setup_dftu_pw.cpp deltaspin_pw.cpp forces_nl.cpp forces_cc.cpp diff --git a/source/source_pw/module_pwdft/dftu_pw.cpp b/source/source_pw/module_pwdft/setup_dftu_pw.cpp similarity index 94% rename from source/source_pw/module_pwdft/dftu_pw.cpp rename to source/source_pw/module_pwdft/setup_dftu_pw.cpp index 667612e23b..42e088c11a 100644 --- a/source/source_pw/module_pwdft/dftu_pw.cpp +++ b/source/source_pw/module_pwdft/setup_dftu_pw.cpp @@ -1,4 +1,4 @@ -#include "source_pw/module_pwdft/dftu_pw.h" +#include "source_pw/module_pwdft/setup_dftu_pw.h" #include "source_lcao/module_dftu/dftu.h" #include "source_io/module_parameter/parameter.h" diff --git a/source/source_pw/module_pwdft/dftu_pw.h b/source/source_pw/module_pwdft/setup_dftu_pw.h similarity index 91% rename from source/source_pw/module_pwdft/dftu_pw.h rename to source/source_pw/module_pwdft/setup_dftu_pw.h index 94e24f31ff..1235d15928 100644 --- a/source/source_pw/module_pwdft/dftu_pw.h +++ b/source/source_pw/module_pwdft/setup_dftu_pw.h @@ -1,5 +1,5 @@ -#ifndef DFTU_PW_H -#define DFTU_PW_H +#ifndef SETUP_DFTU_PW_H +#define SETUP_DFTU_PW_H #include "source_cell/unitcell.h" #include "source_base/matrix.h"