rebase, modify namespace and rename dir

Hamilt and HSolver are not members anymore rename and refactor init_X apply new input
maki49 · Jul 11, 2024 · f01731f · f01731f
1 parent 756af53
commit f01731f
Show file tree

Hide file tree

Showing 62 changed files with 834 additions and 496 deletions.
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -409,12 +409,12 @@
     - [pexsi\_mu\_guard](#pexsi_mu_guard)
     - [pexsi\_elec\_thr](#pexsi_elec_thr)
     - [pexsi\_zero\_thr](#pexsi_zero_thr)
-  - [Beyond DFT](#beyond-dft)
+  - [Linear Response TDDFT](#linear-response-tddft)
     - [xc\_kernel](#xc_kernel)
     - [lr\_solverl](#lr_solver)
     - [lr\_thr](#lr_thr)
     - [nvirt](#nvirt)
-    - [lr_nstates](#lr_nstates)
+    - [lr\_nstates](#lr_nstates)
     - [abs\_wavelen\_range](#abs_wavelen_range)
     - [out\_wfc\_lr](#out_wfc_lr)
 [back to top](#full-list-of-input-keywords)
@@ -3776,7 +3776,7 @@ These variables are used to control the usage of PEXSI (Pole Expansion and Selec
 
 [back to top](#full-list-of-input-keywords)
 
-## Beyond DFT
+## Linear Response TDDFT
 
 These parameters are used to solve the excited states using. e.g. LR-TDDFT.
 

diff --git a/examples/lr-tddft/lcao_H2O/INPUT b/examples/lr-tddft/lcao_H2O/INPUT
@@ -5,7 +5,7 @@ pseudo_dir              ../../../tests/PP_ORB
 orbital_dir                 ../../../tests/PP_ORB
 calculation             scf
 nbands                23
-symmetry               	0
+symmetry               	-1
 
 #Parameters (2.Iteration)
 ecutwfc                 50
@@ -28,7 +28,7 @@ mixing_gg0 0
 lr_nstates 2
 xc_kernel lda
 lr_solver dav
-lr_thr 1e-6
+lr_thr 1e-2
 pw_diag_ndim 2
 
 esolver_type ks-lr

diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt
@@ -16,7 +16,7 @@ add_subdirectory(module_io)
 add_subdirectory(module_relax)
 add_subdirectory(module_ri)
 add_subdirectory(module_parameter)
-add_subdirectory(module_beyonddft)
+add_subdirectory(module_lr)
 
 add_library(
     driver

diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -66,12 +66,12 @@ VPATH=./src_global:\
 ./src_ri:\
 ./module_ri:\
 ./module_parameter:\
-./module_beyonddft:\
-./module_beyonddft/AX:\
-./module_beyonddft/dm_trans:\
-./module_beyonddft/operator_casida:\
-./module_beyonddft/potentials:\
-./module_beyonddft/utils:\
+./module_lr:\
+./module_lr/AX:\
+./module_lr/dm_trans:\
+./module_lr/operator_casida:\
+./module_lr/potentials:\
+./module_lr/utils:\
 ./\
 
 OBJS_ABACUS_PW=${OBJS_MAIN}\

diff --git a/source/driver_run.cpp b/source/driver_run.cpp
@@ -7,6 +7,7 @@
 #include "module_io/print_info.h"
 #include "module_io/winput.h"
 #include "module_md/run_md.h"
+#include "module_parameter/parameter.h"
 
 /**
  * @brief This is the driver function which defines the workflow of ABACUS
@@ -43,8 +44,8 @@ void Driver::driver_run() {
     Check_Atomic_Stru::check_atomic_stru(GlobalC::ucell,
                                          GlobalV::MIN_DIST_COEF);
 
-    //! 2: initialize the ESolver 
-    ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(INPUT, GlobalC::ucell);
+    //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`)
+    ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(INPUT, PARAM.inp, GlobalC::ucell);
 
     //! 3: initialize Esolver and fill json-structure
     p_esolver->before_all_runners(INPUT, GlobalC::ucell);

diff --git a/source/module_beyonddft/lr_spectrum.h b/source/module_beyonddft/lr_spectrum.h
diff --git a/source/module_esolver/esolver.cpp b/source/module_esolver/esolver.cpp
@@ -7,7 +7,7 @@
 #include "esolver_ks_lcaopw.h"
 #include "esolver_ks_lcao.h"
 #include "esolver_ks_lcao_tddft.h"
-#include "module_beyonddft/esolver_lrtd_lcao.h"
+#include "module_lr/esolver_lrtd_lcao.h"
 extern "C"
 {
 #include "module_base/blacs_connector.h"
@@ -122,7 +122,7 @@ std::string determine_type()
 
 
 //Some API to operate E_Solver
-ESolver * init_esolver(Input & input, UnitCell & ucell)
+ESolver* init_esolver(Input& input, const Input_para& input_para, UnitCell& ucell)
 {
 	//determine type of esolver based on INPUT information
 	const std::string esolver_type = determine_type();
@@ -187,9 +187,9 @@ ESolver * init_esolver(Input & input, UnitCell & ucell)
     {
         // use constructor rather than Init function to initialize reference (instead of pointers) to ucell
         if (GlobalV::GAMMA_ONLY_LOCAL)
-            return new ModuleESolver::ESolver_LRTD<double, double>(input, ucell);
-        else if (GlobalV::NSPIN < 4)
-            return new ModuleESolver::ESolver_LRTD<std::complex<double>, double>(input, ucell);
+            return new LR::ESolver_LR<double, double>(input_para, input, ucell);
+        else if (GlobalV::NSPIN < 2)
+            return new LR::ESolver_LR<std::complex<double>, double>(input_para, input, ucell);
         else
             throw std::runtime_error("LR-TDDFT is not implemented for spin polarized case");
     }
@@ -214,25 +214,23 @@ ESolver * init_esolver(Input & input, UnitCell & ucell)
         // force and stress is not needed currently,
         // they will be supported after the analytical gradient
         // of LR-TDDFT is implemented.
-
-        p_esolver->after_all_runners();
-        std::cout << "setting up the esolver for excited state" << std::endl;
+        // after_all_runners() is for output, it is not needed here.
+        std::cout << "Setting up the esolver for excited state..." << std::endl;
         // initialize the 2nd ESolver_LR at the temporary pointer
         ModuleESolver::ESolver* p_esolver_lr = nullptr;
-        if (INPUT.gamma_only)
-            p_esolver_lr = new ModuleESolver::ESolver_LRTD<double, double>(
+        if (GlobalV::GAMMA_ONLY_LOCAL)
+            p_esolver_lr = new LR::ESolver_LR<double, double>(
                 std::move(*dynamic_cast<ModuleESolver::ESolver_KS_LCAO<double, double>*>(p_esolver)),
-                INPUT,
+                input_para,
                 ucell);
         else
-            p_esolver_lr = new ModuleESolver::ESolver_LRTD<std::complex<double>, double>(
+            p_esolver_lr = new LR::ESolver_LR<std::complex<double>, double>(
                 std::move(*dynamic_cast<ModuleESolver::ESolver_KS_LCAO<std::complex<double>, double>*>(p_esolver)),
-                INPUT,
+                input_para,
                 ucell);
         // clean the 1st ESolver_KS and swap the pointer
         ModuleESolver::clean_esolver(p_esolver, false); // do not call Cblacs_exit, remain it for the 2nd ESolver
-        p_esolver = p_esolver_lr;
-        p_esolver_lr = nullptr;
+        return p_esolver_lr;
     }
 #endif
 	else if (esolver_type == "sdft_pw")

diff --git a/source/module_esolver/esolver.h b/source/module_esolver/esolver.h
@@ -4,6 +4,7 @@
 #include "module_base/matrix.h"
 #include "module_cell/unitcell.h"
 #include "module_io/input.h"
+#include "module_parameter/input_parameter.h"
 
 namespace ModuleESolver
 {
@@ -84,7 +85,7 @@ std::string determine_type(void);
  *
  * @return [out] A pointer to an ESolver object that will be initialized.
  */
-ESolver * init_esolver(Input & input, UnitCell & ucell);
+ESolver* init_esolver(Input& input, const Input_para& input_para, UnitCell& ucell);
 
 void clean_esolver(ESolver*& pesolver, const bool lcao_cblacs_exit = false);
 

diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h
@@ -14,11 +14,13 @@
 
 #include <memory>
 
+namespace LR
+{
+    template<typename T, typename TR>
+    class ESolver_LR;
+}
 namespace ModuleESolver
 {
-template<typename T, typename TR>
-class ESolver_LRTD;
-
 template <typename TK, typename TR>
 class ESolver_KS_LCAO : public ESolver_KS<TK> {
   public:
@@ -126,8 +128,8 @@ class ESolver_KS_LCAO : public ESolver_KS<TK> {
     void dpks_cal_projected_DM(
         const elecstate::DensityMatrix<TK, double>* dm) const;
 #endif
-    friend class ESolver_LRTD<double, double>;
-    friend class ESolver_LRTD<std::complex<double>, double>;
+    friend class LR::ESolver_LR<double, double>;
+    friend class LR::ESolver_LR<std::complex<double>, double>;
 };
 } // namespace ModuleESolver
 #endif
diff --git a/source/module_hamilt_lcao/module_gint/gint.h b/source/module_hamilt_lcao/module_gint/gint.h
@@ -92,10 +92,7 @@ class Gint {
     void cpu_force_interface(Gint_inout* inout);
 
     void cpu_force_meta_interface(Gint_inout* inout);
-
-    ///  move operator for the next ESolver to directly use its infomation
-    Gint& operator=(Gint&& rhs);
-
+
     //------------------------------------------------------
     // in gint_vl.cpp
     //------------------------------------------------------

diff --git a/source/module_hamilt_lcao/module_hcontainer/hcontainer.h b/source/module_hamilt_lcao/module_hcontainer/hcontainer.h
@@ -170,11 +170,18 @@ class HContainer
      */
     HContainer(const Parallel_Orbitals* paraV, T* data_pointer = nullptr, const std::vector<int>* ijr_info = nullptr);
 
+    /**
+* @brief set parallel orbital pointer to check parallel information
+*/
     void set_paraV(const Parallel_Orbitals* paraV_in)
     {
         this->paraV = paraV_in;
         for (auto& ap : atom_pairs) ap.set_paraV(paraV_in);
     };
+    /**
+ * @brief get parallel orbital pointer to check parallel information
+ * @return const Parallel_Orbitals* , if return is nullptr, it means HContainer is not in parallel mode
+ */
     const Parallel_Orbitals* get_paraV() const { return this->paraV; };
 
     /**
@@ -426,12 +433,6 @@ class HContainer
         return paraV->get_global_row_size();
     }
 
-    /**
-     * @brief get parallel orbital pointer to check parallel information
-     * @return const Parallel_Orbitals* , if return is nullptr, it means HContainer is not in parallel mode
-     */
-    const Parallel_Orbitals* get_paraV() const{return paraV;}
-
   private:
     // i-j atom pairs, sorted by matrix of (atom_i, atom_j)
     std::vector<AtomPair<T>> atom_pairs;

diff --git a/source/module_io/read_input_item_general.cpp b/source/module_io/read_input_item_general.cpp
@@ -170,12 +170,12 @@ void ReadInput::item_general()
         item.annotation = "the energy solver: ksdft, sdft, ofdft, tddft, lj, dp";
         read_sync_string(esolver_type);
         item.check_value = [](const Input_Item& item, const Parameter& para) {
-            const std::vector<std::string> esolver_types = {"ksdft", "sdft", "ofdft", "tddft", "lj", "dp"};
+            const std::vector<std::string> esolver_types = { "ksdft", "sdft", "ofdft", "tddft", "lj", "dp", "lr", "ks-lr" };
             if (!find_str(esolver_types, para.input.esolver_type))
             {
                 ModuleBase::WARNING_QUIT("ReadInput",
-                                         "esolver_type should be ksdft, sdft, "
-                                         "ofdft, tddft, lj or dp.");
+                    "esolver_type should be ksdft, sdft, "
+                    "ofdft, tddft, lr, ks-lr, lj or dp.");
             }
             if (para.input.esolver_type == "dp")
             {

diff --git a/source/module_io/read_input_item_other.cpp b/source/module_io/read_input_item_other.cpp
@@ -1923,5 +1923,66 @@ void ReadInput::item_others()
         read_sync_double(pexsi_zero_thr);
         this->add_item(item);
     }
+
+    // 25. Linear Response
+    {
+        Input_Item item("lr_nstates");
+        item.annotation = "the number of 2-particle states to be solved";
+        read_sync_int(lr_nstates);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("nvirt");
+        item.annotation = "the number of virtual orbitals to form the 2-particle basis (nocc + nvirt <= nbands)";
+        read_sync_int(nvirt);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("xc_kernel");
+        item.annotation = "exchange correlation (XC) kernel for LR-TDDFT";
+        read_sync_string(xc_kernel);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("lr_solver");
+        item.annotation = "the eigensolver for LR-TDDFT";
+        read_sync_string(lr_solver);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("lr_thr");
+        item.annotation = "convergence threshold of the LR-TDDFT eigensolver";
+        read_sync_double(lr_thr);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("out_wfc_lr");
+        item.annotation = "whether to output the eigenvectors (excitation amplitudes) in the particle-hole basis";
+        read_sync_bool(out_wfc_lr);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("abs_wavelen_range");
+        item.annotation = "the range of wavelength(nm) to output the absorption spectrum ";
+        item.read_value = [](const Input_Item& item, Parameter& para) {
+            size_t count = item.get_size();
+            for (int i = 0; i < count; i++)
+            {
+                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(abs_wavelen_range, 2, 0.0);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("abs_broadening");
+        item.annotation = "the broadening (eta) for LR-TDDFT absorption spectrum";
+        read_sync_double(abs_broadening);
+        this->add_item(item);
+    }
 }
 } // namespace ModuleIO
diff --git a/source/module_beyonddft/AX/AX.h → source/module_lr/AX/AX.h b/source/module_beyonddft/AX/AX.h → source/module_lr/AX/AX.h
@@ -5,7 +5,7 @@
 #ifdef __MPI
 #include "module_basis/module_ao/parallel_2d.h"
 #endif
-namespace hamilt
+namespace LR
 {
     // double
     void  cal_AX_forloop_serial(