Support using hybrid functionals for RT-TDDFT calculations.

docs/advanced/input_files/input-main.md

@@ -4133,11 +4133,12 @@ These variables are used to control berry phase and wannier90 interface paramete
 
 ### out_current
 
-- **Type**: Boolean
+- **Type**: Integer
 - **Description**:
-  - True: Output current.
-  - False: Do not output current.
-- **Default**: False
+  - 0: Do not output current.
+  - 1: Output current using the two-center integral, faster.
+  - 2: Output current using the matrix commutation, more precise.
+- **Default**: 0
 
 ### out_current_k
 

source/Makefile.Objects

@@ -553,6 +553,7 @@ OBJS_IO=input_conv.o\
     write_dipole.o\
     write_init.o\
     td_current_io.o\
+    td_current_io_comm.o\
     write_libxc_r.o\
     output_log.o\
     output_mat_sparse.o\

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -55,7 +55,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::before_all_runners(UnitCell& ucell, cons
     // Run before_all_runners in ESolver_KS_LCAO
     ESolver_KS_LCAO<std::complex<double>, TR>::before_all_runners(ucell, inp);
 
-    td_p = new TD_info(&ucell);
+    td_p = new TD_info(&ucell, this->pv, this->orb_);
     TD_info::td_vel_op = td_p;
     totstep += TD_info::estep_shift;
 
@@ -90,7 +90,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
 
     // Initialize velocity operator for current calculation
-    if (PARAM.inp.td_stype != 1 && TD_info::out_current)
+    if (PARAM.inp.td_stype != 1 && TD_info::out_current == 1)
     {
         // initialize the velocity operator
         velocity_mat = new Velocity_op<TR>(&ucell,
@@ -203,7 +203,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
         }
     }
 
-    if (PARAM.inp.td_stype != 1 && TD_info::out_current)
+    if(PARAM.inp.td_stype != 1 && TD_info::out_current == 1)
     {
         delete velocity_mat;
     }
@@ -296,6 +296,12 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
             srho.begin(is, this->chr, this->pw_rho, ucell.symm);
         }
     }
+#ifdef __EXX
+    if (GlobalC::exx_info.info_ri.real_number)
+        this->exx_nao.exd->exx_hamilt2rho(*this->pelec, this->pv, iter);
+    else
+        this->exx_nao.exc->exx_hamilt2rho(*this->pelec, this->pv, iter);
+#endif
 
     // Calculate delta energy
     this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
@@ -474,6 +480,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::after_scf(UnitCell& ucell, const int ist
     std::cout << " Potential (Ry): " << std::setprecision(15) << this->pelec->f_en.etot << std::endl;
 
     // Output dipole, current, etc.
+    auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, TR>*>(this->p_hamilt);
     ModuleIO::ctrl_output_td<TR>(ucell,
                                  this->chr.rho_save,
                                  this->chr.rhopw,
@@ -485,8 +492,12 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::after_scf(UnitCell& ucell, const int ist
                                  &this->pv,
                                  this->orb_,
                                  this->velocity_mat,
+                                 this->gd,
+                                 hamilt_lcao,
                                  this->RA,
-                                 this->td_p);
+                                 this->td_p,
+                                 this->exx_nao
+                                );
 
     ModuleBase::timer::tick(this->classname, "after_scf");
 }

source/source_io/CMakeLists.txt

