fix xc condition and remove DX symmetrization in cal_force

maki49 · maki49 · commit 7145d5e5984c · 2026-03-04T00:58:16.000+08:00
enable PBE0
diff --git a/source/module_lr/Grad/esolver_lr_grad.cpp b/source/module_lr/Grad/esolver_lr_grad.cpp
@@ -147,14 +147,15 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
         // The imag part will be cancelled in the force calculation, so we use double DM(R) to calculate force. 
         // But complex transition DM(R) is still used in energy density matrix calculation.
         const auto& dm_trans_k = cal_dm_trans_pblas(this->X[ispin].template data<T>() + offset, this->paraX_[ispin], c, this->paraC_, this->nbasis, this->nocc[ispin], this->nvirt[ispin], this->paraMat_);
-        const elecstate::DensityMatrix<T, double>& dm_trans_real =   // D(X), double (FIXME: not enough for periodic system!)
+        auto dm_trans_real =   // D(X), double (FIXME: not enough for periodic system!)
             LR_Util::build_dm_from_dmk<T, double>(dm_trans_k,
                 this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_);
-        const elecstate::DensityMatrix<T, T>& dm_trans =   // D(X) complex
+        LR_Util::transpose_DMR(dm_trans_real, this->ucell.nat); //D(X) is not symmetric, need to transpose for the left side of force calculation
+        auto dm_trans =   // D(X) complex
             LR_Util::build_dm_from_dmk<T, T>(dm_trans_k,
                 this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_);
+        LR_Util::transpose_DMR(dm_trans, this->ucell.nat);
         // LR_Util::print_DMR(dm_trans, "dm_trans of istate " + std::to_string(istate));
-
         // difference density matrix 
         std::vector<ct::Tensor> dm_diff_k = cal_dm_diff_pblas(this->X[ispin].template data<T>() + offset, this->paraX_[ispin], c, this->paraC_, this->nbasis, this->nocc[ispin], this->nvirt[ispin], this->paraMat_);
         // std::cout << "dm_diff_k T(k) before symmetrization, istate " + std::to_string(istate) << std::endl;
@@ -173,7 +174,7 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
         const std::vector<ct::Tensor>& relaxed_diff_dm_k = dm_diff_k + dm_relaxed_k;
         const elecstate::DensityMatrix<T, T>& relaxed_diff_dm =
             LR_Util::build_dm_from_dmk<T, T>(relaxed_diff_dm_k,
-                this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_, /*symmetrize=*/false);
+                this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_);
         // LR_Util::print_DMR(relaxed_diff_dm, "relaxed_diff_dm T+Z (Z symmetrized) of istate " + std::to_string(istate));
 
         // elecstate::DensityMatrix<T, T> relaxed_diff_dm =    // T+D(Z), (R) can be complex
