feat: add BSE module

docs/advanced/input_files/input-main.md

@@ -203,6 +203,8 @@
     - [out\_element\_info](#out_element_info)
     - [restart\_save](#restart_save)
     - [rpa](#rpa)
+    - [rpa_out_vel](#rpa_out_vel)
+    - [rpa_outdir](#rpa_outdir)
     - [out\_pchg](#out_pchg)
     - [out\_wfc\_norm](#out_wfc_norm)
     - [out\_wfc\_re\_im](#out_wfc_re_im)
@@ -542,9 +544,6 @@
     - [pexsi\_elec\_thr](#pexsi_elec_thr)
     - [pexsi\_zero\_thr](#pexsi_zero_thr)
   - [Linear Response TDDFT](#linear-response-tddft)
-    - [ri\_hartree\_benchmark](#ri_hartree_benchmark)
-    - [aims\_nbasis](#aims_nbasis)
-  - [Linear Response TDDFT (Under Development Feature)](#linear-response-tddft-under-development-feature)
     - [xc\_kernel](#xc_kernel)
     - [lr\_init\_xc\_kernel](#lr_init_xc_kernel)
     - [lr\_solver](#lr_solver)
@@ -557,6 +556,23 @@
     - [out\_wfc\_lr](#out_wfc_lr)
     - [abs\_gauge](#abs_gauge)
     - [abs\_broadening](#abs_broadening)
+    - [bse\_tda](#bse_tda)
+    - [bse\_spin\_types](#bse_spin_types)
+    - [bse\_mem\_save](#bse_mem_save)
+    - [bse\_ri\_hartree](#bse_ri_hartree)
+    - [bse\_use\_fine\_kgrid](#bse_use_fine_kgrid)
+    - [out\_bse\_ab](#out_bse_ab)
+    - [bse\_continue](#bse_continue)
+    - [plot\_istate](#plot_istate)
+    - [exciton\_plot\_type](#exciton_plot_type)
+    - [exciton\_plot\_format](#exciton_plot_format)
+    - [exciton\_fixed\_coordinate](#exciton_fixed_coordinate)
+    - [exciton\_slice\_plane](#exciton_slice_plane)
+    - [exciton\_slice\_pos](#exciton_slice_pos)
+    - [exciton\_slice\_npoints](#exciton_slice_npoints)
+    - [exciton\_slice\_scale](#exciton_slice_scale)
+    - [ri\_hartree\_benchmark](#ri_hartree_benchmark)
+    - [aims\_nbasis](#aims_nbasis)
   - [Reduced Density Matrix Functional Theory](#reduced-density-matrix-functional-theory)
     - [rdmft](#rdmft)
     - [rdmft\_power\_alpha](#rdmft_power_alpha)
@@ -2334,11 +2350,23 @@
 ### rpa
 
 - **Type**: Boolean
-- **Description**: Generate output files used in rpa calculations.
+- **Description**: Generate output files used in rpa calculations for LibRPA package.
 
   > Note: If symmetry is set to 1, additional files containing the necessary information for exploiting symmetry in the subsequent rpa calculation will be output: irreducible_sector.txt, symrot_k.txt and symrot_R.txt.
 - **Default**: False
 
+### rpa_out_vel
+
+- **Type**: Boolean
+- **Description**: Whether to output velocity matrix for LibRPA package.
+- **Default**: False
+
+### rpa_outdir
+
+- **Type**: String
+- **Description**: Output directory for LibRPA package.
+- **Default**: "./OUT_rpa/"
+
 ### out_pchg
 
 - **Type**: String
@@ -4798,27 +4826,13 @@
 
 ## Linear Response TDDFT
 
-### ri_hartree_benchmark
-
-- **Type**: String
-- **Description**: Whether to use the RI approximation for the Hartree term in LR-TDDFT for benchmark (with FHI-aims/ABACUS read-in style)
-- **Default**: none
-
-### aims_nbasis
-
-- **Type**: A number(ntype) of Integers
-- **Availability**: *ri_hartree_benchmark = aims*
-- **Description**: Atomic basis set size for each atom type (with the same order as in STRU) in FHI-aims.
-- **Default**: {} (empty list, where ABACUS use its own basis set size)
-
-[back to top](#full-list-of-input-keywords)
-
-## Linear Response TDDFT (Under Development Feature)
+These parameters are used to solve the excited states using. e.g. LR-TDDFT or Bethe-Salpeter Equation.
 
 ### xc_kernel
 
 - **Type**: String
-- **Description**: The exchange-correlation kernel used in the calculation. Currently supported: RPA, LDA, PBE, HSE, HF.
+- **Description**: The exchange-correlation kernel used in the calculation. 
+Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`, `BSE`.
 - **Default**: LDA
 
 ### lr_init_xc_kernel
@@ -4833,10 +4847,12 @@
 ### lr_solver
 
 - **Type**: String
-- **Description**: The method to solve the Casida equation in LR-TDDFT under Tamm-Dancoff approximation (TDA).
-  - dav/dav_subspace/cg: Construct and diagonalize the Hamiltonian matrix iteratively with Davidson/Non-ortho-Davidson/CG algorithm.
-  - lapack: Construct the full matrix and directly diagonalize with LAPACK.
-  - spectrum: Calculate absorption spectrum only without solving Casida equation.
+- **Description**: The method to solve the Casida equation in LR-TDDFT.
+  - `dav`/`dav_subspace`/ `cg`: Construct and diagonalize the Hamiltonian matrix iteratively with Davidson/Non-ortho-Davidson/CG algorithm. Tamm-Dancoff approximation (TDA) only.
+  - `lapack`: Construct the matrix and directly diagonalize with LAPACK. TDA only.
+  - `elpa`: Construct the resonant (and anti-resonant) matrix and diagonalize with ELPA.
+  - `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.
+  - `plot`: Plot the exciton wave function, should identify `plot_istate`.
 - **Default**: dav
 
 ### lr_thr
@@ -4888,19 +4904,150 @@
 ### abs_gauge
 
 - **Type**: String
-- **Description**: Whether to use length or velocity gauge to calculate the absorption spectrum in LR-TDDFT.
-- **Default**: length
+- **Description**: Whether to use `velocity` or `length` formulation to calculate the absorption spectrum. `length` is only suitable for non-periodic system.
+- **Default**: velocity
 
 ### abs_broadening
 
 - **Type**: Real
 - **Description**: The broadening factor for the absorption spectrum calculation.
 - **Default**: 0.01
 
+
+### bse_tda
+
+- **Type**: String
+- **Description**: Whether Tamm-Dancoff Approximation is used (can be 'tda', 'full' or 'both').
+- **Default**: tda
+
+### bse_spin_types
+
+- **Type**: Vector of String
+- **Description**: spin types for close-shell case to be calculated in one task (can be 'singlet', 'triplet', and for test 'rpa', 'ipa').
+- **Defalut**: \{singlet, triplet\}
+
+### bse_mem_save
+
+- **Type**: Boolean
+- **Description**: Whether to save memory by adding V and W to BSE matrix directly. This option is useful when out-of-memory occurs. If this option is on, `bse_ri_hartree` will be on and `bse_continue` will be off automatically.
+- **Default**: false
+
+### bse_ri_hartree
+
+- **Type**: Boolean
+- **Description**: Whether to use RI approximation for Hartree term in BSE.
+- **Default**: true
+
+### bse_use_fine_kgrid
+
+- **Type**: Boolean
+- **Description**: Whether to use a finer k-grid for BSE. If you want to turn it on, file `band_kpath_info`, `band_KS_eigenvector_k_{index}.txt`, `KS_band_spin_{index}.txt` and `GW_band_spin_{index}.txt` should be prepared.
+- **Default**: false
+
+### out_bse_ab
+
+- **Type**: Boolean
+- **Description**: Whether to output the AB matrix to file.
+- **Default**: false
+
+### bse_continue
+
+- **Type**: Integer
+- **Description**: Which step to continue from previous BSE calculation.
+  - 0: new;
+  - 1: continue from A_V;
+  - 2: continue from A_V and A_W;
+  - 3: continue from A_V, A_W and B_V;
+  - 4: continue from A_V, A_W, B_V and B_W
+- **Default**: 0
+
+### plot_istate
+
+- **Type**: Integer
+- **Description**: The index of excited state to be plotted (starting from 0)
+- **Default**: 0
+
+### exciton_plot_type
+
+- **Type**: String
+- **Description**: The exciton density represented when `lr_solver` is `plot`.
+  - `average`: Integrate out the other particle coordinate and plot the average electron and hole densities.
+  - `conditional`: Fix one particle at the Cartesian position specified by the exciton fixed-coordinate parameters and plot the density of the other particle.
+- **Default**: average
+
+### exciton_plot_format
+
+- **Type**: String
+- **Description**: The output format for exciton-density plotting.
+  - `cube`: Write three-dimensional density files in Gaussian cube format.
+  - `slice`: Write two-dimensional cross-section data files.
+  - `both`: Write both cube and slice files.
+  - The default value is `cube`.
+- **Default**: cube
+
+### exciton_fixed_coordinate
+
+- **Type**: Vector of Real (6 values)
+- **Description**: Cartesian coordinates in Bohr used by `exciton_plot_type = conditional`, in the order `hole_x hole_y hole_z electron_x electron_y electron_z`.
+- **Default**: 0.0 0.0 0.0 0.0 0.0 0.0
+- **Unit**: Bohr
+
+### exciton_slice_plane
+
+- **Type**: String
+- **Description**: Pair of lattice-vector directions spanning the cross section.
+  - `ab`: Plane spanned by lattice vectors **a** and **b**.
+  - `bc`: Plane spanned by lattice vectors **b** and **c**.
+  - `ca`: Plane spanned by lattice vectors **c** and **a**.
+- **Default**: ab
+
+### exciton_slice_pos
+
+- **Type**: Real
+- **Description**: Offset of the slice along the direction of the remaining lattice vector: **c** for `ab`, **a** for `bc`, and **b** for `ca`.
+- **Default**: 0.0
+- **Unit**: Bohr
+
+### exciton_slice_npoints
+
+- **Type**: Integer
+- **Description**: Target in-plane grid resolution. The final number of points is adjusted to use a uniform number of points per primitive cell over the plotted BvK range and includes both endpoints.
+- **Default**: 200
+
+### exciton_slice_scale
+
+- **Type**: Real
+- **Description**: Scale factor for the in-plane range relative to the BvK supercell. Values smaller than 1.0 are treated as 1.0.
+- **Default**: 1.3
+
+### ri_hartree_benchmark
+
+- **Type**: String
+- **Description**: Whether to use the localized resolution-of-identity (LRI) approximation for the **Hartree** term of kernel in the $A$ matrix of LR-TDDFT for benchmark (with FHI-aims or another ABACUS calculation). 
+  - `aims`: The `read_file_dir` directory should contain the FHI-aims output files: RI-LVL tensors `Cs_data_0.txt` and `coulomb_mat_0.txt`, and KS eigenstates from FHI-aims: `band_out`and `KS_eigenvectors.out`. The Casida equation will be constructed under FHI-aims' KS eigenpairs.
+    - LRI tensor files (`Cs_data_0.txt` and `coulomb_mat_0.txt`)and Kohn-Sham eigenvalues (`bands_out`): run FHI-aims with periodic boundary conditions and with `total_energy_method rpa` and `output librpa`.
+    - Kohn-Sham eigenstates under aims NAOs (`KS_eigenvectors.out`): run FHI-aims with `output eigenvectors`.
+    - If the number of atomic orbitals of any atom type in FHI-aims is different from that in ABACUS, the `aims_nbasis` should be set.
+  - `abacus`: The `read_file_dir` directory should contain the RI-LVL tensors `Cs` and `Vs` (written by setting `out_ri_cv` to 1). The Casida equation will be constructed under ABACUS' KS eigenpairs, with the only difference that the Hartree term is constructed with RI approximation.
+  - `aims-librpa`: The current directory where you are running task should contain the the LRI tensors `Cs_data_${processor_rank}.txt`, `coulomb_mat${processor_rank}.txt`, and KS eigenstates `band_out`, `KS_eigenvector_${processor_rank}.dat`. All these files can get from FHI-aims task with periodic boundary conditions and with `total_energy_method rpa` and `output librpa`. The `aims_nbasis` should also 0be set.
+  - `abacus-librpa`: The current directory where you are running task should contain the the LRI tensors `Cs_data_${processor_rank}.txt`, `coulomb_mat_${processor_rank}.txt`, and KS eigenstates `band_out`, `KS_eigenvector_${kpoint}.dat`. All these files can get from ABACUS task with `RPA 1`.
+  - `none`: Construct the Hartree term by Poisson equation and grid integration as usual.
+- **Default**: none
+
+### aims_nbasis
+
+- **Type**: A number(ntype) of Integers
+- **Availability**: `ri_hartree_benchmark` = `aims` or `aims-librpa`
+- **Description**: Atomic basis set size for each atom type (with the same order as in `STRU`) in FHI-aims.
+- **Default**: {} (empty list, where ABACUS use its own basis set size)
+
 [back to top](#full-list-of-input-keywords)
 
 ## Reduced Density Matrix Functional Theory
 
+Ab-initio methods and the xc-functional parameters used in RDMFT.
+The physical quantities that RDMFT temporarily expects to output are the kinetic energy, total energy, and 1-RDM of the system in the ground state, etc.
+
 ### rdmft
 
 - **Type**: Boolean

source/source_base/module_external/lapack_connector.h

@@ -318,6 +318,10 @@ void dsysv_(const char* uplo, const int* n, const int* nrhs,
             double* b, const int* ldb,
             double* work, const int* lwork,
             int* info);
+
+// === calculate matrix norm ===
+double dlange_(const char* norm, const int* m, const int* n, const double* A, const int* lda, double* work);
+double zlange_(const char* norm, const int* m, const int* n, const std::complex<double>* A, const int* lda, double* work);
 }  // extern "C"
 
 #ifdef GATHER_INFO
@@ -598,5 +602,18 @@ namespace LapackConnector
         const char trans_changed = change_trans_NC(trans);
         cherk_(&uplo_changed, &trans_changed, &n, &k, &alpha, A, &lda, &beta, C, &ldc);
     }
-} // namespace LapackConnector
-#endif  // LAPACK_CONNECTOR_HPP
+
+    static inline
+    double lange(const char& norm, const int& m, const int& n, const double* A, const int& lda, double* work)
+    {
+        return dlange_(&norm, &m, &n, A, &lda, work);
+    }
+
+    static inline
+    double lange(const char& norm, const int& m, const int& n, const std::complex<double>* A, const int& lda, double* work)
+    {
+        return zlange_(&norm, &m, &n, A, &lda, work);
+    }
+
+};
+#endif  // LAPACKCONNECTOR_HPP

source/source_base/module_external/scalapack_connector.h

@@ -42,11 +42,23 @@ extern "C"
     void pdtran_(const int* m, const int* n,
         const double* alpha, const double* a, const int* ia, const int* ja, const  int* desca,
         const double* beta, double* c, const int* ic, const int* jc, const int* descc);
-
 	void pztranu_(const int *m,const int*n,
         const std::complex<double>* alpha, const std::complex<double>* a, const int* ia, const int* ja, const int* desca,
-		const std::complex<double> *beta ,  std::complex<double> *c , const int *ic ,const int *jc ,const int *descc);
-
+		const std::complex<double> *beta , std::complex<double> *c, const int *ic, const int *jc, const int *descc);
+	void pztranc_(const int *m, const int *n,
+		const std::complex<double> *alpha, const std::complex<double> *a, const int *ia, const int *ja, const int *desca,
+		const std::complex<double> *beta, std::complex<double> *c, const int *ic, const int *jc, const int *descc);
+
+	double pdlange_(const char* norm,
+		const int* m, const int* n,
+		const double* a, const int* ia, const int* ja, const int* desca,
+		double* work);
+
+	double pzlange_(const char* norm,
+		const int* m, const int* n,
+		const std::complex<double>* a, const int* ia, const int* ja, const int* desca,
+		double* work);
+
 	void pzgemv_(
 		const char *transa,
 		const int *M, const int *N,
@@ -144,13 +156,6 @@ extern "C"
                 		const std::complex<double> *beta,
 		const std::complex<double> *c, const int *ic, const int *jc, const int *descc);
 
-    void pztranc_(
-		const int *M, const int *N,
-		const std::complex<double> *alpha,
-		const std::complex<double> *A, const int *IA, const int *JA, const int *DESCA,
-		const std::complex<double> *beta,
-		std::complex<double> *C, const int *IC, const int *JC, const int *DESCC);
-
     void pdgemr2d_(const int *M, const int *N,
 	    double *A, const int *IA, const int *JA, const int *DESCA, 
 		double *B, const int *IB, const int *JB, const int *DESCB,
@@ -394,11 +399,54 @@ class ScalapackConnector
 	static inline
 	void tranu(
 		const int m, const int n,
-		const std::complex<double> alpha , std::complex<double> *a , const int ia , const int ja , const int *desca,
-		const std::complex<double> beta ,  std::complex<double> *c , const int ic , const int jc , const int *descc)
+		const double alpha, double *a, const int ia, const int ja, const int *desca,
+		const double beta, double *c, const int ic, const int jc, const int *descc)
+	{
+		pdtran_(&m, &n, &alpha, a, &ia, &ja, desca, &beta, c, &ic, &jc, descc);
+	}
+	static inline
+	void tranu(
+		const int m, const int n,
+		const std::complex<double> alpha, std::complex<double> *a, const int ia, const int ja, const int *desca,
+		const std::complex<double> beta, std::complex<double> *c, const int ic, const int jc, const int *descc)
 	{
 		pztranu_(&m, &n, &alpha, a, &ia, &ja, desca, &beta, c, &ic, &jc, descc);
 	}
+	static inline
+	void tranc(
+		const int m, const int n,
+		const double alpha, double *a, const int ia, const int ja, const int *desca,
+		const double beta, double *c, const int ic, const int jc, const int *descc)
+	{
+		pdtran_(&m, &n, &alpha, a, &ia, &ja, desca, &beta, c, &ic, &jc, descc);
+	}
+	static inline
+	void tranc(
+		const int m, const int n,
+		const std::complex<double> alpha, std::complex<double> *a, const int ia, const int ja, const int *desca,
+		const std::complex<double> beta, std::complex<double> *c, const int ic, const int jc, const int *descc)
+	{
+		pztranc_(&m, &n, &alpha, a, &ia, &ja, desca, &beta, c, &ic, &jc, descc);
+	}
+
+    static inline
+    double lange(
+        const char norm,
+        const int m, const int n,
+        const double* a, const int ia, const int ja, const int* desca,
+        double* work)
+    {
+        return pdlange_(&norm, &m, &n, a, &ia, &ja, desca, work);
+    }
+    static inline
+    double lange(
+        const char norm,
+        const int m, const int n,
+        const std::complex<double>* a, const int ia, const int ja, const int* desca,
+        double* work)
+    {
+        return pzlange_(&norm, &m, &n, a, &ia, &ja, desca, work);
+    }
 
 	static inline
 	int potrf(char uplo, int na, double* U, int* desc)