fix RHS of Z-eq

Author: maki49

source/module_hamilt_general/operator.h

@@ -27,6 +27,8 @@ enum class calculation_type
     lcao_dftu,
     lcao_sc_lambda,
     lcao_tddft_velocity,
+    lr_dmtrans_vo,
+    lr_dmdiff_vo
 };
 
 // Basic class for operator module,

source/module_lr/Grad/force/lr_force.cpp

@@ -32,7 +32,9 @@ namespace LR
         const Charge chr_diff_relaxed = dm_to_charge(relax_diff_dm);
 
         // 1. local pp (Hellmann-Feynman)(fvl_dvl) + ewald + core correction (+ self-consistent charge)
-        ModuleBase::matrix f_pw = ForcePWTerms<double>()(this->ucell_, chr_diff_relaxed, this->rhopw_, this->locpp_, this->sf_, with_ewald);
+        ModuleBase::matrix f_pw = PARAM.inp.vl_in_h ?
+            ForcePWTerms<double>()(this->ucell_, chr_diff_relaxed, this->rhopw_, this->locpp_, this->sf_, with_ewald) :
+            ModuleBase::matrix(this->ucell_.nat, 3);
 
         // 2. nonlocal pp (Hellmann-Feynman + Pulay)
         ModuleBase::matrix fvnl = cal_force_nonlocal(this->ucell_, this->kvec_d_, this->gd_, this->two_center_bundle_, relax_diff_dm);

source/module_lr/Grad/multipliers/hamilt_zeq_left.h

@@ -46,7 +46,7 @@ namespace LR
             // 1. $2\sum_bX_{ib}K_{ab}[D^X]-2\sum_jX_{ja}K_{ij}[D^X]$
             this->ops = new OperatorLRDiag<T>(eig_ks.c, pX[0], kv.get_nks() / nspin, nocc[0], nvirt[0]);
             // 2. $H_{ia}[D^Z]$, equals to $2K_{ab}[D^Z]$ when $D^Z$ is symmetrized
-            OperatorLRHxc<T>* op_hz = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
+            hamilt::Operator<T>* op_hz = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
                 *this->DM_trans, gint, pot_hxc_gs, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
                 { 0 }, 2.0, ATYPE::CC_vo);
             this->ops->add(op_hz);

source/module_lr/Grad/multipliers/hamilt_zeq_right.h

@@ -51,15 +51,18 @@ namespace LR
                 *this->DM_trans, gint, pot, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
                 { 0 }, -2.0, ATYPE::CXC);
             // 2. $H_{ia}[T]$, equals to $2K_{ab}[T]$ when $T$ is symmetrized
-            OperatorLRHxc<T>* op_ht = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
+            hamilt::Operator<T>* op_ht = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
                 *this->DM_diff, gint, pot_hxc_gs, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
-                { 0 }, T(-2.0), ATYPE::CC_vo);
+                { 0 }, T(-2.0), ATYPE::CC_vo, hamilt::calculation_type::lr_dmdiff_vo);
             this->ops->add(op_ht);
             // 3. $2\sum_{jb,kc} g^{xc}_{ia, jb, kc}X_{jb}X_{kc}$
             this->pot_grad = std::make_shared<PotGradXCLR>(pot.lock()->xc_kernel_components, pot.lock()->get_rho_basis(), ucell, pot.lock()->nrxx);
-            OperatorLRHxc<T>* op_gxc = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
+            // !!op_gxc has some bug
+            hamilt::Operator<T>* op_gxc = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
                 *this->DM_trans, gint, this->pot_grad, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
-                { 0 }, T(-2.0), ATYPE::CC_vo);
+                { 0 }, T(-2.0), ATYPE::CC_vo, hamilt::calculation_type::lr_dmtrans_vo);
+            assert(op_gxc != nullptr);
+            std::cout << "op_gxc=" << op_ht << std::endl;
             this->ops->add(op_gxc);
             // // test: op_ht only 
             // delete this->ops;
@@ -87,8 +90,10 @@ namespace LR
                     const auto psi_ks_is = LR_Util::get_psi_spin(psi_ks, is, this->nk);
 #ifdef __MPI
                     std::vector<ct::Tensor> dm_diff_2d = cal_dm_diff_pblas(X, this->pX[is], psi_ks_is, pc, naos, nocc[is], nvirt[is], pmat);
+                    for (auto& t : dm_diff_2d) LR_Util::matsym(t.data<T>(), naos, pmat);
 #else
                     std::vector<ct::Tensor> dm_diff_2d = cal_dm_diff_blas(X, psi_ks_is, naos, nocc[is], nvirt[is]);