@@ -219,7 +220,7 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
             test_edm_H2<T>(edm_k[0].data<T>(), this->eig_ks.c, c, this->nbasis);
         }
         elecstate::DensityMatrix<T, double> edm_real = LR_Util::build_dm_from_dmk<T, double>(edm_k,
-            this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_);
+            this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_, /*symmetrize=*/true);
         // print edm_real (R)
         if (PARAM.inp.test_force)
         {
diff --git a/source/module_lr/Grad/force/lr_force_test.cpp b/source/module_lr/Grad/force/lr_force_test.cpp
@@ -3,6 +3,10 @@
 #include "pulay_force_hcontainer.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h"   // only for gint terms
 #include "module_lr/utils/lr_util_hcontainer.h"
+#include "module_lr/utils/lr_util_print.h"
+#ifdef __EXX
+#include "module_lr/operator_casida/operator_lr_exx.h"
+#endif
 namespace LR
 {
     template<typename TK>
@@ -34,8 +38,7 @@ namespace LR
         ModuleBase::matrix f_nonortho = cal_force_overlap_edm(edm_gs); // overlap
         this->gint_->reset_DMRGint(1);
 #ifdef __EXX
-        const std::set<std::string> exx_kernel_list = { "hf", "hse" };
-        if (exx_kernel_list.count(PARAM.inp.dft_functional))
+        if (exx_kernel_list().find(PARAM.inp.dft_functional) != exx_kernel_list().end())
         {
             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_);
diff --git a/source/module_lr/Grad/multipliers/cal_multiplier_w_from_z.h b/source/module_lr/Grad/multipliers/cal_multiplier_w_from_z.h
@@ -158,9 +158,12 @@ namespace LR
         cal_dm_diff_relaxed(0, X, Z);  // relaxed difference density matrix T+DZ
         // the 3 terms
         op_ht.act(/*nband=*/1, ld_oo, /*npol=*/1, X, W);    //comment out this line to test H[T+Z]=0
-        if (LR::exx_kernel_list().count(xc_kernel))
+        // std::cout << "W (H[T+Z])) local terms: " << std::endl;
+        // LR_Util::print_value(W, nk, p_occ_occ[0].get_col_size(), p_occ_occ[0].get_row_size());
+        if (LR::exx_kernel_list().count(PARAM.inp.dft_functional))  // H[T+Z] term depends on ground-state kernel (dft_functional)
             op_ht_exx.act(/*nband=*/1, ld_oo, /*npol=*/1, X, W);
-
+        // std::cout << "W (H[T+Z])) local +exx terms: " << std::endl;
+        // LR_Util::print_value(W, nk, p_occ_occ[0].get_col_size(), p_occ_occ[0].get_row_size());
         if (LR_Util::has_local_xc(xc_kernel))
             op_gxc.act(/*nband=*/1, ld_oo, /*npol=*/1, X, W);
 
diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp
@@ -84,7 +84,7 @@ template<typename T, typename TR>
 void LR::ESolver_LR<T, TR>::parameter_check()const
 {
     const std::set<std::string> lr_solvers = { "dav", "lapack" , "spectrum", "dav_subspace", "cg" };
-    const std::set<std::string> xc_kernels = { "rpa", "lda", "pwlda", "pbe", "hf" , "hse" };
+    const std::set<std::string> xc_kernels = { "rpa", "lda", "pwlda", "pbe", "hf" , "hse", "pbe0" };
     if (lr_solvers.find(this->input.lr_solver) == lr_solvers.end())
     {
         throw std::invalid_argument("ESolver_LR: unknown type of lr_solver");
@@ -255,7 +255,7 @@ LR::ESolver_LR<T, TR>::ESolver_LR(ModuleESolver::ESolver_KS_LCAO<T, TR>&& ks_sol
     init_pot(*ks_sol.pelec->charge);
 
 #ifdef __EXX
-    if (xc_kernel == "hf" || xc_kernel == "hse")
+    if (exx_kernel_list().find(xc_kernel) != exx_kernel_list().end())
     {
         // if the same kernel is calculated in the esolver_ks, move it
         std::string dft_functional = LR_Util::tolower(input.dft_functional);
@@ -266,10 +266,10 @@ LR::ESolver_LR<T, TR>::ESolver_LR(ModuleESolver::ESolver_KS_LCAO<T, TR>&& ks_sol
         } else    // construct C, V from scratch
         {
             // 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; }
+            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; }
             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_);
+            this->exx_lri->init(MPI_COMM_WORLD, ucell,this->kv, ks_sol.orb_);   
             this->exx_lri->cal_exx_ions(ucell,input.out_ri_cv);
         }
     }
@@ -442,8 +442,8 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
     // 1. EXX xc_kernel
     // 2. cal_force with ground state with EXX functional
 #ifdef __EXX
-    if (((xc_kernel == "hf" || xc_kernel == "hse") && this->input.lr_solver != "spectrum")
-        || (PARAM.inp.cal_force && (PARAM.inp.dft_functional == "hf" || PARAM.inp.dft_functional == "hse")))
+    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; }
@@ -516,6 +516,7 @@ void LR::ESolver_LR<T, TR>::runner(UnitCell& ucell, const int istep)
             if (input.lr_solver != "lapack") { pre_op.act(1, nloc_per_band, 1, precondition.data(), precondition.data()); }
             // auto spin_types = std::vector<std::string>({ "singlet", "triplet" });
             this->spin_types = { "singlet", "triplet" };
+            // for (int is = 0;is < nspin - 1;++is)
             for (int is = 0;is < nspin;++is)
             {
                 std::cout << "Calculating " << spin_types[is] << " excitations" << std::endl;
@@ -605,6 +606,7 @@ void LR::ESolver_LR<T, TR>::after_all_runners(UnitCell& ucell)
     double lambda_min = std::min(abs_wavelen_range[1], abs_wavelen_range[0]);
     for (int i = 0;i < freq.size();++i) { freq[i] = 91.126664 / (lambda_min + 0.01 * static_cast<double>(i + 1) * lambda_diff); }
     // auto spin_types = (nspin == 2 && !openshell) ? std::vector<std::string>({ "singlet", "triplet" }) : std::vector<std::string>({ "updown" });
+    // for (int is = 0;is < this->X.size() - 1;++is)
     for (int is = 0;is < this->X.size();++is)
     {
         LR_Spectrum<T> spectrum(nspin, this->nbasis, this->nocc, this->nvirt, this->gint_, *this->pw_rho, *this->psi_ks,
diff --git a/source/module_lr/hamilt_casida.h b/source/module_lr/hamilt_casida.h
@@ -99,7 +99,7 @@ namespace LR
                 this->ops->add(lr_hxc);
             }
 #ifdef __EXX
-            if (xc_kernel == "hf" || xc_kernel == "hse")
+            if (exx_kernel_list().find(xc_kernel) != exx_kernel_list().end())
             {   //add Exx operator
                 if (ri_hartree_benchmark != "none" && spin_type == "singlet")
                 {
diff --git a/source/module_lr/hamilt_ulr.hpp b/source/module_lr/hamilt_ulr.hpp
@@ -56,7 +56,7 @@ namespace LR
             this->ops[3]->add(newHxc(1, 1));
 
 #ifdef __EXX
-            if (xc_kernel == "hf" || xc_kernel == "hse")
+            if (exx_kernel_list().find(xc_kernel) != exx_kernel_list().end())
             {
                 std::vector<psi::Psi<T>> psi_ks_spin = { LR_Util::get_psi_spin(psi_ks_in, 0, nk), LR_Util::get_psi_spin(psi_ks_in, 1, nk) };
                 for (int is : {0, 1})
diff --git a/source/module_lr/lr_density.hpp b/source/module_lr/lr_density.hpp
@@ -85,7 +85,7 @@ namespace LR
             // 2. calculate DM(R)
             elecstate::DensityMatrix<T, T> dm_diff=
                 LR_Util::build_dm_from_dmk<T, T>(dm_diff_k,
-                    this->pmat_, this->nk_, this->kv_.kvec_d, this->ucell_, this->gd_, this->orb_cutoff_, /*symmetrize=*/false);
+                    this->pmat_, this->nk_, this->kv_.kvec_d, this->ucell_, this->gd_, this->orb_cutoff_);
             // 3. calculate electron density from DM(R)
             this->dm_to_density(dm_diff, density);
         }
diff --git a/source/module_lr/lr_spectrum.cpp b/source/module_lr/lr_spectrum.cpp
@@ -30,6 +30,7 @@ elecstate::DensityMatrix<T, T> LR::LR_Spectrum<T>::cal_transition_density_matrix
     {
         LR_Util::initialize_DMR(DM_trans, this->pmat, this->ucell, this->gd_, this->orb_cutoff_);
         DM_trans.cal_DMR();
+        LR_Util::swap_atompair_in_DMR(DM_trans, ucell.nat);   // make D(R) consistent with the defination: D(R)[iat1][iat2] = \sum_k c1(k)c2^*(k)exp(-ik(R2-R1))
     }
     return DM_trans;
 }
diff --git a/source/module_lr/operator_casida/operator_lr_exx.h b/source/module_lr/operator_casida/operator_lr_exx.h
@@ -8,7 +8,7 @@
 #include "module_lr/Grad/dm_diff/dm_diff.h"
 namespace LR
 {
-    inline std::set<std::string> exx_kernel_list() { return { "hf", "hse" }; };
+    inline std::set<std::string> exx_kernel_list() { return { "hf", "hse", "pbe0" }; };
     template<typename T = double>
     class OperatorLREXX : public hamilt::Operator<T, base_device::DEVICE_CPU>
     {
diff --git a/source/module_lr/operator_casida/operator_lr_hxc.cpp b/source/module_lr/operator_casida/operator_lr_hxc.cpp
@@ -26,6 +26,7 @@ namespace LR
         const int& lgd = gint->gridt->lgd;
 
         this->DM_trans.cal_DMR();  //DM_trans.get_DMR_vector() is 2d-block parallized
+        LR_Util::swap_atompair_in_DMR(this->DM_trans, ucell.nat);   // make D(R) consistent with the defination: D(R)[iat1][iat2] = \sum_k c1(k)c2^*(k)exp(-ik(R2-R1))
         // LR_Util::print_DMR(DM_trans, ucell.nat, "DMR");
 
         // ========================= begin grid calculation=========================
diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp
@@ -20,7 +20,7 @@ namespace LR
         xc_kernel_components_(rho_basis, ucell, chg_gs, pgrid, nspin_, xc_kernel, lr_init_xc_kernel, (st == SpinType::S2_updown)), //call XC_Functional::set_func_type and libxc
         xc_type_(XCType(XC_Functional::get_func_type()))
     {
-        if (std::set<std::string>({ "lda", "pwlda", "pbe", "hse" }).count(xc_kernel)) { this->set_integral_func(this->spin_type_, this->xc_type_); }
+        if (LR_Util::has_local_xc(xc_kernel)) { this->set_integral_func(this->spin_type_, this->xc_type_); }
     }
 
     void PotHxcLR::cal_v_eff(double** rho, const UnitCell& ucell, ModuleBase::matrix& v_eff, const std::vector<int>& ispin_op) const
diff --git a/source/module_lr/potentials/xc_kernel.cpp b/source/module_lr/potentials/xc_kernel.cpp
@@ -24,7 +24,7 @@ LR::KernelXC::KernelXC(const ModulePW::PW_Basis& rho_basis,
     const std::vector<std::string>& lr_init_xc_kernel,
     const bool openshell) :rho_basis_(rho_basis), openshell_(openshell)
 {
-    if (!std::set<std::string>({ "lda", "pwlda", "pbe", "hse" }).count(kernel_name)) { return; }
+    if (!LR_Util::has_local_xc(kernel_name)) { return; }
     XC_Functional::set_xc_type(kernel_name);    // for hse, (1-alpha) and omega are set here
 
     const int& nrxx = rho_basis.nrxx;
diff --git a/source/module_lr/utils/lr_util.cpp b/source/module_lr/utils/lr_util.cpp
@@ -95,6 +95,42 @@ namespace LR_Util
         pztranc_(&n, &n, &alpha, tmp.data(), &i1, &i1, pmat.desc, &beta, inout, &i1, &i1, pmat.desc);
     }
 
+    template<>
+    void mattrans<double>(const double* in, const int n, const Parallel_2D& pmat, double* out)
+    {
+        std::copy(in, in + pmat.get_local_size(), out);
+        const double alpha = 1.0, beta = 0.0;
+        const int i1 = 1;
+        pdtran_(&n, &n, &alpha, in, &i1, &i1, pmat.desc, &beta, out, &i1, &i1, pmat.desc);
+    }
+    template<>
+    void mattrans<double>(double* inout, const int n, const Parallel_2D& pmat)
+    {
+        std::vector<double> tmp(pmat.get_local_size());
+        std::copy(inout, inout + pmat.get_local_size(), tmp.begin());
+        const double alpha = 1.0, beta = 0.0;
+        const int i1 = 1;
+        pdtran_(&n, &n, &alpha, tmp.data(), &i1, &i1, pmat.desc, &beta, inout, &i1, &i1, pmat.desc);
+    }
+    template<>
+    void mattrans<std::complex<double>>(const std::complex<double>* in, const int n, const Parallel_2D& pmat, std::complex<double>* out)
+    {
+        std::copy(in, in + pmat.get_local_size(), out);
+        const std::complex<double> alpha(1.0, 0.0), beta(0.0, 0.0);
+        const int i1 = 1;
+        pztranc_(&n, &n, &alpha, in, &i1, &i1, pmat.desc, &beta, out, &i1, &i1, pmat.desc);
+    }
+    template<>
+    void mattrans<std::complex<double>>(std::complex<double>* inout, const int n, const Parallel_2D& pmat)
+    {
+        std::vector<std::complex<double>> tmp(pmat.get_local_size());
+        std::copy(inout, inout + pmat.get_local_size(), tmp.begin());
+        const std::complex<double> alpha(1.0, 0.0), beta(0.0, 0.0);
+        const int i1 = 1;
+        pztranc_(&n, &n, &alpha, tmp.data(), &i1, &i1, pmat.desc, &beta, inout, &i1, &i1, pmat.desc);
+    }
+
+
     template<>
     void matantisym<double>(double* inout, const int n, const Parallel_2D& pmat)
     {
diff --git a/source/module_lr/utils/lr_util.h b/source/module_lr/utils/lr_util.h
@@ -29,7 +29,7 @@ namespace LR_Util
     /// @brief check if the xc functional has local xc kernel
     inline bool has_local_xc(const std::string& name)
     {
-        return std::set<std::string>({ "lda", "pwlda", "pbe", "hse" }).count(name);
+        return std::set<std::string>({ "lda", "pwlda", "pbe", "hse", "pbe0" }).count(name);
     }
 
     /// @brief calculate the number of electrons
@@ -80,6 +80,10 @@ namespace LR_Util
     void matsym(const T* in, const int n, const Parallel_2D& pmat, T* out);
     template<typename T>
     void matsym(T* inout, const int n, const Parallel_2D& pmat);
+    template<typename T>
+    void mattrans(const T* in, const int n, const Parallel_2D& pmat, T* out);
+    template<typename T>
+    void mattrans(T* inout, const int n, const Parallel_2D& pmat);
 
     // calculate (A-A^T)/2 (in-place version)
     template<typename T>
diff --git a/source/module_lr/utils/lr_util_hcontainer.h b/source/module_lr/utils/lr_util_hcontainer.h
@@ -176,6 +176,48 @@ namespace LR_Util
     }
 
 
+    template<typename TK, typename TR>
+    void swap_atompair_in_DMR(const elecstate::DensityMatrix<TK, TR>& dm, const int nat)
+    {
+        for (int iat1 = 0; iat1 < nat; ++iat1)
+            for (int iat2 = iat1 + 1; iat2 < nat; ++iat2)
+                for (auto& dr : dm.get_DMR_vector())
+                {
+                    auto ap1 = dr->find_pair(iat1, iat2);
+                    auto ap2 = dr->find_pair(iat2, iat1);
+                    if (ap1 && ap2)
+                        std::swap(ap1, ap2);
+                }
+    }
+
+    template<typename TK>
+    void transpose_DMR(elecstate::DensityMatrix<TK, double>& dm, const int nat)
+    {
+        auto pv = dm.get_paraV_pointer();
+        // 1. transpose dm(k)
+        for (auto& dk : dm.get_DMK_vector())
+            LR_Util::mattrans(dk.data(), pv->get_global_row_size(), *pv);
+
+        // 2. FT
+        dm.cal_DMR();
+        // 3. swap atom pair (iat1, iat2) to (iat2, iat1)
+        swap_atompair_in_DMR(dm, nat);
+    }
+    template<typename TK>
+    void transpose_DMR(elecstate::DensityMatrix<TK, std::complex<double>>& dm, const int nat)
+    {
+        throw std::runtime_error("transpose_DMR is not implemented for complex DMR, due to the lack of minus-sign FT.");
+        auto pv = dm.get_paraV_pointer();
+        // 1. dm(k) dagger
+        for (auto& dk : dm.get_DMK_vector())
+            LR_Util::mattrans(dk.data(), pv->get_global_row_size(), *pv);
+
+        // 2. FT with the minus sign in the exponent (TO DO)
+        dm.cal_DMR();
+        // 3. swap atom pair (iat1, iat2) to (iat2, iat1)
+        swap_atompair_in_DMR(dm, nat);
+    }
+
     template<typename TK, typename TR>
     elecstate::DensityMatrix<TK, TR> build_dm_from_dmk(const std::vector<ct::Tensor>& dmk,
         const Parallel_Orbitals& pmat,
@@ -184,8 +226,9 @@ namespace LR_Util
         const UnitCell& ucell,
         const Grid_Driver& gd,
         const std::vector<double>& orb_cutoff,
-        const bool symmetrize = true,
-        const bool cal_dmr = true)
+        const bool symmetrize = false,
+        const bool cal_dmr = true,
+        const bool transpose = false)
     {
         elecstate::DensityMatrix<TK, TR> dm(&pmat, 1, kvec_d, nk);
         initialize_DMR(dm, pmat, ucell, gd, orb_cutoff);
@@ -197,7 +240,11 @@ namespace LR_Util
         for (int ik = 0; ik < nk; ++ik)
             dm.set_DMK_pointer(ik, dmk[ik].data<TK>());
 
-        if (cal_dmr) { dm.cal_DMR(); }
+        if (cal_dmr)
+        {
+            dm.cal_DMR();
+            LR_Util::swap_atompair_in_DMR(dm, ucell.nat);   // make D(R) consistent with the defination: D(R)[iat1][iat2] = \sum_k c1(k)c2^*(k)exp(-ik(R2-R1))
+        }
         return dm;
     }
 

Original file line number	Diff line number	Diff line change
`@@ -99,7 +99,7 @@ namespace LR`
`99`	`99`	`this->ops->add(lr_hxc);`
`100`	`100`	`}`
`101`	`101`	`#ifdef __EXX`
`102`		`- if (xc_kernel == "hf" \|\| xc_kernel == "hse")`
	`102`	`+ if (exx_kernel_list().find(xc_kernel) != exx_kernel_list().end())`
`103`	`103`	`{ //add Exx operator`
`104`	`104`	`if (ri_hartree_benchmark != "none" && spin_type == "singlet")`
`105`	`105`	`{`
Original file line number	Diff line number	Diff line change
`@@ -56,7 +56,7 @@ namespace LR`
`56`	`56`	`this->ops[3]->add(newHxc(1, 1));`
`57`	`57`
`58`	`58`	`#ifdef __EXX`
`59`		`- if (xc_kernel == "hf" \|\| xc_kernel == "hse")`
	`59`	`+ if (exx_kernel_list().find(xc_kernel) != exx_kernel_list().end())`
`60`	`60`	`{`
`61`	`61`	`std::vector<psi::Psi<T>> psi_ks_spin = { LR_Util::get_psi_spin(psi_ks_in, 0, nk), LR_Util::get_psi_spin(psi_ks_in, 1, nk) };`
`62`	`62`	`for (int is : {0, 1})`
Original file line number	Diff line number	Diff line change
`@@ -85,7 +85,7 @@ namespace LR`
`85`	`85`	`// 2. calculate DM(R)`
`86`	`86`	`elecstate::DensityMatrix<T, T> dm_diff=`
`87`	`87`	`LR_Util::build_dm_from_dmk<T, T>(dm_diff_k,`
`88`		`- this->pmat_, this->nk_, this->kv_.kvec_d, this->ucell_, this->gd_, this->orb_cutoff_, /symmetrize=/false);`
	`88`	`+ this->pmat_, this->nk_, this->kv_.kvec_d, this->ucell_, this->gd_, this->orb_cutoff_);`
`89`	`89`	`// 3. calculate electron density from DM(R)`
`90`	`90`	`this->dm_to_density(dm_diff, density);`
`91`	`91`	`}`
Original file line number	Diff line number	Diff line change
`@@ -30,6 +30,7 @@ elecstate::DensityMatrix<T, T> LR::LR_Spectrum<T>::cal_transition_density_matrix`
`30`	`30`	`{`
`31`	`31`	`LR_Util::initialize_DMR(DM_trans, this->pmat, this->ucell, this->gd_, this->orb_cutoff_);`
`32`	`32`	`DM_trans.cal_DMR();`
	`33`	`+ LR_Util::swap_atompair_in_DMR(DM_trans, ucell.nat); // make D(R) consistent with the defination: D(R)[iat1][iat2] = \sum_k c1(k)c2^*(k)exp(-ik(R2-R1))`
`33`	`34`	`}`
`34`	`35`	`return DM_trans;`
`35`	`36`	`}`
Original file line number	Diff line number	Diff line change
`@@ -8,7 +8,7 @@`
`8`	`8`	`#include "module_lr/Grad/dm_diff/dm_diff.h"`
`9`	`9`	`namespace LR`
`10`	`10`	`{`
`11`		`- inline std::set<std::string> exx_kernel_list() { return { "hf", "hse" }; };`
	`11`	`+ inline std::set<std::string> exx_kernel_list() { return { "hf", "hse", "pbe0" }; };`
`12`	`12`	`template<typename T = double>`
`13`	`13`	`class OperatorLREXX : public hamilt::Operator<T, base_device::DEVICE_CPU>`
`14`	`14`	`{`