@@ -34,6 +34,7 @@ list(APPEND objects
     write_init.cpp
     write_mlkedf_descriptors.cpp
     td_current_io.cpp
+    td_current_io_comm.cpp
     write_libxc_r.cpp
     output_log.cpp
     para_json.cpp

source/source_io/ctrl_output_td.cpp

@@ -20,8 +20,12 @@ void ctrl_output_td(const UnitCell& ucell,
                     const Parallel_Orbitals* pv,
                     const LCAO_Orbitals& orb,
                     const Velocity_op<TR>* velocity_mat,
+                    const Grid_Driver& grid,
+                    hamilt::HamiltLCAO<std::complex<double>, TR>* p_hamilt,
                     Record_adj& RA,
-                    TD_info* td_p)
+                    TD_info* td_p,
+                    const Exx_NAO<std::complex<double>>& exx_nao
+                    )
 {
     ModuleBase::TITLE("ModuleIO", "ctrl_output_td");
 
@@ -46,7 +50,7 @@ void ctrl_output_td(const UnitCell& ucell,
         ModuleBase::WARNING_QUIT("ModuleIO::ctrl_output_td", "Failed to cast ElecState to ElecStateLCAO");
     }
 
-    if (TD_info::out_current)
+    if (TD_info::out_current == 1)
     {
         if (TD_info::out_current_k)
         {
@@ -57,6 +61,10 @@ void ctrl_output_td(const UnitCell& ucell,
             ModuleIO::write_current<TR>(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, RA);
         }
     }
+    else if(TD_info::out_current==2)
+    {
+        ModuleIO::write_current(ucell, grid, istep, psi, pelec, kv, pv, orb, td_p->r_calculator, p_hamilt->getSR(), p_hamilt->getHR(), exx_nao);
+    }
 
     // (3) Output file for restart
     if (PARAM.inp.out_freq_td > 0) // default value of out_freq_td is 0
@@ -88,8 +96,12 @@ template void ctrl_output_td<double>(const UnitCell&,
                                      const Parallel_Orbitals*,
                                      const LCAO_Orbitals&,
                                      const Velocity_op<double>*,
+                                     const Grid_Driver&,
+                                     hamilt::HamiltLCAO<std::complex<double>, double>*,
                                      Record_adj&,
-                                     TD_info*);
+                                     TD_info*,
+                                     const Exx_NAO<std::complex<double>>&
+                                     );
 
 template void ctrl_output_td<std::complex<double>>(const UnitCell&,
                                                    double**,
@@ -102,7 +114,11 @@ template void ctrl_output_td<std::complex<double>>(const UnitCell&,
                                                    const Parallel_Orbitals*,
                                                    const LCAO_Orbitals&,
                                                    const Velocity_op<std::complex<double>>*,
+                                                   const Grid_Driver&,
+                                                   hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>*,
                                                    Record_adj&,
-                                                   TD_info*);
+                                                   TD_info*,
+                                                   const Exx_NAO<std::complex<double>>&
+                                                   );
 
 } // namespace ModuleIO

source/source_io/ctrl_output_td.h

@@ -11,6 +11,11 @@
 #include "source_lcao/module_rt/velocity_op.h"
 #include "source_lcao/record_adj.h"
 #include "source_psi/psi.h"
+#include "source_lcao/hamilt_lcao.h"
+#include "source_lcao/setup_exx.h"
+#ifdef __EXX
+#include <RI/global/Tensor.h>
+#endif
 
 namespace ModuleIO
 {
@@ -27,8 +32,12 @@ void ctrl_output_td(const UnitCell& ucell,
                     const Parallel_Orbitals* pv,
                     const LCAO_Orbitals& orb,
                     const Velocity_op<TR>* velocity_mat,
+                    const Grid_Driver& grid,
+                    hamilt::HamiltLCAO<std::complex<double>, TR>* p_hamilt,
                     Record_adj& RA,
-                    TD_info* td_p);
+                    TD_info* td_p,
+                    const Exx_NAO<std::complex<double>>& exx_nao
+                    );
 
 } // namespace ModuleIO
 

source/source_io/module_parameter/input_parameter.h

@@ -403,7 +403,7 @@ struct Input_para
     int out_wfc_lcao = 0;                 ///< output the wave functions in local basis.
     bool out_dipole = false;              ///< output the dipole or not
     bool out_efield = false;              ///< output the efield or not
-    bool out_current = false;             ///< output the current or not
+    int out_current = 0;                  ///< output the current or not
     bool out_current_k = false;           ///< output tddft current for all k points
     bool out_vecpot = false;              ///< output the vector potential or not
     bool restart_save = false;            ///< restart //Peize Lin add 2020-04-04

source/source_io/read_input_item_exx_dftu.cpp

@@ -164,6 +164,13 @@ void ReadInput::item_exx()
         Input_Item item("exx_separate_loop");
         item.annotation = "if 1, a two-step method is employed, else it will "
                           "start with a GGA-Loop, and then Hybrid-Loop";
+        item.reset_value = [](const Input_Item& item, Parameter& para) {
+            if (para.input.esolver_type == "tddft" && para.input.exx_separate_loop)
+            {
+                GlobalV::ofs_running << "For RT-TDDFT with hybrid functionals, only exx_separate_loop = 0 is supported" << std::endl;
+                para.input.exx_separate_loop = false;
+            }
+        };
         read_sync_bool(input.exx_separate_loop);
         this->add_item(item);
     }

source/source_io/read_input_item_output.cpp