+                    for (auto& t : dm_diff_2d) LR_Util::matsym(t.data<T>(), naos);
 #endif
                     for (int ik = 0;ik < this->nk;++ik) { this->DM_diff->set_DMK_pointer(ik, dm_diff_2d[ik].data<T>()); }
                     // std::cout << "difference density matrix" << std::endl;
@@ -103,7 +108,7 @@ namespace LR
             for (int ib = 0;ib < nband;++ib)
             {
                 const int offset = ib * ld_psi;
-                // this->cal_dm_trans(0, psi + offset);  // transition density matrix, only for test
+                this->cal_dm_trans(0, psi + offset);  // transition density matrix, only for test
                 this->cal_dm_diff(0, psi + offset);  // difference density matrix
                 hamilt::Operator<T>* node(this->ops);
                 while (node != nullptr)

source/module_lr/Grad/multipliers/zeq_solver.hpp

@@ -100,7 +100,7 @@ namespace LR
 #endif 
             gint, pot, pot_hxc_gs, kv, px, pc, pmat, spin_type);
         ModuleBase::timer::tick("Z_vector", "Z_vector_R");
-        ops_R.hPsi(X, R.data<T>(), nk * px[0].get_local_size(), nstates);  // act each operators on X
+        ops_R.hPsi(X, R.data<T>(), nloc_per_band, nstates);  // act each operators on X
         ModuleBase::timer::tick("Z_vector", "Z_vector_R");
 
         // 2. the left-hand side of Z-vector equation

source/module_lr/hamilt_casida.h

@@ -94,7 +94,7 @@ namespace LR
             else
 #endif
             {
-                OperatorLRHxc<T>* lr_hxc = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks_in,
+                hamilt::Operator<T>* lr_hxc = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks_in,
                     *this->DM_trans, gint_in, pot_in, ucell_in, orb_cutoff, gd_in, kv_in, pX_in, pc_in, pmat_in);
                 this->ops->add(lr_hxc);
             }

source/module_lr/operator_casida/operator_lr_hxc.h

@@ -39,14 +39,15 @@ namespace LR
             const Parallel_Orbitals& pmat_in,
             const std::vector<int>& ispin_ks = { 0 },
             const T factor_in = (T)1.0,
-            const MO_TO_AO_TYPE dm_pq_in = MO_TO_AO_TYPE::CC_vo)
-          : nspin(nspin), naos(naos), nocc(nocc), nvirt(nvirt), nk(kv_in.get_nks() / nspin), psi_ks(psi_ks_in),
+            const MO_TO_AO_TYPE dm_pq_in = MO_TO_AO_TYPE::CC_vo,
+            const hamilt::calculation_type cal_type_in = hamilt::calculation_type::lcao_gint)
+            : nspin(nspin), naos(naos), nocc(nocc), nvirt(nvirt), nk(kv_in.get_nks() / nspin), psi_ks(psi_ks_in),
             DM_trans(DM_trans_in), gint(gint_in), pot(pot_in), ucell(ucell_in), orb_cutoff_(orb_cutoff), gd(gd_in),
             kv(kv_in), pX(pX_in), pc(pc_in), pmat(pmat_in), ispin_ks(ispin_ks),
             factor_(factor_in), dm_pq_(dm_pq_in)
         {
           ModuleBase::TITLE("OperatorLRHxc", "OperatorLRHxc");
-          this->cal_type = hamilt::calculation_type::lcao_gint;
+          this->cal_type = cal_type_in;
           this->is_first_node = true;
           this->hR = std::unique_ptr<hamilt::HContainer<T>>(new hamilt::HContainer<T>(&pmat_in));
           LR_Util::initialize_HR<T, T>(*this->hR, ucell_in, gd_in, orb_cutoff);

source/module_lr/utils/lr_util_hcontainer.h

@@ -6,6 +6,7 @@
 #include <ATen/core/tensor.h>
 #include "module_parameter/parameter.h"
 #include "module_io/single_R_io.h"
+#include "module_lr/utils/lr_util.h"
 namespace LR_Util
 {
     template <typename T>
@@ -186,6 +187,8 @@ namespace LR_Util
         initialize_DMR(dm, pmat, ucell, gd, orb_cutoff);
         for (int ik = 0; ik < nk; ++ik)
         {
+            // symmetrize
+            LR_Util::matsym(dmk[ik].data<TK>(), pmat.get_global_row_size(), pmat);
             dm.set_DMK_pointer(ik, dmk[ik].data<TK>());
         }
         if (cal_dmr) { dm.cal_DMR(); }