from-scratch LR

maki49 · Nov 11, 2023 · e7eb223 · e7eb223
1 parent dccd53a
commit e7eb223
Show file tree

Hide file tree

Showing 17 changed files with 438 additions and 114 deletions.
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -363,7 +363,6 @@
     - [sccut](#sccut)
     - [sc\_file](#sc_file)
   - [Beyond DFT](#beyond-dft)
-    - [beyonddft\_method](#beyonddft_method)
     - [nstates](#nstates)
     - [xc\_kernel](#xc_kernel)
     - [lr\_solverl](#lr_solver)
@@ -413,6 +412,8 @@ These variables are used to control general system parameters.
   - tddft: real-time time-dependent density functional theory (TDDFT)
   - lj: Leonard Jones potential
   - dp: DeeP potential, see details in [md.md](../md.md#dpmd)
+  - ks-lr: Kohn-Sham density functional theory + LR-TDDFT
+  - lr: LR-TDDFT with given KS orbitals
 - **Default**: ksdft
 
 ### symmetry
@@ -3338,12 +3339,6 @@ These variables are used to control the usage of deltaspin functionality.
 
 These parameters are used to solve the excited states using. e.g. lr-tddft
 
-### beyonddft_method
-
-- **Type**: String
-- **Description**: The method to solve the excited states. Currently, only `lr-tddft` is supported.
-- **Default**: none
-
 ### nstates
 
 - **Type**: Integer

diff --git a/source/driver_run.cpp b/source/driver_run.cpp
@@ -51,7 +51,18 @@ void Driver::driver_run()
     }
 
     // 4. Initialize Esolver
-    p_esolver->Init(INPUT, GlobalC::ucell);
+#ifdef __LCAO
+    if (GlobalV::ESOLVER_TYPE == "lr")
+        // use constructor rather than Init function to initialize reference (instead of pointers) to ucell
+        if (GlobalV::GAMMA_ONLY_LOCAL)
+            p_esolver = new ModuleESolver::ESolver_LRTD<double, double>(INPUT, GlobalC::ucell);
+        else if (GlobalV::NSPIN < 4)
+            p_esolver = new ModuleESolver::ESolver_LRTD<std::complex<double>, double>(INPUT, GlobalC::ucell);
+        else
+            throw std::runtime_error("LR-TDDFT is not implemented for spin polarized case");
+    else
+#endif
+        p_esolver->Init(INPUT, GlobalC::ucell);
 
     //------------------------------------------------------------
     // This part onward needs to be refactored.
@@ -80,7 +91,7 @@ void Driver::driver_run()
 
 #ifdef __LCAO
     //---------beyond DFT: set up the next ESolver---------
-    if (INPUT.beyonddft_method == "lr-tddft")
+    if (INPUT.esolver_type == "ks-lr")
     {
         std::cout << "setting up the esolver for excited state" << std::endl;
         // ModuleESolver::ESolver_KS_LCAO* p_esolver_lcao_tmp = dynamic_cast<ModuleESolver::ESolver_KS_LCAO<double, double>*>(p_esolver);
@@ -89,11 +100,7 @@ void Driver::driver_run()
             p_esolver_lr = new ModuleESolver::ESolver_LRTD<double, double>(std::move(*dynamic_cast<ModuleESolver::ESolver_KS_LCAO<double, double>*>(p_esolver)), INPUT, GlobalC::ucell);
         else
             p_esolver_lr = new ModuleESolver::ESolver_LRTD<std::complex<double>, double>(std::move(*dynamic_cast<ModuleESolver::ESolver_KS_LCAO<std::complex<double>, double>*>(p_esolver)), INPUT, GlobalC::ucell);
-
-        std::cout << "before clean ks" << std::endl;
         ModuleESolver::clean_esolver(p_esolver);
-        std::cout << "after clean ks" << std::endl;
-
         p_esolver_lr->Run(0, GlobalC::ucell);
 
         std::cout << "before clean lr" << std::endl;
@@ -105,7 +112,6 @@ void Driver::driver_run()
 
     if (INPUT.basis_type == "lcao")
         Cblacs_exit(1); // clean up blacs after all the esolvers are cleaned up without closing MPI
-    std::cout << "befor end" << std::endl;
 #else
     ModuleESolver::clean_esolver(p_esolver);
 #endif

diff --git a/source/module_beyonddft/esolver_lrtd_lcao.h b/source/module_beyonddft/esolver_lrtd_lcao.h
@@ -10,7 +10,6 @@
 #include <memory>
 
 #include "module_esolver/esolver_ks_lcao.h" //for the move constructor
-#include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h"
 #include "module_hamilt_lcao/module_gint/gint_gamma.h"
 #include "module_hamilt_lcao/module_gint/gint_k.h"
 #include "module_hamilt_lcao/module_gint/grid_technique.h"
@@ -39,7 +38,8 @@ namespace ModuleESolver
     public:
         /// @brief  a move constructor from ESolver_KS_LCAO
         ESolver_LRTD(ModuleESolver::ESolver_KS_LCAO<T, TR>&& ks_sol, Input& inp, UnitCell& ucell);
-        ESolver_LRTD() {}
+        /// @brief a from-scratch constructor
+        ESolver_LRTD(Input& inp, UnitCell& ucell);
         ~ESolver_LRTD() {
             delete this->p_hamilt;
             delete this->phsol;
@@ -64,7 +64,6 @@ namespace ModuleESolver
     protected:
 
         const UnitCell& ucell;
-        const Input& input;
 
         hamilt::Hamilt<T>* p_hamilt = nullptr;  //opsd problem first to use base calss
         hsolver::HSolver<T>* phsol = nullptr;
@@ -95,19 +94,14 @@ namespace ModuleESolver
         int nsk = 1; //nspin*nks
         int nspin = 1;
 
-        // basis info (currently use GlobalC)
-        // LCAO_Orbitals orb;
-        // adj info 
-        Record_adj ra;
-        // grid parallel info (no need for 2d-block distribution)?
-        // Grid_Technique gridt;
-        // grid integration method(will be moved to OperatorKernelHxc)
-        ModulePW::PW_Basis_Big bigpw;
         Grid_Technique gt;
         Gint_Gamma gint_g;
         Gint_k gint_k;
         typename TGint<T>::type* gint = nullptr;
 
+        std::string xc_kernel;
+        std::string lr_solver;
+
         void set_gint();
 
         K_Vectors kv;
@@ -119,14 +113,24 @@ namespace ModuleESolver
         /// @brief variables for parallel distribution of matrix in AO representation
         Parallel_Orbitals paraMat_;
 
-        /// move to hsolver::updatePsiK
-        void init_X();
+        /// @brief allocate and initialize X
+        void init_X(const int& nvirt_input);
+        /// @brief allocate and initialize A matrix, density matrix and eignensolver
+        void init_A(hamilt::HContainer<double>* pHR_in, const double lr_thr);
+        /// @brief read in the ground state wave function, band energy and occupation
+        void read_ks_wfc();
+        /// @brief  read in the ground state charge density
+        void read_ks_chg(Charge& chg);
+
+        void init_pot(const Charge& chg_gs);
 
 #ifdef __EXX
         /// Tdata of Exx_LRI is same as T, for the reason, see operator_lr_exx.h
         std::shared_ptr<Exx_LRI<T>> exx_lri = nullptr;
         void move_exx_lri(std::shared_ptr<Exx_LRI<double>>&);
         void move_exx_lri(std::shared_ptr<Exx_LRI<std::complex<double>>>&);
+
+        std::unique_ptr<TwoCenterBundle> two_center_bundle;
 #endif
 
     };