GGA-kernel & refactor & example-Si

unify nks

respective gint for real and imag part in OperatorHxc

refactor print and modify init

scal with 1/Nk and transpose dm_trans

imag part in absorption spectra

hse kernel & IO eigenstate & spectrum-only

fix after rebase

remove a non-used parameter in

fix after rebase 2

fix dipole magnitude

fix k-first in dav_subspace (still bug)

fix V in AX

xc_kernel: replace matrix with vector

refactor XC-kernel and fix the missing e2 (GGA still bug)

fix invalid lr_ethr

fix GGA kernel

make the broadening eta an input parameter

fix after rebase

fix after rebase (dav_subspace still negative eig)

fix exx compile and Makefile

default move Parallel_2D

move LR-init from driver to esolver; rm some move constructors

CMakeLists.txt

@@ -714,7 +714,8 @@ if(ENABLE_LCAO)
     dftu
     hcontainer
     deltaspin
-    numerical_atomic_orbitals)
+    numerical_atomic_orbitals
+    lr)
   if(USE_ELPA)
     target_link_libraries(${ABACUS_BIN_NAME} genelpa)
   endif()

docs/advanced/input_files/input-main.md

@@ -410,9 +410,13 @@
     - [pexsi\_elec\_thr](#pexsi_elec_thr)
     - [pexsi\_zero\_thr](#pexsi_zero_thr)
   - [Beyond DFT](#beyond-dft)
-    - [nstates](#nstates)
     - [xc\_kernel](#xc_kernel)
     - [lr\_solverl](#lr_solver)
+    - [lr\_thr](#lr_thr)
+    - [nvirt](#nvirt)
+    - [lr_nstates](#lr_nstates)
+    - [abs\_wavelen\_range](#abs_wavelen_range)
+    - [out\_wfc\_lr](#out_wfc_lr)
 [back to top](#full-list-of-input-keywords)
 
 ## System variables
@@ -3774,39 +3778,60 @@ These variables are used to control the usage of PEXSI (Pole Expansion and Selec
 
 ## Beyond DFT
 
-These parameters are used to solve the excited states using. e.g. lr-tddft
+These parameters are used to solve the excited states using. e.g. LR-TDDFT.
 
-### nstates
-
-- **Type**: Integer
-- **Description**:  The number of 2-particle states to be solved
-- **Default**: 0
-
-<<<<<<< HEAD
-<<<<<<< HEAD
-[back to top](#full-list-of-input-keywords)
-<<<<<<< HEAD
-=======
->>>>>>> 7d4fe3b32 (add operator=(T&&) in Grid_Technique and its 4 base classes)
-=======
-=======
->>>>>>> 9a97eada1 (psi wrapperL: k1<->bfirst)
 ### xc_kernel
 
 - **Type**: String
-- **Description**: The exchange-correlation kernel used in the calculation. Currently, only `RPA` and `LDA` is supported.
+- **Description**: The exchange-correlation kernel used in the calculation. 
+Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`.
 - **Default**: LDA
 
 ### lr_solver
 
 - **Type**: String
-- **Description**: The diagonalization method for LR-TDDFT.
-- **Default**: LDA
+- **Description**: The method to solve the Casida equation $AX=\Omega X$ in LR-TDDFT under Tamm-Dancoff approximation (TDA), where $A_{ai,bj}=(\epsilon_a-\epsilon_i)\delta_{ij}\delta_{ab}+(ai|f_{Hxc}|bj)+\alpha_{EX}(ab|ij)$ is the particle-hole excitation matrix and $X$ is the transition amplitude.
+  - `dav`: Construct $AX$ and diagonalize the Hamiltonian matrix iteratively with Davidson algorithm.
+  - `lapack`: Construct the full $A$ matrix and directly diagonalize with LAPACK.
+  - `spectrum`: Calculate absorption spectrum only without solving Casida equation. The `OUT.${suffix}/` directory should contain the
+  files for LR-TDDFT eigenstates and eigenvalues, i.e. `Excitation_Energy.dat` and `Excitation_Amplitude_${processor_rank}.dat`
+   output by setting `out_wfc_lr` to true.
+- **Default**: dav
+
+### lr_thr
+
+- **Type**: Real
+- **Description**: The convergence threshold of iterative diagonalization solver fo LR-TDDFT. It is a pure-math number with the same as [pw_diag_thr](#pw_diag_thr), but since the Casida equation is a one-shot eigenvalue problem, it is also the convergence threshold of LR-TDDFT.
+- **Default**: 1e-2
+
+### nvirt
+
+- **Type**: Integer
+- **Description**: The number of virtual orbitals used in the LR-TDDFT calculation.
+- **Default**: 1
+
+### lr_nstates
+
+- **Type**: Integer
+- **Description**:  The number of 2-particle states to be solved
+- **Default**: 0
+
+### abs_wavelen_range
+
+- **Type**: Real Real
+- **Description**: The range of the wavelength for the absorption spectrum calculation.
+- **Default**: 0.0 0.0
+
+### out_wfc_lr
+
+- **Type**: Boolean
+- **Description**: Whether to output the eigenstates (excitation energy) and eigenvectors (excitation amplitude) of the LR-TDDFT calculation.
+The output files are `OUT.${suffix}/Excitation_Energy.dat` and `OUT.${suffix}/Excitation_Amplitude_${processor_rank}.dat`.
+- **Default**: False
+
+### abs_broadening
+- **Type**: Real
+- **Description**: The broadening factor $\eta$ for the absorption spectrum calculation.
+- **Default**: 0.01
 
-<<<<<<< HEAD
-[back to top](#full-list-of-input-keywords)
->>>>>>> ab3506693 (Framework: kernel, pot, operator, hamilt)
->>>>>>> 8316c5c67 (Framework: kernel, pot, operator, hamilt)
-=======
 [back to top](#full-list-of-input-keywords)
->>>>>>> 9a97eada1 (psi wrapperL: k1<->bfirst)

examples/lr-tddft/lcao_H2O/INPUT

@@ -23,8 +23,9 @@ smearing_sigma                   0.02
 #Parameters (5.Mixing)
 mixing_type             pulay
 mixing_beta             0.4
+mixing_gg0 0
 
-nstates 2
+lr_nstates 2
 xc_kernel lda
 lr_solver dav
 lr_thr 1e-6
@@ -37,3 +38,4 @@ out_alllog	1
 
 nvirt 19
 abs_wavelen_range  40 180
+abs_broadening 0.01

examples/lr-tddft/lcao_H2O/STRU

@@ -6,9 +6,6 @@ NUMERICAL_ORBITAL
 H_gga_8au_60Ry_2s1p.orb
 O_gga_7au_60Ry_2s2p1d.orb
 
-NUMERICAL_DESCRIPTOR
-jle.orb
-
 
 LATTICE_CONSTANT
 1

examples/lr-tddft/lcao_Si2/INPUT

@@ -0,0 +1,40 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix                  abacus
+pseudo_dir              ../../../tests/PP_ORB
+orbital_dir                 ../../../tests/PP_ORB
+calculation             scf
+nbands                23
+symmetry               	0
+
+#Parameters (2.Iteration)
+ecutwfc                 50
+scf_thr                     1e-6
+scf_nmax                   100
+
+#Parameters (3.Basis)
+basis_type              lcao 
+gamma_only              0
+
+#Parameters (4.Smearing)
+smearing_method                gaussian
+smearing_sigma                   0.02
+
+#Parameters (5.Mixing)
+mixing_type             pulay
+mixing_beta             0.4
+
+lr_nstates 10 # for test/debug, you can try a smaller one like 2
+xc_kernel lda
+lr_solver dav
+lr_thr 1e-2
+pw_diag_ndim 4
+
+esolver_type ks-lr
+out_alllog	1
+#out_wfc_lcao 1
+#out_chg 1
+
+nvirt 19
+abs_wavelen_range  100 175
+#diago_full_acc 1

examples/lr-tddft/lcao_Si2/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+2 2 2 0 0 0

examples/lr-tddft/lcao_Si2/STRU

@@ -0,0 +1,24 @@
+#This is the atom file containing all the information
+#about the lattice structure.
+
+ATOMIC_SPECIES
+Si 1.000 Si.pz-vbc.UPF 	#Element, Mass, Pseudopotential
+
+NUMERICAL_ORBITAL
+./Si_lda_8.0au_50Ry_2s2p1d
+
+LATTICE_CONSTANT
+10.2  			#Lattice constant
+
+LATTICE_VECTORS
+0.5 0.5 0.0 		#Lattice vector 1
+0.5 0.0 0.5 		#Lattice vector 2
+0.0 0.5 0.5 		#Lattice vector 3
+
+ATOMIC_POSITIONS
+Cartesian 		#Cartesian(Unit is LATTICE_CONSTANT)
+Si 			#Name of element	
+0.0			#Magnetic for this element.
+2			#Number of atoms
+0.00 0.00 0.00 0 0 0	#x,y,z, move_x, move_y, move_z
+0.25 0.25 0.25 1 1 1

source/Makefile.Objects

@@ -66,6 +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:\
 ./\
 
 OBJS_ABACUS_PW=${OBJS_MAIN}\
@@ -106,6 +112,7 @@ ${OBJS_DFTU}\
 ${OBJS_DELTASPIN}\
 ${OBJS_TENSOR}\
 ${OBJS_HSOLVER_PEXSI}\
+${OBJS_LR}\
 
 OBJS_MAIN=main.o\
     driver.o\
@@ -664,3 +671,19 @@ OBJS_TENSOR=tensor.o\
     linalg_op.o\
     cpu_allocator.o\
     refcount.o
+
+    OBJS_LR=lr_util.o\
+    lr_util_hcontainer.o\
+    AX_parallel.o\
+    AX_serial.o\
+    dm_trans_parallel.o\
+    dm_trans_serial.o\
+    dmr_complex.o\
+    operator_lr_hxc.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/driver_run.cpp

@@ -8,13 +8,6 @@
 #include "module_io/winput.h"
 #include "module_md/run_md.h"
 
-#ifdef __LCAO
-#include "module_beyonddft/esolver_lrtd_lcao.hpp"
-#endif
-extern "C"
-{
-#include "module_base/blacs_connector.h"
-}
 /**
  * @brief This is the driver function which defines the workflow of ABACUS
  * calculations. It relies on the class Esolver, which is a class that organizes
@@ -34,10 +27,7 @@ void Driver::driver_run() {
     ModuleBase::TITLE("Driver", "driver_line");
     ModuleBase::timer::tick("Driver", "driver_line");
 
-    //! 1: initialize the ESolver 
-    ModuleESolver::ESolver *p_esolver = ModuleESolver::init_esolver();
-
-    //! 2: setup cell and atom information
+    //! 1: setup cell and atom information
 
     // this warning should not be here, mohan 2024-05-22
 #ifndef __LCAO
@@ -53,18 +43,10 @@ void Driver::driver_run() {
     Check_Atomic_Stru::check_atomic_stru(GlobalC::ucell,
                                          GlobalV::MIN_DIST_COEF);
 
-    //! 3: initialize Esolver and fill json-structure 
-#ifdef __LCAO
-    if (GlobalV::ESOLVER_TYPE == "lr")s
-        // 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
+    //! 2: initialize the ESolver 
+    ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(INPUT, GlobalC::ucell);
+
+    //! 3: initialize Esolver and fill json-structure
     p_esolver->before_all_runners(INPUT, GlobalC::ucell);
 
     // this Json part should be moved to before_all_runners, mohan 2024-05-12
@@ -96,33 +78,8 @@ void Driver::driver_run() {
     //! 5: clean up esolver
     p_esolver->after_all_runners();
 
-#ifdef __LCAO
-    //---------beyond DFT: set up the next ESolver---------
-    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);
-        ModuleESolver::ESolver* p_esolver_lr = nullptr;
-        if (INPUT.gamma_only)
-            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);
-        ModuleESolver::clean_esolver(p_esolver);
-        p_esolver_lr->Run(0, GlobalC::ucell);
-        p_esolver_lr->postprocess();
-
-        std::cout << "before clean lr" << std::endl;
-        ModuleESolver::clean_esolver(p_esolver_lr);
-        std::cout << "after clean lr" << std::endl;
-    } //----------------------beyond DFT------------------------
-    else
-        ModuleESolver::clean_esolver(p_esolver);
-
-    if (INPUT.basis_type == "lcao")
-        Cblacs_exit(1); // clean up blacs after all the esolvers are cleaned up without closing MPI
-#else
     ModuleESolver::clean_esolver(p_esolver);
-#endif
+
     ModuleBase::timer::tick("Driver", "driver_line");
     return;
 }

source/module_basis/module_ao/parallel_2d.cpp

@@ -136,28 +136,3 @@ void Parallel_2D::set_serial(const int mg, const int ng)
     blacs_ctxt = -1;
 #endif
 }
-
-Parallel_2D& Parallel_2D::operator=(Parallel_2D&& rhs)
-{
-    this->nrow = rhs.nrow;
-    this->ncol = rhs.ncol;
-    this->nloc = rhs.nloc;
-    this->nb = rhs.nb;
-    this->dim0 = rhs.dim0;
-    this->dim1 = rhs.dim1;
-    this->coord[0] = rhs.coord[0];
-    this->coord[1] = rhs.coord[1];
-    this->testpb = rhs.testpb;
-    this->local2global_row_ = std::move(rhs.local2global_row_);
-    this->local2global_col_ = std::move(rhs.local2global_col_);
-    this->global2local_row_ = std::move(rhs.global2local_row_);
-    this->global2local_col_ = std::move(rhs.global2local_col_);
-
-#ifdef __MPI
-    this->blacs_ctxt = rhs.blacs_ctxt;
-    this->comm_2D = rhs.comm_2D;
-    for (int i = 0; i < 9; ++i)
-        this->desc[i] = rhs.desc[i];
-#endif
-    return *this;
-}

source/module_basis/module_ao/parallel_2d.h

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