diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp index b3431c0f5d8..aafeb18e6aa 100644 --- a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp +++ b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp @@ -1,6 +1,7 @@ #include "ctrl_scf_lcao.h" // use ctrl_scf_lcao() #include "source_base/formatter.h" +#include "source_base/tool_quit.h" // use ModuleBase::WARNING_QUIT #include "source_estate/elecstate_lcao.h" // use elecstate::ElecState #include "source_hamilt/hamilt.h" // use Hamilt #include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO @@ -55,7 +56,15 @@ void setup_exx_dh_params(ModuleIO::WriteDHParams& dh_params, Exx_NAO void setup_exx_h_params(ModuleIO::WriteHParams& h_params, Exx_NAO& exx_nao) -{} +{ + // Only the gamma-only (TK==double) specialization below actually writes V^EXX(R). + // This generic body is instantiated for the multi-k (TK==std::complex) path, where the + // EXX-H output is unsupported. Reject it explicitly here so the request cannot be silently + // dropped (the WARNING_QUIT inside write_h_exx is unreachable at multi-k). + ModuleBase::WARNING_QUIT("setup_exx_h_params", + "out_mat_h_exx is only supported for gamma-only: the V^EXX(R) " + "output is not available at multi-k. Use gamma_only."); +} template <> void setup_exx_h_params(ModuleIO::WriteHParams& h_params, Exx_NAO& exx_nao) diff --git a/source/source_io/module_dhs/write_dH.cpp b/source/source_io/module_dhs/write_dH.cpp index 37a13058552..ea062bf45a4 100644 --- a/source/source_io/module_dhs/write_dH.cpp +++ b/source/source_io/module_dhs/write_dH.cpp @@ -8,6 +8,9 @@ #include "source_io/module_parameter/parameter.h" #include "source_lcao/module_hcontainer/hcontainer_funcs.h" #include "source_lcao/module_hcontainer/output_hcontainer.h" +#ifdef __EXX +#include "source_hamilt/module_xc/exx_info.h" +#endif #include #include @@ -138,6 +141,20 @@ void write_dH_components(WriteDHParams& params) "nspin=4 (noncollinear) yet; only nspin=1 and nspin=2."); } +#ifdef __EXX + // The EXX interfaces carried by WriteDHParams are gamma-only (see write_dH.h): at multi-k + // dH^EXX would be the derivative with respect to every mirror atom, which this output is + // not meant for. Quit instead of writing a dH sum that silently omits the EXX term. + if (GlobalC::exx_info.info_global.cal_exx && !PARAM.globalv.gamma_only_local + && (PARAM.inp.out_mat_dh[0] || PARAM.inp.out_mat_dh_exx[0])) + { + ModuleBase::WARNING_QUIT("write_dH_components", + "out_mat_dh (the dH sum) and out_mat_dh_exx are only supported for gamma-only when EXX is on. " + "Use gamma_only, or request the individual non-EXX terms (out_mat_dh_t, " + "out_mat_dh_vnl, out_mat_dh_vl, out_mat_dh_vh, out_mat_dh_vxc)."); + } +#endif + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; GlobalV::ofs_running << " | |" << std::endl; GlobalV::ofs_running << " | #Print out dH/dR components# |" << std::endl; diff --git a/source/source_io/module_dhs/write_dH.h b/source/source_io/module_dhs/write_dH.h index 870a19f7e30..0fc944f426e 100644 --- a/source/source_io/module_dhs/write_dH.h +++ b/source/source_io/module_dhs/write_dH.h @@ -54,8 +54,15 @@ struct WriteDHParams std::vector*> dmR; const Charge* chg = nullptr; // ground-state charge for XC Hellmann-Feynman (FDM) #ifdef __EXX - // gamma (TK==double) exx interfaces used by write_dH_exx; exactly one is set depending on - // GlobalC::exx_info.info_ri.real_number (exd: real Hexx, exc: complex Hexx). + // The gamma-only (TK==double) exx interfaces used by write_dH_exx. + // Deliberately NOT templated on TK, for two reasons: + // 1. Physics: at multi-k the derivative would be taken with respect to every mirror + // atom of the periodic images, which is not what this output is used for. + // 2. Cost: templating these pointers on TK would force WriteHParams, WriteDHParams and + // every free function taking them to become templates as well -- a large, purely + // mechanical change for a case nobody needs. + // Multi-k + EXX is therefore rejected up front (see write_dH_components) + // instead of silently producing output with the EXX term missing. Exx_LRI_Interface* exd = nullptr; Exx_LRI_Interface>* exc = 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 ae6cd964ebd..a94e234a9a2 100644 --- a/source/source_io/module_hs/write_H_terms.cpp +++ b/source/source_io/module_hs/write_H_terms.cpp @@ -411,6 +411,8 @@ void write_h_exx(WriteHParams& params) ModuleBase::TITLE("ModuleIO", "write_h_exx"); ModuleBase::timer::start("ModuleIO", "write_h_exx"); + // Multi-k out_mat_h_exx is rejected upstream at the call site (setup_exx_h_params in + // ctrl_scf_lcao.cpp); this function is only reached on the gamma-only path. const UnitCell& ucell = *params.ucell; const Parallel_Orbitals& pv = *params.pv; const K_Vectors& kv = *params.kv; diff --git a/source/source_io/module_hs/write_H_terms.h b/source/source_io/module_hs/write_H_terms.h index 05a88c70a63..7dbede427a8 100644 --- a/source/source_io/module_hs/write_H_terms.h +++ b/source/source_io/module_hs/write_H_terms.h @@ -38,8 +38,12 @@ struct WriteHParams int nat = 0; bool also_hR = false; // H(k) is always written; H(R) (CSR) only when this is true #ifdef __EXX - // gamma (TK==double) exx interfaces used by write_h_exx; exactly one is set depending on - // GlobalC::exx_info.info_ri.real_number (exd: real Hexx, exc: complex Hexx). + // The gamma-only (TK==double) exx interfaces used by write_h_exx. + // Deliberately NOT templated on TK, because it would force WriteHParams, WriteDHParams and + // every free function taking them to become templates as well -- a large, purely + // mechanical change for a case nobody needs. + // Multi-k + EXX is therefore rejected up front (see write_h_exx) + // instead of silently producing output with the EXX term missing. Exx_LRI_Interface* exd = nullptr; Exx_LRI_Interface>* exc = nullptr; #endif diff --git a/source/source_lcao/module_operator_lcao/nonlocal_dh.hpp.bak b/source/source_lcao/module_operator_lcao/nonlocal_dh.hpp.bak deleted file mode 100644 index d247840af0b..00000000000 --- a/source/source_lcao/module_operator_lcao/nonlocal_dh.hpp.bak +++ /dev/null @@ -1,256 +0,0 @@ -#pragma once -#include "nonlocal.h" -#include "operator_force_stress_utils.h" -#include "source_base/timer.h" - -namespace hamilt -{ - -template -void Nonlocal>::cal_dH(std::array*>, 3>& dhR) -{ - ModuleBase::TITLE("Nonlocal", "cal_dH"); - ModuleBase::timer::start("Nonlocal", "cal_dH"); - - const int nat = this->ucell->nat; - assert(static_cast(dhR[0].size()) == nat); - const Parallel_Orbitals* paraV = dhR[0][0]->get_paraV(); - const int npol = this->ucell->get_npol(); - - for (int iat0 = 0; iat0 < nat; iat0++) - { - auto tau0 = this->ucell->get_tau(iat0); - int I0 = 0, T0 = 0; - this->ucell->iat2iait(iat0, &I0, &T0); - - AdjacentAtomInfo adjs; - this->gridD->Find_atom(*this->ucell, tau0, T0, I0, &adjs); - - std::vector is_adj(adjs.adj_num + 1, false); - for (int ad = 0; ad < adjs.adj_num + 1; ++ad) - { - const int T1 = adjs.ntype[ad]; - const int I1 = adjs.natom[ad]; - const int iat1 = this->ucell->itia2iat(T1, I1); - const ModuleBase::Vector3& R_index1 = adjs.box[ad]; - if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0 - < this->orb_cutoff_[T1] + this->ucell->infoNL.Beta[T0].get_rcut_max()) - { - is_adj[ad] = true; - } - } - - for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) - { - if (!is_adj[ad1]) - continue; - const int T1 = adjs.ntype[ad1]; - const int I1 = adjs.natom[ad1]; - const int iat1 = this->ucell->itia2iat(T1, I1); - const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; - - for (int ad2 = 0; ad2 < adjs.adj_num + 1; ++ad2) - { - if (!is_adj[ad2]) - continue; - const int T2 = adjs.ntype[ad2]; - const int I2 = adjs.natom[ad2]; - const int iat2 = this->ucell->itia2iat(T2, I2); - const ModuleBase::Vector3& R_index2 = adjs.box[ad2]; - - if (paraV->get_row_size(iat1) <= 0 || paraV->get_col_size(iat2) <= 0) - { - continue; - } - - ModuleBase::Vector3 dR(R_index2.x - R_index1.x, R_index2.y - R_index1.y, R_index2.z - R_index1.z); - - hamilt::AtomPair ap(iat1, iat2, dR.x, dR.y, dR.z, paraV); - for (int iat = 0; iat < nat; ++iat) - { - for (int d = 0; d < 3; ++d) - dhR[d][iat]->insert_pair(ap); - } - } - } - } - - for (int iat = 0; iat < nat; ++iat) - { - for (int d = 0; d < 3; ++d) - dhR[d][iat]->allocate(nullptr, true); - } - -#pragma omp parallel - { -#pragma omp for schedule(dynamic) - for (int iat0 = 0; iat0 < nat; iat0++) - { - auto tau0 = this->ucell->get_tau(iat0); - int I0 = 0, T0 = 0; - this->ucell->iat2iait(iat0, &I0, &T0); - - AdjacentAtomInfo adjs; - this->gridD->Find_atom(*this->ucell, tau0, T0, I0, &adjs); - - std::vector is_adj(adjs.adj_num + 1, false); - for (int ad = 0; ad < adjs.adj_num + 1; ++ad) - { - const int T1 = adjs.ntype[ad]; - const int I1 = adjs.natom[ad]; - const int iat1 = this->ucell->itia2iat(T1, I1); - const ModuleBase::Vector3& R_index1 = adjs.box[ad]; - if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0 - < this->orb_cutoff_[T1] + this->ucell->infoNL.Beta[T0].get_rcut_max()) - { - is_adj[ad] = true; - } - } - - std::vector>> nlm_iat0(adjs.adj_num + 1); - - for (int ad = 0; ad < adjs.adj_num + 1; ++ad) - { - if (!is_adj[ad]) - continue; - - const int T1 = adjs.ntype[ad]; - const int I1 = adjs.natom[ad]; - const int iat1 = this->ucell->itia2iat(T1, I1); - const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad]; - const Atom* atom1 = &this->ucell->atoms[T1]; - - auto all_indexes = paraV->get_indexes_row(iat1); - auto col_indexes = paraV->get_indexes_col(iat1); - all_indexes.insert(all_indexes.end(), col_indexes.begin(), col_indexes.end()); - std::sort(all_indexes.begin(), all_indexes.end()); - all_indexes.erase(std::unique(all_indexes.begin(), all_indexes.end()), all_indexes.end()); - - for (size_t iw1l = 0; iw1l < all_indexes.size(); iw1l += npol) - { - const int iw1 = all_indexes[iw1l] / npol; - std::vector> nlm; - - OperatorForceStress::OrbitalQuantumNumbers qn1 = OperatorForceStress::get_orbital_qn(*atom1, iw1); - - ModuleBase::Vector3 dtau_at = tau0 - tau1; - this->intor_->snap(T1, qn1.L, qn1.N, qn1.M, T0, dtau_at * this->ucell->lat0, true, nlm); - - const size_t length = nlm[0].size(); - std::vector nlm_target(length * 4); - for (size_t index = 0; index < length; index++) - { - for (int n = 0; n < 4; n++) - nlm_target[index + n * length] = nlm[n][index]; - } - nlm_iat0[ad].insert({all_indexes[iw1l], nlm_target}); - } - } - - for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) - { - if (!is_adj[ad1]) - continue; - const int T1 = adjs.ntype[ad1]; - const int I1 = adjs.natom[ad1]; - const int iat1 = this->ucell->itia2iat(T1, I1); - const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; - - for (int ad2 = 0; ad2 < adjs.adj_num + 1; ++ad2) - { - if (!is_adj[ad2]) - continue; - const int T2 = adjs.ntype[ad2]; - const int I2 = adjs.natom[ad2]; - const int iat2 = this->ucell->itia2iat(T2, I2); - const ModuleBase::Vector3& R_index2 = adjs.box[ad2]; - - ModuleBase::Vector3 dR(R_index2.x - R_index1.x, - R_index2.y - R_index1.y, - R_index2.z - R_index1.z); - - // destination block (iat1,iat2,dR) for the three differentiated atoms: - // iat1 (orbital 1), iat2 (orbital 2), iat0 (projector / Hellmann-Feynman) - hamilt::BaseMatrix* m1[3]; - hamilt::BaseMatrix* m2[3]; - hamilt::BaseMatrix* m0[3]; - for (int d = 0; d < 3; ++d) - { - m1[d] = dhR[d][iat1]->find_matrix(iat1, iat2, dR.x, dR.y, dR.z); - m2[d] = dhR[d][iat2]->find_matrix(iat1, iat2, dR.x, dR.y, dR.z); - m0[d] = dhR[d][iat0]->find_matrix(iat1, iat2, dR.x, dR.y, dR.z); - } - - if (!m1[0] || !m1[1] || !m1[2] || !m2[0] || !m2[1] || !m2[2] || !m0[0] || !m0[1] || !m0[2]) - continue; - - double* p1[3] = {m1[0]->get_pointer(), m1[1]->get_pointer(), m1[2]->get_pointer()}; - double* p2[3] = {m2[0]->get_pointer(), m2[1]->get_pointer(), m2[2]->get_pointer()}; - double* p0[3] = {m0[0]->get_pointer(), m0[1]->get_pointer(), m0[2]->get_pointer()}; - const int col_sz = m1[0]->get_col_size(); - - auto& nlm1_all = nlm_iat0[ad1]; - auto& nlm2_all = nlm_iat0[ad2]; - - auto row_indexes = paraV->get_indexes_row(iat1); - auto col_indexes = paraV->get_indexes_col(iat2); - - for (size_t iw1l = 0; iw1l < row_indexes.size(); iw1l++) - { - auto it1 = nlm1_all.find(row_indexes[iw1l]); - if (it1 == nlm1_all.end()) - continue; - const std::vector& nlm1 = it1->second; - const size_t length = nlm1.size() / 4; - const int iw1_row = paraV->global2local_row(row_indexes[iw1l]); - - for (size_t iw2l = 0; iw2l < col_indexes.size(); iw2l++) - { - auto it2 = nlm2_all.find(col_indexes[iw2l]); - if (it2 == nlm2_all.end()) - continue; - const std::vector& nlm2 = it2->second; - const int iw2_col = paraV->global2local_col(col_indexes[iw2l]); - - // tU = D (orbital 1 moves) - // tV = D (orbital 2 moves) - double tU[3] = {0, 0, 0}; - double tV[3] = {0, 0, 0}; - - for (int no = 0; no < this->ucell->atoms[T0].ncpp.non_zero_count_soc[0]; no++) - { - const int p1_idx = this->ucell->atoms[T0].ncpp.index1_soc[0][no]; - const int p2_idx = this->ucell->atoms[T0].ncpp.index2_soc[0][no]; - const double* tmp_d = nullptr; - this->ucell->atoms[T0].ncpp.get_d(0, p1_idx, p2_idx, tmp_d); - for (int d = 0; d < 3; ++d) - { - tU[d] += nlm1[p1_idx + length * (d + 1)] * nlm2[p2_idx] * (*tmp_d); - tV[d] += nlm1[p1_idx] * nlm2[p2_idx + length * (d + 1)] * (*tmp_d); - } - } - - const int idx = iw1_row * col_sz + iw2_col; - // d/dtau_iat1, d/dtau_iat2, and (translational invariance) d/dtau_iat0 - // dtau=-- - // =- for Hellmann-Feynman terms - for (int d = 0; d < 3; ++d) - { -#pragma omp atomic - p1[d][idx] -= tU[d]; -#pragma omp atomic - p2[d][idx] -= tV[d]; -#pragma omp atomic - p0[d][idx] += tU[d] + tV[d]; - } - } - } - } - } - } - } - - ModuleBase::timer::end("Nonlocal", "cal_dH"); -} - -} // namespace hamilt