fix segfault and non-hermitian in op_lr_exx

Author: maki49

source/module_lr/dm_trans/dm_trans_parallel.cpp

@@ -109,7 +109,7 @@ std::vector<container::Tensor> cal_dm_trans_pblas(const std::complex<double>* X_
         //     &beta, dm_trans[isk].data<std::complex<double>>(), &i1, &i1, pmat.desc);
 
         // ============== [C_virt * X * C_occ^\dagger]^T=============
-        // ============== = [C_occ^* * X^T * C_virt^T]^T=============
+        // ============== = [C_occ^* * X^T * C_virt^T]=============
         // 1. X*C_occ^\dagger
         char transa = 'N';
         char transb = 'C';

source/module_lr/operator_casida/operator_lr_exx.cpp

@@ -27,6 +27,9 @@ namespace LR
     void OperatorLREXX<double>::cal_DM_onebase(const int io, const int iv, const int ik) const
     {
         ModuleBase::TITLE("OperatorLREXX", "cal_DM_onebase");
+        // NOTICE: DM_onebase will be passed into `cal_energy` interface and conjugated by "zdotc". 
+        // So the formula should be the same as RHS. instead of LHS of the A-matrix, 
+        // i.e. c1v · conj(c2o) ·  e^{-ik(R2-R1)}
         assert(ik == 0);
         for (auto cell : this->BvK_cells)
         {
@@ -42,8 +45,12 @@ namespace LR
                             const int nw2 = aims_nbasis.empty() ? ucell.atoms[it2].nw : aims_nbasis[it2];
                             for (int iw1 = 0;iw1 < nw1;++iw1)
                                 for (int iw2 = 0;iw2 < nw2;++iw2)
-                                    if (this->pmat.in_this_processor(ucell.itiaiw2iwt(it1, ia1, iw1), ucell.itiaiw2iwt(it2, ia2, iw2)))
-                                        D2d(iw1, iw2) = this->psi_ks_full(ik, io, ucell.itiaiw2iwt(it1, ia1, iw1)) * this->psi_ks_full(ik, nocc + iv, ucell.itiaiw2iwt(it2, ia2, iw2));
+                                {
+                                    const int iwt1 = ucell.itiaiw2iwt(it1, ia1, iw1);
+                                    const int iwt2 = ucell.itiaiw2iwt(it2, ia2, iw2);
+                                    if (this->pmat.in_this_processor(iwt1, iwt2))
+                                        D2d(iw1, iw2) = this->psi_ks_full(ik, io, iwt1) * this->psi_ks_full(ik, nocc + iv, iwt2);
+                                }
                         }
         }
     }
