Feature: LR-TDDFT for open-shell systems

docs/advanced/input_files/input-main.md

@@ -3957,6 +3957,13 @@ Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`.
 - **Description**:  The number of 2-particle states to be solved
 - **Default**: 0
 
+### lr_unrestricted
+- **Type**: Boolean
+- **Description**: Whether to use unrestricted construction for LR-TDDFT (the matrix size will be doubled).
+  - True:  Always use unrestricted LR-TDDFT. 
+  - False: Use unrestricted LR-TDDFT only when the system is open-shell.
+- **Default**: False
+
 ### abs_wavelen_range
 
 - **Type**: Real Real

source/Makefile.Objects

@@ -729,7 +729,6 @@ OBJS_TENSOR=tensor.o\
     operator_lr_exx.o\
     kernel_xc.o\
     pot_hxc_lrtd.o\
-    hsolver_lrtd.o\
     lr_spectrum.o\
     hamilt_casida.o\
     esolver_lrtd_lcao.o\

source/module_basis/module_ao/parallel_2d.h

@@ -15,6 +15,7 @@ class Parallel_2D
     ~Parallel_2D() = default;
 
     Parallel_2D& operator=(Parallel_2D&& rhs) = default;
+    Parallel_2D(Parallel_2D&& rhs) = default;
 
     /// number of local rows
     int get_row_size() const

source/module_esolver/esolver.cpp

@@ -187,13 +187,8 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell)
     else if (esolver_type == "lr_lcao")
     {
         // use constructor rather than Init function to initialize reference (instead of pointers) to ucell
-        if (PARAM.globalv.gamma_only_local){
-            return new LR::ESolver_LR<double, double>(inp, ucell);
-        } else if (PARAM.inp.nspin < 2) {
-            return new LR::ESolver_LR<std::complex<double>, double>(inp, ucell);
-        } else {
-            throw std::runtime_error("LR-TDDFT is not implemented for spin polarized case");
-}
+        if (PARAM.globalv.gamma_only_local) { return new LR::ESolver_LR<double, double>(inp, ucell); }
+        else { return new LR::ESolver_LR<std::complex<double>, double>(inp, ucell); }
     }
     else if (esolver_type == "ksdft_lr_lcao")
     {

source/module_io/read_input_item_tddft.cpp

@@ -327,6 +327,12 @@ void ReadInput::item_lr_tddft()
         read_sync_bool(input.out_wfc_lr);
         this->add_item(item);
     }
+    {
+        Input_Item item("lr_unrestricted");
+        item.annotation = "Whether to use unrestricted construction for LR-TDDFT";
+        read_sync_bool(input.lr_unrestricted);
+        this->add_item(item);
+    }
     {
         Input_Item item("abs_wavelen_range");
         item.annotation = "the range of wavelength(nm) to output the absorption spectrum ";
@@ -337,10 +343,6 @@ void ReadInput::item_lr_tddft()
                 para.input.abs_wavelen_range.push_back(std::stod(item.str_values[i]));
             }
             };
-        item.check_value = [](const Input_Item& item, const Parameter& para) {
-            auto& awr = para.input.abs_wavelen_range;
-            if (awr.size() < 2) { ModuleBase::WARNING_QUIT("ReadInput", "abs_wavelen_range must have two values"); }
-            };
         sync_doublevec(input.abs_wavelen_range, 2, 0.0);
         this->add_item(item);
     }

source/module_io/test/read_input_ptest.cpp

@@ -419,6 +419,7 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_EQ(param.inp.xc_kernel, "LDA");
     EXPECT_EQ(param.inp.lr_solver, "dav");
     EXPECT_DOUBLE_EQ(param.inp.lr_thr, 1e-2);
+    EXPECT_FALSE(param.inp.lr_unrestricted);
     EXPECT_FALSE(param.inp.out_wfc_lr);
     EXPECT_EQ(param.inp.abs_wavelen_range.size(), 2);
     EXPECT_DOUBLE_EQ(param.inp.abs_wavelen_range[0], 0.0);

source/module_lr/AX/AX.h

@@ -13,25 +13,25 @@ namespace LR
         const psi::Psi<double>& c,
         const int& nocc,
         const int& nvirt,
-        psi::Psi<double>& AX_istate);
+        double* AX_istate);
     void cal_AX_blas(
         const std::vector<container::Tensor>& V_istate,
         const psi::Psi<double>& c,
         const int& nocc,
         const int& nvirt,
-        psi::Psi<double>& AX_istate,
+        double* AX_istate,
         const bool add_on = true);
 #ifdef __MPI
     void cal_AX_pblas(
         const std::vector<container::Tensor>& V_istate,
         const Parallel_2D& pmat,
         const psi::Psi<double>& c,
         const Parallel_2D& pc,
-        int naos,
-        int nocc,
-        int nvirt,
-        Parallel_2D& pX,
-        psi::Psi<double>& AX_istate,
+        const int& naos,
+        const int& nocc,
+        const int& nvirt,
+        const Parallel_2D& pX,
+        double* AX_istate,
         const bool add_on=true);
 #endif
     // complex
@@ -40,13 +40,13 @@ namespace LR
         const psi::Psi<std::complex<double>>& c,
         const int& nocc,
         const int& nvirt,
-        psi::Psi<std::complex<double>>& AX_istate);
+        std::complex<double>* AX_istate);
     void cal_AX_blas(
         const std::vector<container::Tensor>& V_istate,
         const psi::Psi<std::complex<double>>& c,
         const int& nocc,
         const int& nvirt,
-        psi::Psi<std::complex<double>>& AX_istate,
+        std::complex<double>* AX_istate,
         const bool add_on = true);
 
 #ifdef __MPI
@@ -55,11 +55,11 @@ namespace LR
         const Parallel_2D& pmat,
         const psi::Psi<std::complex<double>>& c,
         const Parallel_2D& pc,
-        int naos,
-        int nocc,
-        int nvirt,
-        Parallel_2D& pX,
-        psi::Psi<std::complex<double>>& AX_istate,
+        const int& naos,
+        const int& nocc,
+        const int& nvirt,
+        const Parallel_2D& pX,
+        std::complex<double>* AX_istate,
         const bool add_on = true);
 #endif
 }

source/module_lr/AX/AX_parallel.cpp

@@ -14,28 +14,25 @@ namespace LR
         const Parallel_2D& pmat,
         const psi::Psi<double>& c,
         const Parallel_2D& pc,
-        int naos,
-        int nocc,
-        int nvirt,
-        Parallel_2D& pX,
-        psi::Psi<double>& AX_istate,
+        const int& naos,
+        const int& nocc,
+        const int& nvirt,
+        const Parallel_2D& pX,
+        double* AX_istate,
         const bool add_on)
     {
         ModuleBase::TITLE("hamilt_lrtd", "cal_AX_pblas");
-        assert(pmat.comm() == pc.comm());
-        assert(pmat.blacs_ctxt == pc.blacs_ctxt);
-
-        if (pX.comm() != pmat.comm() || pX.blacs_ctxt != pmat.blacs_ctxt)
-            LR_Util::setup_2d_division(pX, pmat.get_block_size(), nvirt, nocc, pmat.blacs_ctxt);
-        else assert(pX.get_local_size() > 0 && AX_istate.get_nbasis() == pX.get_local_size());
+        assert(pmat.comm() == pc.comm() && pmat.comm() == pX.comm());
+        assert(pmat.blacs_ctxt == pc.blacs_ctxt && pmat.blacs_ctxt == pX.blacs_ctxt);
+        assert(pX.get_local_size() > 0);
 
         const int nks = V_istate.size();
 
         Parallel_2D pVc;        // for intermediate Vc
         LR_Util::setup_2d_division(pVc, pmat.get_block_size(), naos, nocc, pmat.blacs_ctxt);
         for (int isk = 0;isk < nks;++isk)
         {
-            AX_istate.fix_k(isk);
+            const int ax_start = isk * pX.get_local_size();
             c.fix_k(isk);
 
             //Vc
@@ -60,7 +57,7 @@ namespace LR
             pdgemm_(&transa, &transb, &nvirt, &nocc, &naos,
                 &alpha, c.get_pointer(), &i1, &ivirt, pc.desc,
                 Vc.data<double>(), &i1, &i1, pVc.desc,
-                &beta, AX_istate.get_pointer(), &i1, &i1, pX.desc);
+                &beta, AX_istate + ax_start, &i1, &i1, pX.desc);
 
         }
     }
@@ -70,28 +67,25 @@ namespace LR
         const Parallel_2D& pmat,
         const psi::Psi<std::complex<double>>& c,
         const Parallel_2D& pc,
-        int naos,
-        int nocc,
-        int nvirt,
-        Parallel_2D& pX,
-        psi::Psi<std::complex<double>>& AX_istate,
+        const int& naos,
+        const int& nocc,
+        const int& nvirt,
+        const Parallel_2D& pX,
+        std::complex<double>* AX_istate,
         const bool add_on)
     {
         ModuleBase::TITLE("hamilt_lrtd", "cal_AX_plas");
-        assert(pmat.comm() == pc.comm());
-        assert(pmat.blacs_ctxt == pc.blacs_ctxt);
-
-        if (pX.comm() != pmat.comm() || pX.blacs_ctxt != pmat.blacs_ctxt)
-            LR_Util::setup_2d_division(pX, pmat.get_block_size(), nvirt, nocc, pmat.blacs_ctxt);
-        else assert(pX.get_local_size() > 0 && AX_istate.get_nbasis() == pX.get_local_size());
+        assert(pmat.comm() == pc.comm() && pmat.comm() == pX.comm());
+        assert(pmat.blacs_ctxt == pc.blacs_ctxt && pmat.blacs_ctxt == pX.blacs_ctxt);
+        assert(pX.get_local_size() > 0);
 
         const int nks = V_istate.size();
 
         Parallel_2D pVc;        // for intermediate Vc
         LR_Util::setup_2d_division(pVc, pmat.get_block_size(), naos, nocc, pmat.blacs_ctxt);
         for (int isk = 0;isk < nks;++isk)
         {
-            AX_istate.fix_k(isk);
+            const int ax_start = isk * pX.get_local_size();
             c.fix_k(isk);
 
             //Vc
@@ -116,7 +110,7 @@ namespace LR
             pzgemm_(&transa, &transb, &nvirt, &nocc, &naos,
                 &alpha, c.get_pointer(), &i1, &ivirt, pc.desc,
                 Vc.data<std::complex<double>>(), &i1, &i1, pVc.desc,
-                &beta, AX_istate.get_pointer(), &i1, &i1, pX.desc);
+                &beta, AX_istate + ax_start, &i1, &i1, pX.desc);
         }
     }
 }

source/module_lr/AX/AX_serial.cpp

@@ -9,18 +9,17 @@ namespace LR
         const psi::Psi<double>& c,
         const int& nocc,
         const int& nvirt,
-        psi::Psi<double>& AX_istate)
+        double* AX_istate)
     {
         ModuleBase::TITLE("hamilt_lrtd", "cal_AX_forloop");
         const int nks = V_istate.size();
         int naos = c.get_nbasis();
-        AX_istate.fix_k(0);
-        ModuleBase::GlobalFunc::ZEROS(AX_istate.get_pointer(), nks * nocc * nvirt);
+        ModuleBase::GlobalFunc::ZEROS(AX_istate, nks * nocc * nvirt);
 
         for (int isk = 0;isk < nks;++isk)
         {
             c.fix_k(isk);
-            AX_istate.fix_k(isk);
+            const int ax_start = isk * nocc * nvirt;
             for (int i = 0;i < nocc;++i)
             {
                 for (int a = 0;a < nvirt;++a)
@@ -29,7 +28,7 @@ namespace LR
                     {
                         for (int mu = 0;mu < naos;++mu)
                         {
-                            AX_istate(i * nvirt + a) += c(nocc + a, mu) * V_istate[isk].data<double>()[nu * naos + mu] * c(i, nu);
+                            AX_istate[ax_start + i * nvirt + a] += c(nocc + a, mu) * V_istate[isk].data<double>()[nu * naos + mu] * c(i, nu);
                         }
                     }
                 }
@@ -41,18 +40,17 @@ namespace LR
         const psi::Psi<std::complex<double>>& c,
         const int& nocc,
         const int& nvirt,
-        psi::Psi<std::complex<double>>& AX_istate)
+        std::complex<double>* AX_istate)
     {
         ModuleBase::TITLE("hamilt_lrtd", "cal_AX_forloop");
         const int nks = V_istate.size();
         int naos = c.get_nbasis();
-        AX_istate.fix_k(0);
-        ModuleBase::GlobalFunc::ZEROS(AX_istate.get_pointer(), nks * nocc * nvirt);
+        ModuleBase::GlobalFunc::ZEROS(AX_istate, nks * nocc * nvirt);
 
         for (int isk = 0;isk < nks;++isk)
         {
             c.fix_k(isk);
-            AX_istate.fix_k(isk);
+            const int ax_start = isk * nocc * nvirt;
             for (int i = 0;i < nocc;++i)
             {
                 for (int a = 0;a < nvirt;++a)
@@ -61,7 +59,7 @@ namespace LR
                     {
                         for (int mu = 0;mu < naos;++mu)
                         {
-                            AX_istate(i * nvirt + a) += std::conj(c(nocc + a, mu)) * V_istate[isk].data<std::complex<double>>()[nu * naos + mu] * c(i, nu);
+                            AX_istate[ax_start + i * nvirt + a] += std::conj(c(nocc + a, mu)) * V_istate[isk].data<std::complex<double>>()[nu * naos + mu] * c(i, nu);
                         }
                     }
                 }
@@ -74,7 +72,7 @@ namespace LR
         const psi::Psi<double>& c,
         const int& nocc,
         const int& nvirt,
-        psi::Psi<double>& AX_istate,
+        double* AX_istate,
         const bool add_on)
     {
         ModuleBase::TITLE("hamilt_lrtd", "cal_AX_blas");
@@ -84,7 +82,7 @@ namespace LR
         for (int isk = 0;isk < nks;++isk)
         {
             c.fix_k(isk);
-            AX_istate.fix_k(isk);
+            const int ax_start = isk * nocc * nvirt;
 
             // Vc[naos*nocc]
             container::Tensor Vc(DAT::DT_DOUBLE, DEV::CpuDevice, { nocc, naos });// (Vc)^T
@@ -101,15 +99,15 @@ namespace LR
             //AX_istate=c^TVc (nvirt major)
             dgemm_(&transa, &transb, &nvirt, &nocc, &naos, &alpha,
                 c.get_pointer(nocc), &naos, Vc.data<double>(), &naos, &beta,
-                AX_istate.get_pointer(), &nvirt);
+                AX_istate + ax_start, &nvirt);
         }
     }
     void cal_AX_blas(
         const std::vector<container::Tensor>& V_istate,
         const psi::Psi<std::complex<double>>& c,
         const int& nocc,
         const int& nvirt,
-        psi::Psi<std::complex<double>>& AX_istate,
+        std::complex<double>* AX_istate,
         const bool add_on)
     {
         ModuleBase::TITLE("hamilt_lrtd", "cal_AX_blas");
@@ -119,7 +117,7 @@ namespace LR
         for (int isk = 0;isk < nks;++isk)
         {
             c.fix_k(isk);
-            AX_istate.fix_k(isk);
+            const int ax_start = isk * nocc * nvirt;
 
             // Vc[naos*nocc] (V is hermitian)
             container::Tensor Vc(DAT::DT_COMPLEX_DOUBLE, DEV::CpuDevice, { nocc, naos });// (Vc)^T
@@ -136,7 +134,7 @@ namespace LR
             //AX_istate=c^\dagger Vc (nvirt major)
             zgemm_(&transa, &transb, &nvirt, &nocc, &naos, &alpha,
                 c.get_pointer(nocc), &naos, Vc.data<std::complex<double>>(), &naos, &beta,
-                AX_istate.get_pointer(), &nvirt);
+                AX_istate + ax_start, &nvirt);
         }
     }
 }