Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Main routines for GW + Bethe-Salpeter for computing electronic excitations
10 : !> \par History
11 : !> 04.2024 created [Maximilian Graml]
12 : ! **************************************************************************************************
13 :
14 : MODULE bse_main
15 :
16 : USE bse_full_diag, ONLY: create_A,&
17 : create_B,&
18 : create_hermitian_form_of_ABBA,&
19 : diagonalize_A,&
20 : diagonalize_C
21 : USE bse_iterative, ONLY: do_subspace_iterations,&
22 : fill_local_3c_arrays
23 : USE bse_print, ONLY: print_BSE_start_flag
24 : USE bse_util, ONLY: adapt_BSE_input_params,&
25 : deallocate_matrices_bse,&
26 : estimate_BSE_resources,&
27 : mult_B_with_W,&
28 : truncate_BSE_matrices
29 : USE cp_fm_types, ONLY: cp_fm_release,&
30 : cp_fm_type
31 : USE cp_log_handling, ONLY: cp_get_default_logger,&
32 : cp_logger_type
33 : USE cp_output_handling, ONLY: debug_print_level
34 : USE group_dist_types, ONLY: group_dist_d1_type
35 : USE input_constants, ONLY: bse_abba,&
36 : bse_both,&
37 : bse_fulldiag,&
38 : bse_iterdiag,&
39 : bse_tda
40 : USE kinds, ONLY: dp
41 : USE message_passing, ONLY: mp_para_env_type
42 : USE mp2_types, ONLY: mp2_type
43 : USE qs_environment_types, ONLY: qs_environment_type
44 : #include "./base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 :
48 : PRIVATE
49 :
50 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_main'
51 :
52 : PUBLIC :: start_bse_calculation
53 :
54 : CONTAINS
55 :
56 : ! **************************************************************************************************
57 : !> \brief Main subroutine managing BSE calculations
58 : !> \param fm_mat_S_ia_bse ...
59 : !> \param fm_mat_S_ij_bse ...
60 : !> \param fm_mat_S_ab_bse ...
61 : !> \param fm_mat_Q_static_bse_gemm ...
62 : !> \param Eigenval ...
63 : !> \param Eigenval_scf ...
64 : !> \param homo ...
65 : !> \param virtual ...
66 : !> \param dimen_RI ...
67 : !> \param dimen_RI_red ...
68 : !> \param bse_lev_virt ...
69 : !> \param gd_array ...
70 : !> \param color_sub ...
71 : !> \param mp2_env ...
72 : !> \param qs_env ...
73 : !> \param mo_coeff ...
74 : !> \param unit_nr ...
75 : ! **************************************************************************************************
76 32 : SUBROUTINE start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
77 : fm_mat_Q_static_bse_gemm, &
78 : Eigenval, Eigenval_scf, &
79 32 : homo, virtual, dimen_RI, dimen_RI_red, bse_lev_virt, &
80 32 : gd_array, color_sub, mp2_env, qs_env, mo_coeff, unit_nr)
81 :
82 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_ij_bse, &
83 : fm_mat_S_ab_bse
84 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_Q_static_bse_gemm
85 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
86 : INTENT(IN) :: Eigenval, Eigenval_scf
87 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
88 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, bse_lev_virt
89 : TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
90 : INTEGER, INTENT(IN) :: color_sub
91 : TYPE(mp2_type) :: mp2_env
92 : TYPE(qs_environment_type), POINTER :: qs_env
93 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
94 : INTEGER, INTENT(IN) :: unit_nr
95 :
96 : CHARACTER(LEN=*), PARAMETER :: routineN = 'start_bse_calculation'
97 :
98 : INTEGER :: handle, homo_red, virtual_red
99 : LOGICAL :: my_do_abba, my_do_fulldiag, &
100 : my_do_iterat_diag, my_do_tda
101 : REAL(KIND=dp) :: diag_runtime_est
102 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Eigenval_reduced
103 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_abQ_bse_local, B_bar_iaQ_bse_local, &
104 32 : B_bar_ijQ_bse_local, B_iaQ_bse_local
105 : TYPE(cp_fm_type) :: fm_A_BSE, fm_B_BSE, fm_C_BSE, fm_inv_sqrt_A_minus_B, fm_mat_S_ab_trunc, &
106 : fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, &
107 : fm_sqrt_A_minus_B
108 : TYPE(cp_logger_type), POINTER :: logger
109 : TYPE(mp_para_env_type), POINTER :: para_env
110 :
111 32 : CALL timeset(routineN, handle)
112 :
113 32 : para_env => fm_mat_S_ia_bse%matrix_struct%para_env
114 :
115 32 : my_do_fulldiag = .FALSE.
116 32 : my_do_iterat_diag = .FALSE.
117 32 : my_do_tda = .FALSE.
118 32 : my_do_abba = .FALSE.
119 : !Method: Iterative or full diagonalization
120 32 : SELECT CASE (mp2_env%bse%bse_diag_method)
121 : CASE (bse_iterdiag)
122 0 : my_do_iterat_diag = .TRUE.
123 : !MG: Basics of the Davidson solver are implemented, but not rigorously checked.
124 0 : CPABORT("Iterative BSE not yet implemented")
125 : CASE (bse_fulldiag)
126 32 : my_do_fulldiag = .TRUE.
127 : END SELECT
128 : !Approximation: TDA and/or full ABBA matrix
129 48 : SELECT CASE (mp2_env%bse%flag_tda)
130 : CASE (bse_tda)
131 16 : my_do_tda = .TRUE.
132 : CASE (bse_abba)
133 16 : my_do_abba = .TRUE.
134 : CASE (bse_both)
135 0 : my_do_tda = .TRUE.
136 32 : my_do_abba = .TRUE.
137 : END SELECT
138 :
139 32 : CALL print_BSE_start_flag(my_do_tda, my_do_abba, unit_nr)
140 :
141 : ! Link BSE debug flag against debug print level
142 32 : logger => cp_get_default_logger()
143 32 : IF (logger%iter_info%print_level == debug_print_level) THEN
144 0 : mp2_env%bse%bse_debug_print = .TRUE.
145 : END IF
146 :
147 32 : CALL fm_mat_S_ia_bse%matrix_struct%para_env%sync()
148 : ! We apply the BSE cutoffs using the DFT Eigenenergies
149 : ! Reduce matrices in case of energy cutoff for occupied and unoccupied in A/B-BSE-matrices
150 : CALL truncate_BSE_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
151 : fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
152 : Eigenval_scf(:, 1, 1), Eigenval(:, 1, 1), Eigenval_reduced, &
153 : homo(1), virtual(1), dimen_RI, unit_nr, &
154 : bse_lev_virt, &
155 : homo_red, virtual_red, &
156 32 : mp2_env)
157 : ! \bar{B}^P_rs = \sum_R W_PR B^R_rs where B^R_rs = \sum_T [1/sqrt(v)]_RT (T|rs)
158 : ! r,s: MO-index, P,R,T: RI-index
159 : ! B: fm_mat_S_..., W: fm_mat_Q_...
160 : CALL mult_B_with_W(fm_mat_S_ij_trunc, fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, &
161 : fm_mat_S_bar_ij_bse, fm_mat_Q_static_bse_gemm, &
162 32 : dimen_RI_red, homo_red, virtual_red)
163 :
164 32 : IF (my_do_iterat_diag) THEN
165 : CALL fill_local_3c_arrays(fm_mat_S_ab_trunc, fm_mat_S_ia_trunc, &
166 : fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
167 : B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
168 : B_iaQ_bse_local, dimen_RI_red, homo_red, virtual_red, &
169 : gd_array, color_sub, para_env)
170 : END IF
171 :
172 32 : CALL adapt_BSE_input_params(homo_red, virtual_red, unit_nr, mp2_env)
173 :
174 32 : IF (my_do_fulldiag) THEN
175 : ! Quick estimate of memory consumption and runtime of diagonalizations
176 : CALL estimate_BSE_resources(homo_red, virtual_red, unit_nr, my_do_abba, &
177 32 : para_env, diag_runtime_est)
178 : ! Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W)
179 : ! A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
180 : ! ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
181 : ! α is a spin-dependent factor
182 : ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
183 : ! W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
184 : CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_bar_ij_bse, fm_mat_S_ab_trunc, &
185 : fm_A_BSE, Eigenval_reduced, unit_nr, &
186 : homo_red, virtual_red, dimen_RI, mp2_env, &
187 32 : para_env)
188 32 : IF (my_do_abba) THEN
189 : ! Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W)
190 : ! B_ia,jb = α * v_ia,jb - W_ib,aj
191 : ! α is a spin-dependent factor
192 : ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
193 : ! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
194 : CALL create_B(fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, fm_B_BSE, &
195 16 : homo_red, virtual_red, dimen_RI, unit_nr, mp2_env)
196 : ! Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
197 : ! (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
198 : ! We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
199 : ! of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
200 : CALL create_hermitian_form_of_ABBA(fm_A_BSE, fm_B_BSE, fm_C_BSE, &
201 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
202 16 : homo_red, virtual_red, unit_nr, mp2_env, diag_runtime_est)
203 : END IF
204 32 : CALL cp_fm_release(fm_B_BSE)
205 32 : IF (my_do_tda) THEN
206 : ! Solving the hermitian eigenvalue equation A X^n = Ω^n X^n
207 : CALL diagonalize_A(fm_A_BSE, homo_red, virtual_red, homo(1), &
208 16 : unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
209 : END IF
210 : ! Release to avoid faulty use of changed A matrix
211 32 : CALL cp_fm_release(fm_A_BSE)
212 32 : IF (my_do_abba) THEN
213 : ! Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
214 : ! Here, the eigenvectors Z^n relate to X^n via
215 : ! Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
216 : CALL diagonalize_C(fm_C_BSE, homo_red, virtual_red, homo(1), &
217 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
218 16 : unit_nr, diag_runtime_est, mp2_env, qs_env, mo_coeff)
219 : END IF
220 : ! Release to avoid faulty use of changed C matrix
221 32 : CALL cp_fm_release(fm_C_BSE)
222 : END IF
223 :
224 : CALL deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
225 : fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
226 32 : fm_mat_Q_static_bse_gemm, mp2_env)
227 32 : DEALLOCATE (Eigenval_reduced)
228 32 : IF (my_do_iterat_diag) THEN
229 : ! Contains untested Block-Davidson algorithm
230 : CALL do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
231 : B_iaQ_bse_local, homo(1), virtual(1), mp2_env%bse%bse_spin_config, unit_nr, &
232 0 : Eigenval(:, 1, 1), para_env, mp2_env)
233 : ! Deallocate local 3c-B-matrices
234 0 : DEALLOCATE (B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, B_iaQ_bse_local)
235 : END IF
236 :
237 32 : IF (unit_nr > 0) THEN
238 16 : WRITE (unit_nr, '(T2,A4,T7,A53)') 'BSE|', 'The BSE was successfully calculated. Have a nice day!'
239 : END IF
240 :
241 32 : CALL timestop(handle)
242 :
243 64 : END SUBROUTINE start_bse_calculation
244 :
245 : END MODULE bse_main
|