@@ -52,6 +59,9 @@ namespace LR
     void OperatorLREXX<std::complex<double>>::cal_DM_onebase(const int io, const int iv, const int ik) const
     {
         ModuleBase::TITLE("OperatorLREXX", "cal_DM_onebase");
+        // NOTICE: DM_onebase will be passed into `cal_energy` interface and conjugated by "zdotc". 
+        // So the formula should be the same as RHS. instead of LHS of the A-matrix, 
+        // i.e. c1v · conj(c2o) ·  e^{-ik(R2-R1)}
         for (auto cell : this->BvK_cells)
         {
             std::complex<double> frac = RI::Global_Func::convert<std::complex<double>>(std::exp(
@@ -68,8 +78,12 @@ namespace LR
                             const int nw2 = aims_nbasis.empty() ? ucell.atoms[it2].nw : aims_nbasis[it2];
                             for (int iw1 = 0;iw1 < nw1;++iw1)
                                 for (int iw2 = 0;iw2 < nw2;++iw2)
-                                    if (this->pmat.in_this_processor(ucell.itiaiw2iwt(it1, ia1, iw1), ucell.itiaiw2iwt(it2, ia2, iw2)))
-                                        D2d(iw1, iw2) = frac * std::conj(this->psi_ks_full(ik, io, ucell.itiaiw2iwt(it1, ia1, iw1))) * this->psi_ks_full(ik, nocc + iv, ucell.itiaiw2iwt(it2, ia2, iw2));
+                                {
+                                    const int iwt1 = ucell.itiaiw2iwt(it1, ia1, iw1);
+                                    const int iwt2 = ucell.itiaiw2iwt(it2, ia2, iw2);
+                                    if (this->pmat.in_this_processor(iwt1, iwt2))
+                                        D2d(iw1, iw2) = frac * std::conj(this->psi_ks_full(ik, io, iwt2)) * this->psi_ks_full(ik, nocc + iv, iwt1);
+                                }
                         }
         }
     }
@@ -78,9 +92,6 @@ namespace LR
     void OperatorLREXX<T>::act(const int nbands, const int nbasis, const int npol, const T* psi_in, T* hpsi, const int ngk_ik, const bool is_first_node)const
     {
         ModuleBase::TITLE("OperatorLREXX", "act");
-
-        const int& nk = this->kv.get_nks() / this->nspin;
-
         // convert parallel info to LibRI interfaces
         std::vector<std::tuple<std::set<TA>, std::set<TA>>> judge = RI_2D_Comm::get_2D_judge(this->pmat);
 
@@ -120,12 +131,12 @@ namespace LR
                 {
                     const int xstart_bk = ik * pX.get_local_size();
                     this->cal_DM_onebase(io, iv, ik);       //set Ds_onebase for all e-h pairs (not only on this processor)
-                    // LR_Util::print_CV(Ds_onebase[is], "Ds_onebase of occ " + std::to_string(io) + ", virtual " + std::to_string(iv) + " in OperatorLREXX", 1e-10);
+                    // LR_Util::print_CV(Ds_onebase, "Ds_onebase of occ " + std::to_string(io) + ", virtual " + std::to_string(iv) + " in OperatorLREXX", 1e-10);
                     const T& ene = 2 * alpha * //minus for exchange(but here plus is right, why?), 2 for Hartree to Ry
                         lri->exx_lri.post_2D.cal_energy(this->Ds_onebase, lri->Hexxs[0]);
                     if (this->pX.in_this_processor(iv, io))
                     {
-                        hpsi[xstart_bk + ik * pX.get_local_size() + this->pX.global2local_col(io) * this->pX.get_row_size() + this->pX.global2local_row(iv)] += ene;
+                        hpsi[xstart_bk + this->pX.global2local_col(io) * this->pX.get_row_size() + this->pX.global2local_row(iv)] += ene;
                     }
                 }
             }

source/module_lr/operator_casida/operator_lr_exx.h

@@ -32,7 +32,7 @@ namespace LR
             const Parallel_Orbitals& pmat_in,
             const double& alpha = 1.0,
             const std::vector<int>& aims_nbasis = {})
-            : nspin(nspin), naos(naos), nocc(nocc), nvirt(nvirt),
+            : nspin(nspin), naos(naos), nocc(nocc), nvirt(nvirt), nk(kv_in.get_nks() / nspin),
             psi_ks(psi_ks_in), DM_trans(DM_trans_in), exx_lri(exx_lri_in), kv(kv_in),
             pX(pX_in), pc(pc_in), pmat(pmat_in), ucell(ucell_in), alpha(alpha),
             aims_nbasis(aims_nbasis)
@@ -42,8 +42,11 @@ namespace LR
             this->is_first_node = false;
 
             // reduce psi_ks for later use
-            this->psi_ks_full.resize(this->kv.get_nks(), nocc + nvirt, this->naos);
-            LR_Util::gather_2d_to_full(this->pc, this->psi_ks.get_pointer(), this->psi_ks_full.get_pointer(), false, this->naos, nocc + nvirt);
+            this->psi_ks_full.resize(this->nk, nocc + nvirt, this->naos);
+            for (int ik = 0;ik < nk;++ik)
+            {
+                LR_Util::gather_2d_to_full(this->pc, &this->psi_ks(ik, 0, 0), &this->psi_ks_full(ik, 0, 0), false, this->naos, nocc + nvirt);
+            }
 
             // get cells in BvK supercell
             const TC period = RI_Util::get_Born_vonKarmen_period(kv_in);
@@ -63,12 +66,13 @@ namespace LR
                          const int ngk_ik = 0,
                          const bool is_first_node = false) const override;
 
-      private:
+    private:
         //global sizes
-        const int& nspin;
-        const int& naos;
-        const int& nocc;
-        const int& nvirt;
+        const int nspin = 1;
+        const int naos = 1;
+        const int nocc = 1;
+        const int nvirt = 1;
+        const int nk = 1;  ///< number of k-points
         const double alpha = 1.0;   //(allow non-ref constant)
         const bool cal_dm_trans = false;
         const bool tdm_sym = false; ///< whether transition density matrix is symmetric

tests/integrate/392_NO_GO_LR_HF/result.ref

@@ -1 +1 @@
-totexcitationenergyref 3.067082
+totexcitationenergyref 3.067058