small fixes

Author: maki49

source/source_esolver/esolver_lr_lcao_tddft.cpp

@@ -28,6 +28,42 @@
 // gradient
 #include "source_lcao/module_lr/Grad/multipliers/zeq_solver.h"
 
+#ifdef __EXX
+namespace
+{
+    /// Screening of the Coulomb operator that `Exx_LRI` is built with: `hse` is erfc-screened,
+    /// `hf` and `pbe0` are bare.
+    ///
+    /// One `Exx_LRI` carries ONE screening, and it may be needed for two different reasons: the
+    /// LR kernel (when `xc_kernel` is a hybrid) and the ground-state force (when
+    /// `dft_functional` is a hybrid). Keying the choice off `xc_kernel` alone -- which is what
+    /// this used to do -- silently produced *unscreened* exchange whenever the object existed
+    /// only for the force, e.g. `dft_functional hse` with `xc_kernel lda` or `rpa`.
+    Conv_Coulomb_Pot_K::Ccp_Type exx_ccp_type(const std::string& name)
+    {
+        return (name == "hse") ? Conv_Coulomb_Pot_K::Ccp_Type::Erfc
+                               : Conv_Coulomb_Pot_K::Ccp_Type::Hf;
+    }
+
+    /// Which functional the single `Exx_LRI` must follow. Prefer the LR kernel, since that one
+    /// enters the eigenproblem; fall back to the ground-state functional, which is the only
+    /// reason the object exists when the kernel is local.
+    std::string exx_source(const std::string& xc_kernel, const std::string& dft_functional)
+    {
+        const bool k = LR::exx_kernel_list().count(xc_kernel) > 0;
+        const bool g = LR::exx_kernel_list().count(dft_functional) > 0;
+        if (k && g && xc_kernel != dft_functional)
+        {
+            GlobalV::ofs_running << " WARNING: xc_kernel (" << xc_kernel << ") and dft_functional ("
+                << dft_functional << ") are two DIFFERENT hybrids. A single Exx_LRI carries one"
+                " screening, so only " << xc_kernel << "'s is used; the ground-state EXX force"
+                " will be inconsistent." << std::endl;
+        }
+        return k ? xc_kernel : dft_functional;
+    }
+}
+#endif
+
 #ifdef __EXX
 template<>
 void ModuleESolver::ESolver_LR<double>::move_exx_lri(std::shared_ptr<Exx_LRI<double>>& exx_ks)
@@ -349,9 +385,7 @@ void ModuleESolver::ESolver_LR<T, TR>::initialize_from_ks_(ModuleESolver::ESolve
             this->move_exx_lri(ks_sol.exx_nao.exc->exx_ptr);
         } else    // construct C, V from scratch
         {
-            // set ccp_type according to the xc_kernel
-            if (xc_kernel == "hf" || xc_kernel == "pbe0") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; }
-            else if (xc_kernel == "hse") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Erfc; }
+            exx_info.info_global.ccp_type = exx_ccp_type(exx_source(xc_kernel, dft_functional));
             exx_info.sync_from_global();
             this->exx_lri = std::make_shared<Exx_LRI<T>>(exx_info.info_ri);
             this->exx_lri->init(MPI_COMM_WORLD, ucell,this->kv, ks_sol.orb_);   
@@ -495,9 +529,8 @@ void ModuleESolver::ESolver_LR<T, TR>::initialize_from_unitcell_(UnitCell& ucell
     if (((exx_kernel_list().count(xc_kernel)) && this->input.lr_solver != "spectrum")
         || (PARAM.inp.cal_force && (exx_kernel_list().count(PARAM.inp.dft_functional) )))
     {
-        // set ccp_type according to the xc_kernel
-        if (xc_kernel == "hf") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; }
-        else if (xc_kernel == "hse") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Erfc; }
+        exx_info.info_global.ccp_type =
+            exx_ccp_type(exx_source(xc_kernel, LR_Util::tolower(PARAM.inp.dft_functional)));
         this->exx_lri = std::make_shared<Exx_LRI<T>>(exx_info.info_ri);
         this->exx_lri->init(MPI_COMM_WORLD, ucell,this->kv, orb);
         this->exx_lri->cal_exx_ions(ucell,input.out_ri_cv);

source/source_lcao/module_lr/Grad/force/cal_hs_grad.h

@@ -1,4 +1,5 @@
 #pragma once
+#include <cassert>
 #include "source_cell/module_neighbor/sltk_grid_driver.h"
 #include "source_cell/unitcell.h"
 #include "source_basis/module_ao/parallel_orbitals.h"
@@ -55,14 +56,7 @@ hamilt::HContainer<TR> build_hcontainer_local_op(const UnitCell& ucell, const Gr
     return hcontainer;
 }
 