@@ -443,7 +443,7 @@ void ReadInput::item_output()
     {
         Input_Item item("out_current");
         item.annotation = "output current or not";
-        read_sync_bool(input.out_current);
+        read_sync_int(input.out_current);
         this->add_item(item);
     }
     {

source/source_io/td_current_io.h

@@ -6,6 +6,10 @@
 #include "source_estate/module_dm/density_matrix.h"
 #include "source_psi/psi.h"
 #include "source_lcao/module_rt/velocity_op.h"
+#include "source_lcao/setup_exx.h"
+#ifdef __EXX
+#include <RI/global/Tensor.h>
+#endif
 
 namespace ModuleIO
 {
@@ -33,7 +37,22 @@ void write_current(const UnitCell& ucell,
                 const LCAO_Orbitals& orb,
                 const Velocity_op<TR>* cal_current,
                 Record_adj& ra);
-
+/// @brief func to output current calculated using i[r,H] directly
+template <typename TR>
+void write_current(
+    const UnitCell& ucell,
+    const Grid_Driver& GridD,
+    const int istep,
+    const psi::Psi<std::complex<double>>* psi,
+    const elecstate::ElecState* pelec,
+    const K_Vectors& kv,
+    const Parallel_Orbitals* pv,
+    const LCAO_Orbitals& orb,
+    cal_r_overlap_R& r_calculator,
+    const hamilt::HContainer<TR>* sR,
+    const hamilt::HContainer<TR>* hR,
+    const Exx_NAO<std::complex<double>>& exx_nao
+);
 /// @brief calculate sum_n[𝜌_(𝑛𝑘,𝜇𝜈)] for current calculation
 void cal_tmp_DM_k(const UnitCell& ucell,
                 elecstate::DensityMatrix<std::complex<double>, double>& DM_real,
@@ -47,7 +66,61 @@ void cal_tmp_DM(const UnitCell& ucell,
                 elecstate::DensityMatrix<std::complex<double>, double>& DM_real,
                 elecstate::DensityMatrix<std::complex<double>, double>& DM_imag,
                 const int nspin);
+void set_rR_from_hR(const UnitCell& ucell,
+                    const Grid_Driver& GridD,
+                    const LCAO_Orbitals& orb,
+                    const Parallel_Orbitals* pv,
+                    cal_r_overlap_R& r_calculator,
+                    const hamilt::HContainer<std::complex<double>>* hR,
+                    ModuleBase::Vector3<hamilt::HContainer<double>*>& rR);
+template <typename TR>
+void sum_HR(
+    const UnitCell& ucell,
+    const Parallel_Orbitals& pv,
+    const K_Vectors& kv,
+    const hamilt::HContainer<TR>* hR,
+    hamilt::HContainer<std::complex<double>>* full_hR,
+    const Exx_NAO<std::complex<double>>& exx_nao
+);
+
+template <typename Tadd, typename Tfull>
+void add_HR(const hamilt::HContainer<Tadd>* hR, hamilt::HContainer<Tfull>* full_hR);
+
+void init_from_adj(const UnitCell& ucell,
+                   const Grid_Driver& GridD,
+                   const LCAO_Orbitals& orb,
+                   const Parallel_Orbitals* pv,
+                   std::vector<AdjacentAtomInfo>& adjs_all,
+                   ModuleBase::Vector3<hamilt::HContainer<double>*>& rR);
+template <typename TR, typename TA>
+void init_from_hR(const hamilt::HContainer<TR>* hR, hamilt::HContainer<TA>* aimR);
+template <typename TR>
+void cal_velocity_basis_k(const UnitCell& ucell,
+                          const LCAO_Orbitals& orb,
+                          const Parallel_Orbitals* pv,
+                          const K_Vectors& kv,
+                          const ModuleBase::Vector3<hamilt::HContainer<double>*>& rR,
+                          const hamilt::HContainer<TR>& sR,
+                          const hamilt::HContainer<std::complex<double>>& hR,
+                          std::vector<ModuleBase::Vector3<std::complex<double>*>>& velocity_basis_k);
 
+void cal_velocity_matrix(const psi::Psi<std::complex<double>>* psi,
+                         const Parallel_Orbitals* pv,
+                         const K_Vectors& kv,
+                         const std::vector<ModuleBase::Vector3<std::complex<double>*>>& velocity_basis_k,
+                         std::vector<std::array<ModuleBase::ComplexMatrix, 3>>& velocity_k);
+template <typename TR>
+void cal_current_comm_k(const UnitCell& ucell,
+                        const Grid_Driver& GridD,
+                        const LCAO_Orbitals& orb,
+                        const Parallel_Orbitals* pv,
+                        const K_Vectors& kv,
+                        cal_r_overlap_R& r_calculator,
+                        const hamilt::HContainer<TR>& sR,
+                        const hamilt::HContainer<std::complex<double>>& hR,
+                        const psi::Psi<std::complex<double>>* psi,
+                        const elecstate::ElecState* pelec,
+                        std::vector<ModuleBase::Vector3<double>>& current_k);
 #endif // __LCAO
 } // namespace ModuleIO
 #endif