restrict DM size within orb_rcut

Author: maki49

source/module_lr/esolver_lrtd_lcao.cpp

@@ -219,6 +219,7 @@ LR::ESolver_LR<T, TR>::ESolver_LR(ModuleESolver::ESolver_KS_LCAO<T, TR>&& ks_sol
     if (std::is_same<T, double>::value) { this->gint_g_ = std::move(ks_sol.GG); }
     else { this->gint_k_ = std::move(ks_sol.GK); }
     this->set_gint();
+    this->gint_->reset_DMRGint(1);
 
     // move pw basis
     delete this->pw_rho;    // newed in ESolver_FP::ESolver_FP
@@ -349,7 +350,6 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
         PARAM.inp.test_atom_input);
     this->set_gint();
     this->gint_->gridt = &this->gt_;
-    this->gint_->reset_DMRGint(1);
 
     // (3) Periodic condition search for each grid.
     double dr_uniform = 0.001;
@@ -505,7 +505,7 @@ void LR::ESolver_LR<T, TR>::after_all_runners()
     //cal spectrum
     std::vector<double> freq(100);
     std::vector<double> abs_wavelen_range({ 20, 200 });//default range
-    if (input.abs_wavelen_range.size() == 2 && std::abs(input.abs_wavelen_range[1] - input.abs_wavelen_range[0]) > 0.02)
+    if (input.abs_wavelen_range.size() >= 2 && std::abs(input.abs_wavelen_range[1] - input.abs_wavelen_range[0]) > 0.02)
     {
         abs_wavelen_range = input.abs_wavelen_range;
     }
@@ -516,7 +516,8 @@ void LR::ESolver_LR<T, TR>::after_all_runners()
     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,
-            this->ucell, this->kv, this->paraX_, this->paraC_, this->paraMat_,
+            this->ucell, this->kv, GlobalC::GridD, this->orb_cutoff_,
+            this->paraX_, this->paraC_, this->paraMat_,
             &this->pelec->ekb.c[is * nstates], this->X[is].template data<T>(), nstates, openshell);
         spectrum.transition_analysis(spin_types[is]);
         spectrum.optical_absorption(freq, input.abs_broadening, spin_types[is]);

source/module_lr/hamilt_casida.h

@@ -46,7 +46,8 @@ namespace LR
             if (ri_hartree_benchmark != "aims") { assert(aims_nbasis.empty()); }
             // always use nspin=1 for transition density matrix
             this->DM_trans = LR_Util::make_unique<elecstate::DensityMatrix<T, T>>(&pmat_in, 1, kv_in.kvec_d, nk);
-            this->DM_trans->init_DMR(&gd_in, &ucell_in);
+            LR_Util::initialize_DMR(*this->DM_trans, pmat_in, ucell_in, gd_in, orb_cutoff);
+            // this->DM_trans->init_DMR(&gd_in, &ucell_in); // too large due to not restricted by orb_cutoff
 
             // add the diag operator  (the first one)
             this->ops = new OperatorLRDiag<T>(eig_ks.c, pX[0], nk, nocc[0], nvirt[0]);

source/module_lr/hamilt_ulr.hpp

@@ -40,11 +40,8 @@ namespace LR
         {
             ModuleBase::TITLE("HamiltULR", "HamiltULR");
             this->DM_trans = LR_Util::make_unique<elecstate::DensityMatrix<T, T>>(&pmat_in, 1, kv_in.kvec_d, nk);
-            this->DM_trans->init_DMR(&gd_in, &ucell_in);
-
-            // how to change the index of eig_ks and psi_ks_In?
-            // modify the interface of opetators to support different left- and right- spin-pairs
-
+            LR_Util::initialize_DMR(*this->DM_trans, pmat_in, ucell_in, gd_in, orb_cutoff);
+            // this->DM_trans->init_DMR(&gd_in, &ucell_in); // too large due to not restricted by orb_cutoff
             this->ops.resize(4);
 
             this->ops[0] = new OperatorLRDiag<T>(eig_ks.c, pX_in[0], nk, nocc[0], nvirt[0]);

source/module_lr/lr_spectrum.cpp

@@ -24,14 +24,14 @@ inline void check_sum_rule(const double& osc_tot)
 }
 
 template<>