-// ModuleBase::Vector3<double> operator*(const ModuleBase::Vector3<double>& row_vec, ModuleBase::Matrix3& mat3)
-// {
-//     return ModuleBase::Vector3<double>(row_vec.x * mat3.e11 + row_vec.y * mat3.e21 + row_vec.z * mat3.e31,
-//                                        row_vec.x * mat3.e12 + row_vec.y * mat3.e22 + row_vec.z * mat3.e32,
-//                                        row_vec.x * mat3.e13 + row_vec.y * mat3.e23 + row_vec.z * mat3.e33);
-// }
-
-
+/// @brief  Calculate <dPhi|phi> or <dPhi|T|phi> by 2-center integration
 inline std::vector<hamilt::HContainer<double>> cal_hs_grad(const char job,
     const UnitCell& ucell,
     const Parallel_Orbitals& pv,
@@ -107,9 +101,17 @@ inline std::vector<hamilt::HContainer<double>> cal_hs_grad(const char job,
 
             for (int iR = 0;iR < nR;++iR)
             {
-                ModuleBase::Vector3<double> R(*it++, *it++, *it++); // int to double
+                // Read the three components in separate statements. The order in which function
+                // arguments are evaluated is UNSPECIFIED in C++, so `Vector3<double>(*it++, *it++, *it++)`
+                // may store the R triple permuted and silently corrupt periodic systems.
+                const int Rx = *it++;
+                const int Ry = *it++;
+                const int Rz = *it++;
+                const ModuleBase::Vector3<double> R(Rx, Ry, Rz);    // int to double
                 ModuleBase::Vector3<double> relative_position = (tau1 - tau0 + R * ucell.latvec) * ucell.lat0;
                 hamilt::BaseMatrix<double>* dHS_block = dHS[ixyz].find_matrix(iat0, iat1, R.x, R.y, R.z);
+                // `ijr_info` came from this very container, so a miss means the indices are wrong.
+                assert(dHS_block != nullptr);
 
                 // OMP can be used here
                 for (int lw0 = 0;lw0 < row_indexes.size();lw0 += npol)    // spin 1-3 of dHS is not needed at nspin=4

source/source_lcao/module_lr/Grad/force/lr_force.h

@@ -9,6 +9,11 @@ using TAC = std::pair<int, std::array<int, 3>>;
 #endif
 namespace LR
 {
+    /// `dm_gs` carries the ground-state occupations, at nspin=1 they are 2 for fully occupied
+    /// bands, so any term that contracts against ONE channel of `dm_gs` needs this factor.
+    /// Closed-shell bookkeeping only: the open-shell path sums the two channels explicitly and
+    /// must not apply it.
+    inline double gs_dm_channel_factor() { return (PARAM.inp.nspin == 1) ? 0.5 : 1.0; }
 
     template<typename TK>
     class LR_Force

source/source_lcao/module_lr/Grad/force/lr_force_test.cpp

@@ -40,11 +40,14 @@ namespace LR
         {
             const auto& Ds_gs = LR_Util::get_exx_Ds_gs(dm_gs, ucell_, kv, pv_);
             const auto& Ds_gs_2 = LR_Util::get_exx_Ds_gs(dm_gs, ucell_, kv, pv_);
+            // at nspin=1 there is only one channel, and `get_exx_Ds_gs` already returns 0.5*D
+            // (= D_up), so both halves below reuse channel 0.
+            const int is_2nd = (Ds_gs.size() > 1) ? 1 : 0;
             ModuleBase::matrix f_gs_exx(ucell_.nat, 3);
             // test the two function using the two spin channels respectively
             // 0.5 is from dE = 0.5 dTr[D(HD)]. No 0.5 in excited-state calculateion of dTr[(T+Z)(HD)]
             f_gs_exx += cal_force_exx_gs_dm_relaxed_diff(Ds_gs.at(0), Ds_gs_2.at(0), alpha_, std::to_string(0)) * 0.5;    // test passed， = 0.5 groud-state EXX force
-            f_gs_exx += cal_force_exx_dm_trans(Ds_gs.at(1), alpha_, std::to_string(1)) * 0.5;
+            f_gs_exx += cal_force_exx_dm_trans(Ds_gs.at(is_2nd), alpha_, std::to_string(is_2nd)) * 0.5;
             if (PARAM.inp.test_force)
                 ModuleIO::print_force(GlobalV::ofs_running, ucell_, "EXX GS FORCE reproduce (eV/Angstrom)", f_gs_exx, false);
             f_gs_hf_pulay += f_gs_exx;

source/source_lcao/module_lr/operator_casida/operator_lr_exx.cpp

@@ -49,8 +49,9 @@ namespace LR
         }
         case MO_TO_AO_TYPE::CC_oo:
         {
+            // Cv -> CoX^T, i.e. C_o [C_oX^T]
             DMBand<double>(ucell, pmat, this->kv.kvec_c, this->BvK_cells, this->psi_ks_full, this->psi_ks_full)
-                .cal_dm_band(io, io, ik, this->Ds_onebase, 1.0, this->aims_nbasis, this->aims_nbasis);
+                .cal_dm_band(io, iv, ik, this->Ds_onebase, 1.0);
             break;
         }
         case MO_TO_AO_TYPE::CXC_o: