fix the asymmetry of exx force with 2 different DM

Author: maki49

source/module_base/matrix.cpp

@@ -154,6 +154,16 @@ void matrix::create( const int nrow, const int ncol, const bool flag_zero )
 	}
 }
 
+/* Unary minus*/
+matrix operator-(const matrix& m1)
+{
+    matrix tm(m1);
+    const int size = m1.nr * m1.nc;
+    for (int i = 0; i < size; i++)
+        tm.c[i] = -tm.c[i];
+    return tm;
+}
+
 /* Adding matrices, as a friend */
 matrix operator+(const matrix &m1, const matrix &m2)
 {

source/module_base/matrix.h

@@ -83,7 +83,7 @@ class matrix
 	using type=double;					// Peiae Lin add 2022.08.08 for template
 };
 
-
+matrix operator-(const matrix& m1); // unary minus
 matrix operator+(const matrix &m1, const matrix &m2);
 matrix operator-(const matrix &m1, const matrix &m2);
 matrix operator*(const matrix &m1, const matrix &m2);

source/module_lr/Grad/esolver_lr_grad.cpp

@@ -158,20 +158,23 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
         // 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;
-        // LR_Util::print_value(dm_diff_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
-        for (auto& d : dm_diff_k) { LR_Util::matsym(d.data<T>(), this->nbasis, this->paraMat_); }   // symmetrize
+        std::cout << "dm_diff_k T(k) before symmetrization, istate " + std::to_string(istate) << std::endl;
+        LR_Util::print_value(dm_diff_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
+        // for (auto& d : dm_diff_k) { LR_Util::matsym(d.data<T>(), this->nbasis, this->paraMat_); }   // symmetrize
         // std::cout << "dm_diff_k T(k) after symmetrization, istate " + std::to_string(istate) << std::endl;
         // LR_Util::print_value(dm_diff_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
 
         const std::vector<ct::Tensor>& dm_relaxed_k = cal_dm_trans_pblas(Z.template data<T>() + offset, this->paraX_[ispin], c, this->paraC_, this->nbasis, this->nocc[ispin], this->nvirt[ispin], this->paraMat_);
-        // std::cout << "dm_relaxed_k Z(k) before symmetrization, istate " + std::to_string(istate) << std::endl;
-        // LR_Util::print_value(dm_relaxed_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
+        std::cout << "dm_relaxed_k Z(k) before symmetrization, istate " + std::to_string(istate) << std::endl;
+        LR_Util::print_value(dm_relaxed_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
         for (auto& d : dm_relaxed_k) { LR_Util::matsym(d.data<T>(), this->nbasis, this->paraMat_); }    // symmetrize
-        // std::cout << "dm_relaxed_k Z(k) after symmetrization, istate " + std::to_string(istate) << std::endl;
-        // LR_Util::print_value(dm_relaxed_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
+        std::cout << "dm_relaxed_k Z(k) after symmetrization, istate " + std::to_string(istate) << std::endl;
+        LR_Util::print_value(dm_relaxed_k[0].data<T>(), this->paraMat_.get_col_size(), this->paraMat_.get_row_size());
         // relaxed difference density matrix
         const std::vector<ct::Tensor>& relaxed_diff_dm_k = dm_diff_k + dm_relaxed_k;
+        const elecstate::DensityMatrix<T, T>& diff_dm =
+            LR_Util::build_dm_from_dmk<T, T>(dm_diff_k,
+                this->paraMat_, this->nk, this->kv.kvec_d, this->ucell, this->gd, this->orb_cutoff_);
         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_);
@@ -241,6 +244,20 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
         if (PARAM.inp.test_force)
             ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "OVERLAP-EDM FORCE (eV/Angstrom)", force_overlap_edm, false);
 
+        if (PARAM.inp.test_force)
+        {
+            // test H[T] force (Z=0), non-EXX part
+            elecstate::DensityMatrix<T, double> diff_dm_real(&this->paraMat_, 1, this->kv.kvec_d, this->nk);
+            LR_Util::initialize_DMR(diff_dm_real, this->paraMat_, this->ucell, this->gd, this->orb_cutoff_);
+            LR_Util::get_DMR_real_imag_part(diff_dm, diff_dm_real, 'R');
+
+            GlobalV::ofs_running << "========== [TEST H_GS-(T) force (Z=0), non-EXX part] ===========" << std::endl;
+            ModuleBase::matrix force_hamiltgs_diff = lr_force.cal_force_hamilt_gs_dm_relaxed_diff(diff_dm_real, dm_gs, /*with_ewald=*/false);
+            ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "H_GS-T FORCE (without EXX) (eV/Angstrom)", force_hamiltgs_diff, false);
+            GlobalV::ofs_running << "========== [\\TEST H_GS-(T) force (Z=0), non-EXX part] ===========" << std::endl;
+        }
+
+
 #ifdef __EXX
         const double& alpha = this->exx_info.info_global.hybrid_alpha;
 
@@ -253,15 +270,26 @@ std::vector<ModuleBase::matrix> LR::ESolver_LR<T, TR>::cal_force(const int ispin
             force_hxc_dmtrans += force_exx_dmtrans;
 
         }
+
         if (LR::exx_kernel_list().count(PARAM.inp.dft_functional))
         {
             const auto& Ds_gs = LR_Util::get_exx_Ds_spin1(dm_gs, this->ucell, this->kv, this->paraMat_);    // returns 0.5*D[0]
             const auto& Ds_relaxed_diff = LR_Util::get_exx_Ds_spin1(relaxed_diff_dm, this->ucell, this->kv, this->paraMat_);   // returns 0.5*D[0]
             // LR_Util::print_CV(Ds_relaxed_diff, "Ds_relaxed_diff for EXX force");
-            ModuleBase::matrix force_exx_gs_diff = lr_force.cal_force_exx_gs_dm_relaxed_diff(Ds_gs, Ds_relaxed_diff, alpha * 4.0);  // cancel the two 0.5s in Ds
+            ModuleBase::matrix force_exx_gs_relaxed_diff = lr_force.cal_force_exx_gs_dm_relaxed_diff(Ds_gs, Ds_relaxed_diff, alpha * 4.0);  // cancel the two 0.5s in Ds
+            if (PARAM.inp.test_force)
+                ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "EXX GS-(T+Z) FORCE (eV/Angstrom)", force_exx_gs_relaxed_diff, false);
+            force_hamiltgs_relaxed_diff += force_exx_gs_relaxed_diff;
+
             if (PARAM.inp.test_force)
-                ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "EXX GS-(T+Z) FORCE (eV/Angstrom)", force_exx_gs_diff, false);
-            force_hamiltgs_relaxed_diff += force_exx_gs_diff;
+            {
+                // test H[T] force (Z=0), EXX part
+                const auto& Ds_diff = LR_Util::get_exx_Ds_spin1(diff_dm, this->ucell, this->kv, this->paraMat_);   // returns 0.5*D[0]
+                GlobalV::ofs_running << "========== [TEST H_GS-(T) force (Z=0), EXX part] ===========" << std::endl;
+                ModuleBase::matrix force_exx_gs_diff = lr_force.cal_force_exx_gs_dm_relaxed_diff(Ds_gs, Ds_diff, alpha * 4.0);  // cancel the two 0.5s in Ds
+                ModuleIO::print_force(GlobalV::ofs_running, this->ucell, "H_GS-T EXX FORCE (Z=0) (eV/Angstrom)", force_exx_gs_diff, false);
+                GlobalV::ofs_running << "========== [\\TEST H_GS-(T) force (Z=0), EXX part] ===========" << std::endl;
+            }
         }
 #endif
         forces[istate] = force_hxc_dmtrans + force_hamiltgs_relaxed_diff + force_overlap_edm;
@@ -319,7 +347,7 @@ void LR::ESolver_LR<T, TR>::test_force()
     if (this->nbasis == 2 && ucell.nat == 2)
     {
         // lr_force.cal_H2_sz_center2_deriv(orb_cutoff_, kv);  // for gradient
-        // lr_force.cal_H2_sz_center4(orb_cutoff_, kv, /*is_grad=*/false);  // for Coulomb energy
+        // lr_force.cal_H2_sz_center4(orb_cutoff_, kv, /*is_grad=*/false);  // for 4-center integrals
         // lr_force.cal_H2_sz_center4(orb_cutoff_, kv, /*is_grad=*/true);   // for gradient
         // exit(0);
     }

source/module_lr/Grad/force/lr_force.cpp

@@ -138,8 +138,9 @@ namespace LR
         for (std::size_t idim = 0; idim < 3; ++idim)
             for (const auto& force_item : exx_lri_kernel.force[idim])
                 f_exx_dmtrans(force_item.first, idim) = std::real(force_item.second);
-        const double fac = -2.0 * alpha; //-2 is the same as post_process_Hexx (which didn't act on Hs)
-        return f_exx_dmtrans * fac;
+        const double fac = -2.0 * alpha; //-2 is the same as post_process_Hexx, Hartree to Ry (which didn't act on Hs)
+        return f_exx_dmtrans * fac; // dm_trans (DX) already contain the spin channel (sqrt(2) times of up/down channel DX)
+        // return f_exx_dmtrans * fac * 2; // 2 is the same in post_process_Eexx at nspin=1 ( up->up + down->down, 2 spin-conserving transitions)
     }
 
     template<typename TK>
@@ -151,19 +152,55 @@ namespace LR
     {
         ModuleBase::matrix f_exx_gs_diff(this->ucell_.nat, 3);
         auto& exx_lri_kernel = this->exx_lri_.lock()->get();
-        exx_lri_kernel.set_Ds(dm_gs, this->exx_lri_.lock()->get_info().dm_threshold, spin_suffix);
-        exx_lri_kernel.cal_Hs({ "", "", spin_suffix });  // using dm_gs
-        // // test: print Hs here
-        // LR_Util::print_CV(exx_lri_kernel.Hs, "EXX Hs from gs dm in cal_force_exx_gs_dm_relaxed_diff");
-        // auto* lr_ptr = dynamic_cast<RI::LR<int, int, 3, TK>*>(&exx_lri_kernel);  // wrong: 
-        // assert(lr_ptr != nullptr);
         RI::LR<int, int, 3, TK> lr_exx_kernel(std::move(exx_lri_kernel));
-        lr_exx_kernel.cal_force(relaxed_diff_dm, { "", "", spin_suffix , "", "" }); // using relaxed_diff_dm
-        exx_lri_kernel = std::move(lr_exx_kernel);        // move back 
-        for (std::size_t idim = 0; idim < 3; ++idim)
-            for (const auto& force_item : exx_lri_kernel.force[idim])
-                f_exx_gs_diff(force_item.first, idim) = std::real(force_item.second);
-        const double fac = -2.0 * alpha; //-2 is the same as post_process_Hexx (a.u. to Ry, which didn't act on Hs)
+
+        auto add_force_from_kernel = [&]() {
+            for (std::size_t idim = 0; idim < 3; ++idim)
+                for (const auto& force_item : lr_exx_kernel.force[idim])
+                    f_exx_gs_diff(force_item.first, idim) += std::real(force_item.second);
+            };
+
+        auto transpose_dm = [](const std::map<int, std::map<TAC, RI::Tensor<TK>>>& dm)
+            -> std::map<int, std::map<TAC, RI::Tensor<TK>>>
+            {
+                std::map<int, std::map<TAC, RI::Tensor<TK>>> dm_transpose;
+                for (const auto& pair0 : dm)
+                    for (const auto& pair1 : pair0.second)
+                    {
+                        const int& iat0 = pair0.first;
+                        const int& iat1 = pair1.first.first;
+                        const auto& R = pair1.first.second;
+                        dm_transpose[iat1][{iat0, { -R[0], -R[1], -R[2] }}] = pair1.second.transpose();
+                    }
+                return dm_transpose;
+            };
+
+        // `cal_force` calculates Pulay term. 
+        // If D_IJ = D_KL(H - F = Pulay), it caluclates 0.5 * d(ik | jl).
+        // Multiply spin factor (outside) on it gives the final result.
+        // The spin factor is not hard-coded in this function.
+
+        // 1. Pulay term
+        lr_exx_kernel.set_Ds(dm_gs, this->exx_lri_.lock()->get_info().dm_threshold, spin_suffix);
+        lr_exx_kernel.cal_Hs({ "", "", spin_suffix });  // using dm_gs as D_KL
+        lr_exx_kernel.cal_force(relaxed_diff_dm, { "", "", spin_suffix , "", "" }); // using relaxed_diff_dm as D_IJ
+        add_force_from_kernel();
+
+        // 2. Hellmann-Feynman term
+        const auto& dm_gs_transpose = transpose_dm(dm_gs);
+        const auto& relaxed_diff_dm_transpose = transpose_dm(relaxed_diff_dm);
+        lr_exx_kernel.set_Ds(relaxed_diff_dm_transpose, this->exx_lri_.lock()->get_info().dm_threshold, spin_suffix);
+        lr_exx_kernel.cal_Hs({ "", "", spin_suffix });  // using relaxed_diff_dm as D_KL
+        lr_exx_kernel.cal_force(dm_gs_transpose, { "", "", spin_suffix , "", "" }); // using dm_gs as D_IJ
+        add_force_from_kernel();
+
+        // move back 
+        exx_lri_kernel = std::move(lr_exx_kernel);
+
+        // -2 * 0.5 * alpha
+        // -2 is the same as post_process_Hexx (a.u. to Ry, which didn't act on Hs)
+        // 0.5 is the 2-electron integral prefactor
+        const double fac = -alpha;
         return f_exx_gs_diff * fac;
     }
 #endif