-void LR::LR_Spectrum<double>::oscillator_strength()
+void LR::LR_Spectrum<double>::oscillator_strength(Grid_Driver& gd, const std::vector<double>& orb_cutoff)
 {
     ModuleBase::TITLE("LR::LR_Spectrum", "oscillator_strength");
     std::vector<double>& osc = this->oscillator_strength_;  // unit: Ry
     osc.resize(nstate, 0.0);
     double osc_tot = 0.0;
     elecstate::DensityMatrix<double, double> DM_trans(&this->pmat, 1, this->kv.kvec_d, this->nk);
-    DM_trans.init_DMR(&GlobalC::GridD, &this->ucell);
+    LR_Util::initialize_DMR(DM_trans, this->pmat, this->ucell, gd, orb_cutoff);
     this->transition_dipole_.resize(nstate, ModuleBase::Vector3<double>(0.0, 0.0, 0.0));
     for (int istate = 0;istate < nstate;++istate)
     {
@@ -80,16 +80,16 @@ void LR::LR_Spectrum<double>::oscillator_strength()
 }
 
 template<>
-void LR::LR_Spectrum<std::complex<double>>::oscillator_strength()
+void LR::LR_Spectrum<std::complex<double>>::oscillator_strength(Grid_Driver& gd, const std::vector<double>& orb_cutoff)
 {
     ModuleBase::TITLE("LR::LR_Spectrum", "oscillator_strength");
     std::vector<double>& osc = this->oscillator_strength_;  // unit: Ry
     osc.resize(nstate, 0.0);
     double osc_tot = 0.0;
     elecstate::DensityMatrix<std::complex<double>, std::complex<double>> DM_trans(&this->pmat, 1, this->kv.kvec_d, this->nk);
-    DM_trans.init_DMR(&GlobalC::GridD, &this->ucell);
+    LR_Util::initialize_DMR(DM_trans, this->pmat, this->ucell, gd, orb_cutoff);
     elecstate::DensityMatrix<std::complex<double>, double> DM_trans_real_imag(&this->pmat, 1, this->kv.kvec_d, this->nk);
-    DM_trans_real_imag.init_DMR(&GlobalC::GridD, &this->ucell);
+    LR_Util::initialize_DMR(DM_trans_real_imag, this->pmat, this->ucell, gd, orb_cutoff);
 
     this->transition_dipole_.resize(nstate, ModuleBase::Vector3<std::complex<double>>(0.0, 0.0, 0.0));
     for (int istate = 0;istate < nstate;++istate)

source/module_lr/lr_spectrum.h

@@ -12,7 +12,8 @@ namespace LR
     public:
         LR_Spectrum(const int& nspin_global, const int& naos, const std::vector<int>& nocc, const std::vector<int>& nvirt,
             typename TGint<T>::type* gint, const ModulePW::PW_Basis& rho_basis, psi::Psi<T>& psi_ks_in,
-            const UnitCell& ucell, const K_Vectors& kv_in, const std::vector<Parallel_2D>& pX_in, const Parallel_2D& pc_in, const Parallel_Orbitals& pmat_in,
+            const UnitCell& ucell, const K_Vectors& kv_in, Grid_Driver& gd, const std::vector<double>& orb_cutoff,
+            const std::vector<Parallel_2D>& pX_in, const Parallel_2D& pc_in, const Parallel_Orbitals& pmat_in,
             const double* eig, const T* X, const int& nstate, const bool& openshell) :
             nspin_x(openshell ? 2 : 1), naos(naos), nocc(nocc), nvirt(nvirt), nk(kv_in.get_nks() / nspin_global),
             gint(gint), rho_basis(rho_basis), ucell(ucell), kv(kv_in),
@@ -22,15 +23,15 @@ namespace LR
             gdim(nk* std::inner_product(nocc.begin(), nocc.end(), nvirt.begin(), 0))
         {
             for (int is = 0;is < nspin_global;++is) { psi_ks.emplace_back(LR_Util::get_psi_spin(psi_ks_in, is, nk)); }
-            this->oscillator_strength();
+            this->oscillator_strength(gd, orb_cutoff);
         };
         /// @brief calculate the optical absorption spectrum
         void optical_absorption(const std::vector<double>& freq, const double eta, const std::string& spintype);
         /// @brief print out the transition dipole moment and the main contributions to the transition amplitude
         void transition_analysis(const std::string& spintype);
     private:
         /// $$2/3\Omega\sum_{ia\sigma} |\braket{\psi_{i}|\mathbf{r}|\psi_{a}} |^2\int \rho_{\alpha\beta}(\mathbf{r}) \mathbf{r} d\mathbf{r}$$
-        void oscillator_strength();
+        void oscillator_strength(Grid_Driver& gd, const std::vector<double>& orb_cutoff);
         const int nspin_x = 1;   ///< 1 for singlet/triplet, 2 for updown(openshell)
         const int naos = 1;
         const std::vector<int>& nocc;

source/module_lr/operator_casida/operator_lr_hxc.cpp

@@ -94,7 +94,7 @@ namespace LR
         elecstate::DensityMatrix<std::complex<double>, double> DM_trans_real_imag(&pmat, 1, kv.kvec_d, kv.get_nks() / nspin);
         DM_trans_real_imag.init_DMR(*this->hR);
         hamilt::HContainer<double> HR_real_imag(GlobalC::ucell, &this->pmat);
-        this->initialize_HR(HR_real_imag, ucell, gd);
+        LR_Util::initialize_HR<std::complex<double>, double>(HR_real_imag, ucell, gd, orb_cutoff_);
 
         auto dmR_to_hR = [&, this](const char& type) -> void
             {

source/module_lr/operator_casida/operator_lr_hxc.h

@@ -40,7 +40,8 @@ namespace LR
             this->cal_type = hamilt::calculation_type::lcao_gint;
             this->is_first_node = true;
             this->hR = std::unique_ptr<hamilt::HContainer<T>>(new hamilt::HContainer<T>(&pmat_in));
-            this->initialize_HR(*this->hR, ucell_in, gd_in);
+            LR_Util::initialize_HR<T, T>(*this->hR, ucell_in, gd_in, orb_cutoff);
+            assert(&pmat_in == this->hR->get_paraV());
         };
         ~OperatorLRHxc() { };
 
@@ -49,33 +50,6 @@ namespace LR
         virtual void act(const int nbands, const int nbasis, const int npol, const T* psi_in, T* hpsi, const int ngk_ik = 0)const override;
 
     private:
-        template<typename TR>   //T=double, TR=double; T=std::complex<double>, TR=std::complex<double>/double
-        void initialize_HR(hamilt::HContainer<TR>& hR, const UnitCell& ucell, Grid_Driver& gd) const
-        {
-            for (int iat1 = 0; iat1 < ucell.nat; iat1++)
-            {
-                auto tau1 = ucell.get_tau(iat1);
-                int T1, I1;
-                ucell.iat2iait(iat1, &I1, &T1);
-                AdjacentAtomInfo adjs;
-                gd.Find_atom(ucell, tau1, T1, I1, &adjs);
-                for (int ad = 0; ad < adjs.adj_num + 1; ++ad)
-                {
-                    const int T2 = adjs.ntype[ad];
-                    const int I2 = adjs.natom[ad];
-                    int iat2 = this->ucell.itia2iat(T2, I2);
-                    if (pmat.get_row_size(iat1) <= 0 || pmat.get_col_size(iat2) <= 0) { continue; }
-                    const ModuleBase::Vector3<int>& R_index = adjs.box[ad];
-                    if (ucell.cal_dtau(iat1, iat2, R_index).norm() * this->ucell.lat0 >= orb_cutoff_[T1] + orb_cutoff_[T2]) { continue; }
-                    hamilt::AtomPair<TR> tmp(iat1, iat2, R_index.x, R_index.y, R_index.z, &pmat);
-                    hR.insert_pair(tmp);
-                }
-            }
-            hR.allocate(nullptr, true);
-            hR.set_paraV(&pmat);
-            if (std::is_same<T, double>::value) { hR.fix_gamma(); }
-        }
-
         void grid_calculation(const int& nbands)const;
 
         //global sizes

source/module_lr/utils/lr_util_hcontainer.h

@@ -44,4 +44,39 @@ namespace LR_Util
         hamilt::HContainer<std::complex<double>>& HR,
         const int& nat,
         const char& type = 'R');
+
+    template<typename T, typename TR>
+    void initialize_HR(hamilt::HContainer<TR>& hR, const UnitCell& ucell, Grid_Driver& gd, const std::vector<double>& orb_cutoff)
+    {
+        const auto& pmat = *hR.get_paraV();
+        for (int iat1 = 0; iat1 < ucell.nat; iat1++)
+        {
+            auto tau1 = ucell.get_tau(iat1);
+            int T1, I1;
+            ucell.iat2iait(iat1, &I1, &T1);
+            AdjacentAtomInfo adjs;
+            gd.Find_atom(ucell, tau1, T1, I1, &adjs);
+            for (int ad = 0; ad < adjs.adj_num + 1; ++ad)
+            {
+                const int T2 = adjs.ntype[ad];
+                const int I2 = adjs.natom[ad];
+                int iat2 = ucell.itia2iat(T2, I2);
+                if (pmat.get_row_size(iat1) <= 0 || pmat.get_col_size(iat2) <= 0) { continue; }
+                const ModuleBase::Vector3<int>& R_index = adjs.box[ad];
+                if (ucell.cal_dtau(iat1, iat2, R_index).norm() * ucell.lat0 >= orb_cutoff[T1] + orb_cutoff[T2]) { continue; }
+                hamilt::AtomPair<TR> tmp(iat1, iat2, R_index.x, R_index.y, R_index.z, &pmat);
+                hR.insert_pair(tmp);
+            }
+        }
+        hR.allocate(nullptr, true);
+        // hR.set_paraV(&pmat);
+        if (std::is_same<T, double>::value) { hR.fix_gamma(); }
+    }
+    template<typename T, typename TR>
+    void initialize_DMR(elecstate::DensityMatrix<T, TR>& dm, const Parallel_Orbitals& pmat, const UnitCell& ucell, Grid_Driver& gd, const std::vector<double>& orb_cutoff)
+    {
+        hamilt::HContainer<TR> hR_tmp(&pmat);
+        initialize_HR<T, TR>(hR_tmp, ucell, gd, orb_cutoff);
+        dm.init_DMR(hR_tmp);
+    }
 }