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 Iterative routines for GW + Bethe-Salpeter for computing electronic excitations
10 : !> \par History
11 : !> 04.2017 created [Jan Wilhelm]
12 : !> 11.2023 Davidson solver implemented [Maximilian Graml]
13 : ! **************************************************************************************************
14 : MODULE bse_iterative
15 : USE cp_fm_types, ONLY: cp_fm_get_info,&
16 : cp_fm_type
17 : USE group_dist_types, ONLY: get_group_dist,&
18 : group_dist_d1_type
19 : USE input_constants, ONLY: bse_singlet,&
20 : bse_triplet
21 : USE kinds, ONLY: dp
22 : USE message_passing, ONLY: mp_para_env_type,&
23 : mp_request_type
24 : USE mp2_types, ONLY: integ_mat_buffer_type,&
25 : mp2_type
26 : USE physcon, ONLY: evolt
27 : USE rpa_communication, ONLY: communicate_buffer
28 : #include "./base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 :
32 : PRIVATE
33 :
34 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_iterative'
35 :
36 : PUBLIC :: fill_local_3c_arrays, do_subspace_iterations
37 :
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief ...
42 : !> \param B_bar_ijQ_bse_local ...
43 : !> \param B_abQ_bse_local ...
44 : !> \param B_bar_iaQ_bse_local ...
45 : !> \param B_iaQ_bse_local ...
46 : !> \param homo ...
47 : !> \param virtual ...
48 : !> \param bse_spin_config ...
49 : !> \param unit_nr ...
50 : !> \param Eigenval ...
51 : !> \param para_env ...
52 : !> \param mp2_env ...
53 : ! **************************************************************************************************
54 0 : SUBROUTINE do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
55 : B_iaQ_bse_local, homo, virtual, bse_spin_config, unit_nr, &
56 0 : Eigenval, para_env, mp2_env)
57 :
58 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_bar_ijQ_bse_local, B_abQ_bse_local, &
59 : B_bar_iaQ_bse_local, B_iaQ_bse_local
60 : INTEGER :: homo, virtual, bse_spin_config, unit_nr
61 : REAL(KIND=dp), DIMENSION(:) :: Eigenval
62 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
63 : TYPE(mp2_type) :: mp2_env
64 :
65 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_subspace_iterations'
66 :
67 : CHARACTER(LEN=10) :: bse_davidson_abort_cond_string, &
68 : success_abort_string
69 : INTEGER :: bse_davidson_abort_cond, davidson_converged, fac_max_z_space, handle, i_iter, &
70 : j_print, local_RI_size, num_add_start_z_space, num_davidson_iter, num_en_unconverged, &
71 : num_exact_en_unconverged, num_exc_en, num_max_z_space, num_new_t, num_res_unconverged, &
72 : num_Z_vectors, num_Z_vectors_init
73 : LOGICAL :: bse_full_diag_debug
74 : REAL(kind=dp) :: eps_exc_en, eps_res, max_en_diff, &
75 : max_res_norm, z_space_energy_cutoff
76 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: En_diffs, En_diffs_exact, Full_exc_spectrum, &
77 0 : Res_norms, Subspace_full_eigenval, Subspace_new_eigenval, Subspace_prev_eigenval
78 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: AZ_reshaped, M_ia_tmp, M_ji_tmp, RI_vector, &
79 0 : Subspace_new_eigenvec, Subspace_residuals_reshaped, Z_vectors_reshaped
80 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: AZ, BZ, Subspace_add_dir, &
81 0 : Subspace_ritzvec, W_vectors, Z_vectors
82 :
83 0 : CALL timeset(routineN, handle)
84 :
85 : !MG to del
86 : !Debug flag for exact diagonalization (only using lapack!!!)
87 0 : bse_full_diag_debug = .TRUE.
88 0 : num_en_unconverged = -1
89 0 : num_res_unconverged = -1
90 0 : num_exact_en_unconverged = -1
91 :
92 0 : bse_davidson_abort_cond = mp2_env%ri_g0w0%davidson_abort_cond
93 0 : num_exc_en = mp2_env%ri_g0w0%num_exc_en
94 0 : num_add_start_z_space = mp2_env%ri_g0w0%num_add_start_z_space
95 0 : fac_max_z_space = mp2_env%ri_g0w0%fac_max_z_space
96 0 : num_new_t = mp2_env%ri_g0w0%num_new_t
97 0 : num_davidson_iter = mp2_env%ri_g0w0%num_davidson_iter
98 0 : eps_res = mp2_env%ri_g0w0%eps_res
99 0 : eps_exc_en = mp2_env%ri_g0w0%eps_exc_en
100 0 : z_space_energy_cutoff = mp2_env%ri_g0w0%z_space_energy_cutoff
101 :
102 0 : num_Z_vectors_init = num_exc_en + num_add_start_z_space
103 :
104 0 : IF (unit_nr > 0) THEN
105 0 : WRITE (unit_nr, *) "bse_spin_config", bse_spin_config
106 0 : WRITE (unit_nr, *) "num_exc_en", num_exc_en
107 0 : WRITE (unit_nr, *) "num_add_start_z_space", num_add_start_z_space
108 0 : WRITE (unit_nr, *) "num_Z_vectors_init", num_Z_vectors_init
109 0 : WRITE (unit_nr, *) "fac_max_z_space", fac_max_z_space
110 0 : WRITE (unit_nr, *) "num_new_t", num_new_t
111 0 : WRITE (unit_nr, *) "eps_res", eps_res
112 0 : WRITE (unit_nr, *) "num_davidson_iter", num_davidson_iter
113 0 : WRITE (unit_nr, *) "eps_exc_en", eps_exc_en
114 0 : WRITE (unit_nr, *) "bse_davidson_abort_cond", bse_davidson_abort_cond
115 0 : WRITE (unit_nr, *) "z_space_energy_cutoff", z_space_energy_cutoff
116 0 : WRITE (unit_nr, *) "Printing B_bar_iaQ_bse_local of shape", SHAPE(B_bar_iaQ_bse_local)
117 : END IF
118 :
119 0 : local_RI_size = SIZE(B_iaQ_bse_local, 3)
120 :
121 0 : num_Z_vectors = num_Z_vectors_init
122 0 : num_max_z_space = num_Z_vectors_init*fac_max_z_space
123 :
124 : !Check input parameters and correct them if necessary
125 0 : IF (num_new_t > num_Z_vectors_init) THEN
126 0 : num_new_t = num_Z_vectors_init
127 0 : IF (unit_nr > 0) THEN
128 : CALL cp_warn(__LOCATION__, "Number of added directions has to be smaller/equals than "// &
129 0 : "initial dimension. Corrected num_new_t accordingly.")
130 : END IF
131 : END IF
132 0 : IF (unit_nr > 0) THEN
133 0 : WRITE (unit_nr, *) "Between BSE correction Warnings"
134 : END IF
135 : !If initial number is too big, already the first iteration causes trouble in LAPACK diagonal. (DORGQR)
136 0 : IF (2*num_Z_vectors_init > homo*virtual) THEN
137 : CALL cp_abort(__LOCATION__, "Initial dimension was too large and could not be corrected. "// &
138 0 : "Choose another num_exc_en and num_add_start_z_space or adapt your basis set.")
139 : END IF
140 0 : IF (num_max_z_space .GE. homo*virtual) THEN
141 0 : fac_max_z_space = homo*virtual/num_Z_vectors_init
142 0 : num_max_z_space = num_Z_vectors_init*fac_max_z_space
143 :
144 0 : IF (fac_max_z_space == 0) THEN
145 : CALL cp_abort(__LOCATION__, "Maximal dimension was too large and could not be corrected. "// &
146 0 : "Choose another fac_max_z_space and num_Z_vectors_init or adapt your basis set.")
147 : ELSE
148 0 : IF (unit_nr > 0) THEN
149 : CALL cp_warn(__LOCATION__, "Maximal dimension of Z space has to be smaller than homo*virtual. "// &
150 0 : "Corrected fac_max_z_space accordingly.")
151 : END IF
152 : END IF
153 : END IF
154 :
155 0 : DO i_iter = 1, num_davidson_iter
156 0 : IF (unit_nr > 0) THEN
157 0 : WRITE (unit_nr, *) "Allocating Z_vec,AZ,BZ with dimensions (homo,virt,num_Z)", homo, virtual, num_Z_vectors
158 0 : WRITE (unit_nr, *) 'ProcNr', para_env%mepos, 'you really enter here for i_iter', i_iter
159 : END IF
160 0 : ALLOCATE (Z_vectors(homo, virtual, num_Z_vectors))
161 0 : Z_vectors = 0.0_dp
162 :
163 : !Dellocation procedures are a bit intricate, W_/Z_vectors and eigenvalues are needed for the next iteration,
164 : ! therefore we have to deallocate them separately from the other quantities
165 0 : IF (i_iter == 1) THEN
166 0 : CALL initial_guess_Z_vectors(Z_vectors, Eigenval, num_Z_vectors, homo, virtual)
167 0 : ALLOCATE (Subspace_prev_eigenval(num_exc_en))
168 0 : Subspace_prev_eigenval = 0.0_dp
169 : ELSE
170 0 : Z_vectors(:, :, :) = W_vectors(:, :, :)
171 0 : DEALLOCATE (W_vectors)
172 : END IF
173 0 : IF (unit_nr > 0) THEN
174 0 : WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "Allocated/rewritten Z arrays"
175 : END IF
176 :
177 : CALL create_bse_work_arrays(AZ, Z_vectors_reshaped, AZ_reshaped, BZ, M_ia_tmp, M_ji_tmp, &
178 : RI_vector, Subspace_new_eigenval, Subspace_full_eigenval, Subspace_new_eigenvec, &
179 : Subspace_residuals_reshaped, Subspace_ritzvec, Subspace_add_dir, W_vectors, &
180 0 : homo, virtual, num_Z_vectors, local_RI_size, num_new_t)
181 0 : IF (unit_nr > 0) THEN
182 0 : WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "Allocated Work arrays"
183 : END IF
184 :
185 : CALL compute_AZ(AZ, Z_vectors, B_iaQ_bse_local, B_bar_ijQ_bse_local, B_abQ_bse_local, &
186 : M_ia_tmp, RI_vector, Eigenval, homo, virtual, num_Z_vectors, local_RI_size, &
187 : para_env, bse_spin_config, z_space_energy_cutoff, i_iter, bse_full_diag_debug, &
188 0 : Full_exc_spectrum, unit_nr)
189 :
190 : !MG: functionality of BZ not checked (issue with fm_mat_Q_static_bse_gemm in rpa_util needs to be checked!)
191 : !CALL compute_BZ(BZ, Z_vectors, B_iaQ_bse_local, B_bar_iaQ_bse_local, &
192 : ! M_ji_tmp, homo, virtual, num_Z_vectors, local_RI_size, &
193 : ! para_env)
194 :
195 0 : IF (unit_nr > 0) THEN
196 0 : WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "Computed AZ"
197 : END IF
198 :
199 : !MG to check: Reshaping correct?
200 0 : AZ_reshaped(:, :) = RESHAPE(AZ, (/homo*virtual, num_Z_vectors/))
201 0 : Z_vectors_reshaped(:, :) = RESHAPE(Z_vectors, (/homo*virtual, num_Z_vectors/))
202 :
203 : ! Diagonalize M and extract smallest eigenvalues/corresponding eigenvector
204 : CALL compute_diagonalize_ZAZ(AZ_reshaped, Z_vectors_reshaped, num_Z_vectors, Subspace_new_eigenval, &
205 0 : Subspace_new_eigenvec, num_new_t, Subspace_full_eigenval, para_env, unit_nr)
206 0 : IF (unit_nr > 0) THEN
207 0 : WRITE (unit_nr, *) "Eigenval (eV) in iter=", i_iter, " is:", Subspace_new_eigenval(:6)*evolt
208 : END IF
209 :
210 : ! Threshold in energies
211 : CALL check_en_convergence(Subspace_full_eigenval, Subspace_prev_eigenval, eps_exc_en, num_en_unconverged, &
212 0 : num_exc_en, max_en_diff, En_diffs)
213 0 : IF (unit_nr > 0) THEN
214 0 : WRITE (unit_nr, *) "Largest change of desired exc ens =", max_en_diff
215 : END IF
216 : ! Compute residuals
217 : CALL compute_residuals(AZ_reshaped, Z_vectors_reshaped, Subspace_new_eigenval, Subspace_new_eigenvec, &
218 0 : Subspace_residuals_reshaped, homo, virtual, num_new_t, num_Z_vectors, Subspace_ritzvec)
219 :
220 : !Abort, if residuals are small enough w.r.t threshold
221 : CALL check_res_convergence(Subspace_residuals_reshaped, num_new_t, eps_res, num_res_unconverged, &
222 0 : i_iter, max_res_norm, unit_nr, Res_norms)
223 :
224 0 : davidson_converged = -1
225 0 : IF (num_res_unconverged == 0 .AND. bse_davidson_abort_cond /= 0) THEN
226 0 : davidson_converged = 1
227 0 : success_abort_string = "RESIDUALS"
228 0 : ELSE IF (num_en_unconverged == 0 .AND. (bse_davidson_abort_cond /= 1)) THEN
229 0 : davidson_converged = 1
230 0 : success_abort_string = "ENERGIES"
231 0 : ELSE IF (i_iter == num_davidson_iter) THEN
232 0 : davidson_converged = -100
233 0 : success_abort_string = "-----"
234 : ELSE
235 : davidson_converged = -1
236 : END IF
237 :
238 0 : IF (bse_davidson_abort_cond == 0) THEN
239 0 : bse_davidson_abort_cond_string = "ENERGY"
240 0 : ELSE IF (bse_davidson_abort_cond == 1) THEN
241 0 : bse_davidson_abort_cond_string = "RESIDUAL"
242 : ELSE
243 0 : bse_davidson_abort_cond_string = "EITHER"
244 : END IF
245 :
246 0 : IF (davidson_converged == 1) THEN
247 : CALL success_message(Subspace_full_eigenval, num_new_t, eps_res, num_res_unconverged, &
248 : bse_spin_config, unit_nr, num_exc_en, num_Z_vectors_init, &
249 : num_davidson_iter, i_iter, num_Z_vectors, num_max_z_space, max_res_norm, &
250 : max_en_diff, num_en_unconverged, bse_davidson_abort_cond_string, &
251 0 : eps_exc_en, success_abort_string, z_space_energy_cutoff)
252 :
253 : !Deallocate matrices, which are otherwise not cleared due to exiting the loop
254 0 : DEALLOCATE (AZ, BZ, &
255 0 : Z_vectors, M_ia_tmp, M_ji_tmp, RI_vector, Subspace_prev_eigenval, &
256 0 : Subspace_new_eigenval, Subspace_new_eigenvec, Subspace_residuals_reshaped, &
257 0 : Subspace_add_dir, AZ_reshaped, Z_vectors_reshaped, Subspace_ritzvec, Subspace_full_eigenval)
258 :
259 0 : EXIT
260 0 : ELSE IF (davidson_converged < -1) THEN
261 : CALL print_davidson_parameter(i_iter, num_davidson_iter, num_Z_vectors, num_res_unconverged, max_res_norm, &
262 : eps_res, num_en_unconverged, max_en_diff, eps_exc_en, num_exc_en, &
263 : num_Z_vectors_init, num_max_z_space, num_new_t, unit_nr, &
264 0 : success_abort_string, bse_davidson_abort_cond_string, z_space_energy_cutoff)
265 :
266 : CALL cp_abort(__LOCATION__, "BSE/TDA-Davidson did not converge using "// &
267 0 : bse_davidson_abort_cond_string//" threshold condition!")
268 : END IF
269 :
270 : ! Calculate and add next orthonormal vector and update num_Z_vectors
271 : CALL compute_new_directions(homo, virtual, Subspace_residuals_reshaped, Subspace_new_eigenval, Eigenval, &
272 0 : num_new_t, Subspace_add_dir)
273 :
274 : !If exact-diag: compute difference to exact eigenvalues
275 : IF (bse_full_diag_debug) THEN
276 0 : ALLOCATE (En_diffs_exact(num_exc_en))
277 0 : num_exact_en_unconverged = 0
278 0 : DO j_print = 1, num_exc_en
279 0 : En_diffs_exact(j_print) = ABS(Subspace_full_eigenval(j_print) - Full_exc_spectrum(j_print))
280 0 : IF (En_diffs_exact(j_print) > eps_exc_en) num_exact_en_unconverged = num_exact_en_unconverged + 1
281 : END DO
282 : END IF
283 :
284 : !Check dimensions and orthonormalize vector system, depending on dimensionality
285 : CALL check_Z_space_dimension(W_vectors, Z_vectors, Subspace_add_dir, Subspace_ritzvec, &
286 0 : num_Z_vectors, num_new_t, num_max_z_space, homo, virtual, i_iter, unit_nr)
287 :
288 : !Copy eigenvalues for threshold
289 0 : Subspace_prev_eigenval(:) = Subspace_full_eigenval(:num_exc_en)
290 :
291 0 : DEALLOCATE (AZ, & !BZ,
292 0 : Z_vectors, M_ia_tmp, M_ji_tmp, RI_vector, &
293 0 : Subspace_new_eigenval, Subspace_new_eigenvec, Subspace_residuals_reshaped, &
294 0 : Subspace_add_dir, AZ_reshaped, Z_vectors_reshaped, Subspace_ritzvec, Subspace_full_eigenval, &
295 0 : Res_norms, En_diffs)
296 :
297 : IF (bse_full_diag_debug) THEN
298 0 : DEALLOCATE (En_diffs_exact)
299 : END IF
300 :
301 : !Orthonorm:
302 0 : CALL orthonormalize_W(W_vectors, num_Z_vectors, homo, virtual)
303 :
304 : END DO
305 :
306 0 : CALL timestop(handle)
307 :
308 0 : END SUBROUTINE
309 :
310 : ! **************************************************************************************************
311 : !> \brief ...
312 : !> \param W_vectors ...
313 : !> \param Z_vectors ...
314 : !> \param Subspace_add_dir ...
315 : !> \param Subspace_ritzvec ...
316 : !> \param num_Z_vectors ...
317 : !> \param num_new_t ...
318 : !> \param num_max_z_space ...
319 : !> \param homo ...
320 : !> \param virtual ...
321 : !> \param i_iter ...
322 : !> \param unit_nr ...
323 : ! **************************************************************************************************
324 0 : SUBROUTINE check_Z_space_dimension(W_vectors, Z_vectors, Subspace_add_dir, Subspace_ritzvec, &
325 : num_Z_vectors, num_new_t, num_max_z_space, homo, virtual, i_iter, unit_nr)
326 :
327 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: W_vectors, Z_vectors, Subspace_add_dir, &
328 : Subspace_ritzvec
329 : INTEGER :: num_Z_vectors, num_new_t, &
330 : num_max_z_space, homo, virtual, &
331 : i_iter, unit_nr
332 :
333 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_Z_space_dimension'
334 :
335 : INTEGER :: handle
336 :
337 0 : CALL timeset(routineN, handle)
338 :
339 0 : IF (num_Z_vectors + num_new_t .LE. num_max_z_space) THEN
340 0 : W_vectors(:, :, :num_Z_vectors) = Z_vectors(:, :, :)
341 0 : W_vectors(:, :, num_Z_vectors + 1:) = Subspace_add_dir
342 0 : num_Z_vectors = num_Z_vectors + num_new_t
343 : ELSE
344 0 : IF (unit_nr > 0) THEN
345 0 : WRITE (unit_nr, *) "Resetting dimension in i_iter=", i_iter
346 : END IF
347 0 : DEALLOCATE (W_vectors)
348 0 : ALLOCATE (W_vectors(homo, virtual, 2*num_new_t))
349 0 : W_vectors(:, :, :num_new_t) = Subspace_ritzvec(:, :, :)
350 0 : W_vectors(:, :, num_new_t + 1:) = Subspace_add_dir
351 0 : num_Z_vectors = 2*num_new_t
352 : END IF
353 :
354 0 : CALL timestop(handle)
355 :
356 0 : END SUBROUTINE
357 :
358 : ! **************************************************************************************************
359 : !> \brief ...
360 : !> \param AZ ...
361 : !> \param Z_vectors_reshaped ...
362 : !> \param AZ_reshaped ...
363 : !> \param BZ ...
364 : !> \param M_ia_tmp ...
365 : !> \param M_ji_tmp ...
366 : !> \param RI_vector ...
367 : !> \param Subspace_new_eigenval ...
368 : !> \param Subspace_full_eigenval ...
369 : !> \param Subspace_new_eigenvec ...
370 : !> \param Subspace_residuals_reshaped ...
371 : !> \param Subspace_ritzvec ...
372 : !> \param Subspace_add_dir ...
373 : !> \param W_vectors ...
374 : !> \param homo ...
375 : !> \param virtual ...
376 : !> \param num_Z_vectors ...
377 : !> \param local_RI_size ...
378 : !> \param num_new_t ...
379 : ! **************************************************************************************************
380 0 : SUBROUTINE create_bse_work_arrays(AZ, Z_vectors_reshaped, AZ_reshaped, BZ, M_ia_tmp, M_ji_tmp, &
381 : RI_vector, Subspace_new_eigenval, Subspace_full_eigenval, Subspace_new_eigenvec, &
382 : Subspace_residuals_reshaped, Subspace_ritzvec, Subspace_add_dir, W_vectors, &
383 : homo, virtual, num_Z_vectors, local_RI_size, num_new_t)
384 :
385 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: AZ
386 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Z_vectors_reshaped, AZ_reshaped
387 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: BZ
388 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: M_ia_tmp, M_ji_tmp, RI_vector
389 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_new_eigenval, &
390 : Subspace_full_eigenval
391 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_new_eigenvec, &
392 : Subspace_residuals_reshaped
393 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Subspace_ritzvec, Subspace_add_dir, &
394 : W_vectors
395 : INTEGER :: homo, virtual, num_Z_vectors, &
396 : local_RI_size, num_new_t
397 :
398 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_bse_work_arrays'
399 :
400 : INTEGER :: handle
401 :
402 0 : CALL timeset(routineN, handle)
403 :
404 0 : ALLOCATE (AZ(homo, virtual, num_Z_vectors))
405 0 : AZ = 0.0_dp
406 :
407 0 : ALLOCATE (Z_vectors_reshaped(homo*virtual, num_Z_vectors))
408 0 : Z_vectors_reshaped = 0.0_dp
409 :
410 0 : ALLOCATE (AZ_reshaped(homo*virtual, num_Z_vectors))
411 0 : AZ_reshaped = 0.0_dp
412 :
413 0 : ALLOCATE (BZ(homo, virtual, num_Z_vectors))
414 0 : BZ = 0.0_dp
415 :
416 0 : ALLOCATE (M_ia_tmp(homo, virtual))
417 0 : M_ia_tmp = 0.0_dp
418 :
419 0 : ALLOCATE (M_ji_tmp(homo, homo))
420 0 : M_ji_tmp = 0.0_dp
421 :
422 0 : ALLOCATE (RI_vector(local_RI_size, num_Z_vectors))
423 0 : RI_vector = 0.0_dp
424 :
425 0 : ALLOCATE (Subspace_new_eigenval(num_new_t))
426 0 : Subspace_new_eigenval = 0.0_dp
427 :
428 0 : ALLOCATE (Subspace_full_eigenval(num_Z_vectors))
429 0 : Subspace_full_eigenval = 0.0_dp
430 :
431 0 : ALLOCATE (Subspace_new_eigenvec(num_Z_vectors, num_new_t))
432 0 : Subspace_new_eigenvec = 0.0_dp
433 :
434 0 : ALLOCATE (subspace_residuals_reshaped(homo*virtual, num_new_t))
435 0 : subspace_residuals_reshaped = 0.0_dp
436 :
437 0 : ALLOCATE (Subspace_ritzvec(homo, virtual, num_new_t))
438 0 : Subspace_ritzvec = 0.0_dp
439 :
440 0 : ALLOCATE (Subspace_add_dir(homo, virtual, num_new_t))
441 0 : Subspace_add_dir = 0.0_dp
442 :
443 0 : ALLOCATE (W_vectors(homo, virtual, num_Z_vectors + num_new_t))
444 0 : W_vectors = 0.0_dp
445 :
446 0 : CALL timestop(handle)
447 :
448 0 : END SUBROUTINE
449 :
450 : ! **************************************************************************************************
451 : !> \brief ...
452 : !> \param Subspace_full_eigenval ...
453 : !> \param num_new_t ...
454 : !> \param eps_res ...
455 : !> \param num_res_unconverged ...
456 : !> \param bse_spin_config ...
457 : !> \param unit_nr ...
458 : !> \param num_exc_en ...
459 : !> \param num_Z_vectors_init ...
460 : !> \param num_davidson_iter ...
461 : !> \param i_iter ...
462 : !> \param num_Z_vectors ...
463 : !> \param num_max_z_space ...
464 : !> \param max_res_norm ...
465 : !> \param max_en_diff ...
466 : !> \param num_en_unconverged ...
467 : !> \param bse_davidson_abort_cond_string ...
468 : !> \param eps_exc_en ...
469 : !> \param success_abort_string ...
470 : !> \param z_space_energy_cutoff ...
471 : ! **************************************************************************************************
472 0 : SUBROUTINE success_message(Subspace_full_eigenval, num_new_t, eps_res, num_res_unconverged, &
473 : bse_spin_config, unit_nr, num_exc_en, num_Z_vectors_init, &
474 : num_davidson_iter, i_iter, num_Z_vectors, num_max_z_space, max_res_norm, &
475 : max_en_diff, num_en_unconverged, bse_davidson_abort_cond_string, &
476 : eps_exc_en, success_abort_string, z_space_energy_cutoff)
477 :
478 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_full_eigenval
479 : INTEGER :: num_new_t
480 : REAL(KIND=dp) :: eps_res
481 : INTEGER :: num_res_unconverged, bse_spin_config, unit_nr, num_exc_en, num_Z_vectors_init, &
482 : num_davidson_iter, i_iter, num_Z_vectors, num_max_z_space
483 : REAL(KIND=dp) :: max_res_norm, max_en_diff
484 : INTEGER :: num_en_unconverged
485 : CHARACTER(LEN=10) :: bse_davidson_abort_cond_string
486 : REAL(KIND=dp) :: eps_exc_en
487 : CHARACTER(LEN=10) :: success_abort_string
488 : REAL(KIND=dp) :: z_space_energy_cutoff
489 :
490 : CHARACTER(LEN=*), PARAMETER :: routineN = 'success_message'
491 :
492 : CHARACTER(LEN=10) :: multiplet
493 : INTEGER :: handle, i
494 : REAL(KIND=dp) :: alpha
495 :
496 0 : CALL timeset(routineN, handle)
497 :
498 : !Prepare variables for printing
499 0 : SELECT CASE (bse_spin_config)
500 : CASE (bse_singlet)
501 0 : alpha = 2.0_dp
502 0 : multiplet = "Singlet"
503 : CASE (bse_triplet)
504 0 : alpha = 0.0_dp
505 0 : multiplet = "Triplet"
506 : END SELECT
507 :
508 0 : IF (unit_nr > 0) THEN
509 0 : WRITE (unit_nr, *) ' '
510 0 : WRITE (unit_nr, '(T3,A)') '******************************************************************************'
511 0 : WRITE (unit_nr, '(T3,A)') '** **'
512 0 : WRITE (unit_nr, '(T3,A)') '** BSE-TDA EXCITONIC ENERGIES **'
513 0 : WRITE (unit_nr, '(T3,A)') '** **'
514 0 : WRITE (unit_nr, '(T3,A)') '******************************************************************************'
515 0 : WRITE (unit_nr, '(T3,A)') ' '
516 0 : WRITE (unit_nr, '(T3,A)') ' '
517 0 : WRITE (unit_nr, '(T3,A)') ' The excitation energies are calculated by iteratively diagonalizing: '
518 0 : WRITE (unit_nr, '(T3,A)') ' '
519 0 : WRITE (unit_nr, '(T3,A)') ' A_iajb = (E_a-E_i) delta_ij delta_ab + alpha * v_iajb - W_ijab '
520 0 : WRITE (unit_nr, '(T3,A)') ' '
521 : WRITE (unit_nr, '(T3,A48,A7,A12,F3.1)') &
522 0 : ' The spin-dependent factor for the requested ', multiplet, " is alpha = ", alpha
523 0 : WRITE (unit_nr, '(T3,A)') ' '
524 : WRITE (unit_nr, '(T3,A16,T50,A22)') &
525 0 : ' Excitonic level', 'Excitation energy (eV)'
526 : !prints actual energies values
527 0 : DO i = 1, num_exc_en
528 0 : WRITE (unit_nr, '(T3,I16,T50,F22.3)') i, Subspace_full_eigenval(i)*evolt
529 : END DO
530 :
531 0 : WRITE (unit_nr, '(T3,A)') ' '
532 :
533 : !prints parameters of Davidson algorithm
534 : CALL print_davidson_parameter(i_iter, num_davidson_iter, num_Z_vectors, num_res_unconverged, max_res_norm, &
535 : eps_res, num_en_unconverged, max_en_diff, eps_exc_en, num_exc_en, &
536 : num_Z_vectors_init, num_max_z_space, num_new_t, unit_nr, &
537 0 : success_abort_string, bse_davidson_abort_cond_string, z_space_energy_cutoff)
538 :
539 : !Insert warning if energies are not converged (could probably be the case if one uses residual threshold)
540 0 : IF (num_en_unconverged > 0) THEN
541 0 : WRITE (unit_nr, '(T3,A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
542 0 : WRITE (unit_nr, '(T3,A2,T79,A2)') '!!', "!!"
543 0 : WRITE (unit_nr, '(T3,A2,T8,A65,T79,A2)') '!!', "THERE ARE UNCONVERGED ENERGIES PRINTED OUT, SOMETHING WENT WRONG!", "!!"
544 0 : WRITE (unit_nr, '(T3,A2,T79,A2)') '!!', "!!"
545 0 : WRITE (unit_nr, '(T3,A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
546 : END IF
547 : END IF
548 :
549 0 : CALL timestop(handle)
550 :
551 0 : END SUBROUTINE
552 :
553 : ! **************************************************************************************************
554 : !> \brief ...
555 : !> \param i_iter ...
556 : !> \param num_davidson_iter ...
557 : !> \param num_Z_vectors ...
558 : !> \param num_res_unconverged ...
559 : !> \param max_res_norm ...
560 : !> \param eps_res ...
561 : !> \param num_en_unconverged ...
562 : !> \param max_en_diff ...
563 : !> \param eps_exc_en ...
564 : !> \param num_exc_en ...
565 : !> \param num_Z_vectors_init ...
566 : !> \param num_max_z_space ...
567 : !> \param num_new_t ...
568 : !> \param unit_nr ...
569 : !> \param success_abort_string ...
570 : !> \param bse_davidson_abort_cond_string ...
571 : !> \param z_space_energy_cutoff ...
572 : ! **************************************************************************************************
573 0 : SUBROUTINE print_davidson_parameter(i_iter, num_davidson_iter, num_Z_vectors, num_res_unconverged, max_res_norm, &
574 : eps_res, num_en_unconverged, max_en_diff, eps_exc_en, num_exc_en, &
575 : num_Z_vectors_init, num_max_z_space, num_new_t, unit_nr, &
576 : success_abort_string, bse_davidson_abort_cond_string, z_space_energy_cutoff)
577 :
578 : INTEGER :: i_iter, num_davidson_iter, &
579 : num_Z_vectors, num_res_unconverged
580 : REAL(KIND=dp) :: max_res_norm, eps_res
581 : INTEGER :: num_en_unconverged
582 : REAL(KIND=dp) :: max_en_diff, eps_exc_en
583 : INTEGER :: num_exc_en, num_Z_vectors_init, &
584 : num_max_z_space, num_new_t, unit_nr
585 : CHARACTER(LEN=10) :: success_abort_string, &
586 : bse_davidson_abort_cond_string
587 : REAL(KIND=dp) :: z_space_energy_cutoff
588 :
589 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_davidson_parameter'
590 :
591 : INTEGER :: handle
592 :
593 0 : CALL timeset(routineN, handle)
594 :
595 0 : WRITE (unit_nr, '(T3,A)') '******************************************************************************'
596 : WRITE (unit_nr, '(T3,A2,T15,A49,T79,A2)') &
597 0 : '**', "Parameters of the BSE-Davidson solver:", "**"
598 : WRITE (unit_nr, '(T3,A2,T79,A2)') &
599 0 : '**', "**"
600 : WRITE (unit_nr, '(T3,A2,T79,A2)') &
601 0 : '**', "**"
602 : WRITE (unit_nr, '(T3,A2,T10,A16,I5,A12,I5,A8,T79,A2)') &
603 0 : '**', "Converged after ", i_iter, " of maximal ", num_davidson_iter, " cycles,", "**"
604 : WRITE (unit_nr, '(T3,A2,T20,A11,A9,A7,A8,A20,T79,A2)') &
605 0 : '**', "because of ", success_abort_string, " using ", &
606 0 : bse_davidson_abort_cond_string, " threshold condition", "**"
607 : WRITE (unit_nr, '(T3,A2,T79,A2)') &
608 0 : '**', "**"
609 : WRITE (unit_nr, '(T3,A2,T10,A32,T65,I11,T79,A2)') &
610 0 : '**', "The Z space has at the end dim. ", num_Z_vectors, "**"
611 : WRITE (unit_nr, '(T3,A2,T10,A45,T65,I11,T79,A2)') &
612 0 : '**', "Number of unconverged residuals in subspace: ", num_res_unconverged, "**"
613 : WRITE (unit_nr, '(T3,A2,T10,A35,T65,E11.4,T79,A2)') &
614 0 : '**', "largest unconverged residual (eV): ", max_res_norm*evolt, "**"
615 : WRITE (unit_nr, '(T3,A2,T10,A45,T65,E11.4,T79,A2)') &
616 0 : '**', "threshold for convergence of residuals (eV): ", eps_res*evolt, "**"
617 : WRITE (unit_nr, '(T3,A2,T10,A45,T65,I11,T79,A2)') &
618 0 : '**', "Number of desired, but unconverged energies: ", num_en_unconverged, "**"
619 : WRITE (unit_nr, '(T3,A2,T10,A44,T65,E11.4,T79,A2)') &
620 0 : '**', "largest unconverged energy difference (eV): ", max_en_diff*evolt, "**"
621 : WRITE (unit_nr, '(T3,A2,T10,A44,T65,E11.4,T79,A2)') &
622 0 : '**', "threshold for convergence of energies (eV): ", eps_exc_en*evolt, "**"
623 : WRITE (unit_nr, '(T3,A2,T10,A40,T65,I11,T79,A2)') &
624 0 : '**', "number of computed excitation energies: ", num_exc_en, "**"
625 :
626 0 : IF (z_space_energy_cutoff > 0) THEN
627 : WRITE (unit_nr, '(T3,A2,T10,A37,T65,E11.4,T79,A2)') &
628 0 : '**', "cutoff for excitation energies (eV): ", z_space_energy_cutoff*evolt, "**"
629 : END IF
630 :
631 : WRITE (unit_nr, '(T3,A2,T10,A36,T65,I11,T79,A2)') &
632 0 : '**', "number of Z space at the beginning: ", num_Z_vectors_init, "**"
633 : WRITE (unit_nr, '(T3,A2,T10,A30,T65,I11,T79,A2)') &
634 0 : '**', "maximal dimension of Z space: ", num_max_z_space, "**"
635 : WRITE (unit_nr, '(T3,A2,T10,A31,T65,I11,T79,A2)') &
636 0 : '**', "added directions per iteration: ", num_new_t, "**"
637 : WRITE (unit_nr, '(T3,A2,T79,A2)') &
638 0 : '**', "**"
639 : WRITE (unit_nr, '(T3,A2,T79,A2)') &
640 0 : '**', "**"
641 0 : WRITE (unit_nr, '(T3,A)') '******************************************************************************'
642 0 : WRITE (unit_nr, '(T3,A)') ' '
643 :
644 0 : CALL timestop(handle)
645 :
646 0 : END SUBROUTINE
647 :
648 : ! **************************************************************************************************
649 : !> \brief ...
650 : !> \param Subspace_full_eigenval ...
651 : !> \param Subspace_prev_eigenval ...
652 : !> \param eps_exc_en ...
653 : !> \param num_en_unconverged ...
654 : !> \param num_exc_en ...
655 : !> \param max_en_diff ...
656 : !> \param En_diffs ...
657 : ! **************************************************************************************************
658 0 : SUBROUTINE check_en_convergence(Subspace_full_eigenval, Subspace_prev_eigenval, eps_exc_en, num_en_unconverged, &
659 : num_exc_en, max_en_diff, En_diffs)
660 :
661 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_full_eigenval, &
662 : Subspace_prev_eigenval
663 : REAL(KIND=dp) :: eps_exc_en
664 : INTEGER :: num_en_unconverged, num_exc_en
665 : REAL(KIND=dp) :: max_en_diff
666 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: En_diffs
667 :
668 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_en_convergence'
669 :
670 : INTEGER :: handle, mu_l
671 :
672 0 : CALL timeset(routineN, handle)
673 :
674 0 : num_en_unconverged = 0
675 0 : ALLOCATE (En_diffs(num_exc_en))
676 0 : DO mu_l = 1, num_exc_en
677 0 : En_diffs(mu_l) = ABS(Subspace_full_eigenval(mu_l) - Subspace_prev_eigenval(mu_l))
678 0 : IF (En_diffs(mu_l) > eps_exc_en) num_en_unconverged = num_en_unconverged + 1
679 : END DO
680 0 : max_en_diff = MAXVAL(En_diffs)
681 :
682 0 : CALL timestop(handle)
683 :
684 0 : END SUBROUTINE
685 :
686 : ! **************************************************************************************************
687 : !> \brief ...
688 : !> \param Subspace_residuals_reshaped ...
689 : !> \param num_new_t ...
690 : !> \param eps_res ...
691 : !> \param num_res_unconverged ...
692 : !> \param i_iter ...
693 : !> \param max_res_norm ...
694 : !> \param unit_nr ...
695 : !> \param Res_norms ...
696 : ! **************************************************************************************************
697 0 : SUBROUTINE check_res_convergence(Subspace_residuals_reshaped, num_new_t, eps_res, num_res_unconverged, &
698 : i_iter, max_res_norm, unit_nr, Res_norms)
699 :
700 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_residuals_reshaped
701 : INTEGER :: num_new_t
702 : REAL(KIND=dp) :: eps_res
703 : INTEGER :: num_res_unconverged, i_iter
704 : REAL(KIND=dp) :: max_res_norm
705 : INTEGER :: unit_nr
706 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Res_norms
707 :
708 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_res_convergence'
709 :
710 : INTEGER :: handle, mu_l
711 :
712 0 : CALL timeset(routineN, handle)
713 :
714 0 : num_res_unconverged = 0
715 0 : ALLOCATE (Res_norms(num_new_t))
716 0 : DO mu_l = 1, num_new_t
717 0 : Res_norms(mu_l) = NORM2(Subspace_residuals_reshaped(:, mu_l))
718 0 : IF (Res_norms(mu_l) > eps_res) THEN
719 0 : num_res_unconverged = num_res_unconverged + 1
720 0 : IF (unit_nr > 0) THEN
721 0 : WRITE (unit_nr, *) "Unconverged res in i_iter=", i_iter, "is:", Res_norms(mu_l)
722 : END IF
723 : END IF
724 : END DO
725 0 : max_res_norm = MAXVAL(Res_norms)
726 0 : IF (unit_nr > 0) THEN
727 0 : WRITE (unit_nr, *) "Maximal unconverged res (of ", num_res_unconverged, &
728 0 : " unconverged res in this step) in i_iter=", i_iter, "is:", max_res_norm
729 : END IF
730 :
731 0 : CALL timestop(handle)
732 :
733 0 : END SUBROUTINE
734 :
735 : ! **************************************************************************************************
736 : !> \brief ...
737 : !> \param W_vectors ...
738 : !> \param num_Z_vectors ...
739 : !> \param homo ...
740 : !> \param virtual ...
741 : ! **************************************************************************************************
742 0 : SUBROUTINE orthonormalize_W(W_vectors, num_Z_vectors, homo, virtual)
743 :
744 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: W_vectors
745 : INTEGER :: num_Z_vectors, homo, virtual
746 :
747 : CHARACTER(LEN=*), PARAMETER :: routineN = 'orthonormalize_W'
748 :
749 : INTEGER :: handle, info_dor, info_orth, LWORK_dor, &
750 : LWORK_W
751 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Tau_W, WORK_W, WORK_W_dor
752 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: W_vectors_reshaped
753 :
754 0 : CALL timeset(routineN, handle)
755 :
756 0 : ALLOCATE (W_vectors_reshaped(homo*virtual, num_Z_vectors))
757 0 : W_vectors_reshaped(:, :) = RESHAPE(W_vectors, (/homo*virtual, num_Z_vectors/))
758 :
759 0 : ALLOCATE (Tau_W(MIN(homo*virtual, num_Z_vectors)))
760 0 : Tau_W = 0.0_dp
761 :
762 0 : ALLOCATE (WORK_W(1))
763 0 : WORK_W = 0.0_dp
764 :
765 0 : ALLOCATE (WORK_W_dor(1))
766 0 : WORK_W_dor = 0.0_dp
767 :
768 0 : CALL DGEQRF(homo*virtual, num_Z_vectors, W_vectors_reshaped, homo*virtual, Tau_W, WORK_W, -1, info_orth)
769 0 : LWORK_W = INT(WORK_W(1))
770 0 : DEALLOCATE (WORK_W)
771 0 : ALLOCATE (WORK_W(LWORK_W))
772 0 : WORK_W = 0.0_dp
773 0 : CALL DGEQRF(homo*virtual, num_Z_vectors, W_vectors_reshaped, homo*virtual, Tau_W, WORK_W, LWORK_W, info_orth)
774 0 : IF (info_orth /= 0) THEN
775 0 : CPABORT("QR Decomp Step 1 doesnt work")
776 : END IF
777 : CALL DORGQR(homo*virtual, num_Z_vectors, MIN(homo*virtual, num_Z_vectors), W_vectors_reshaped, homo*virtual, &
778 0 : Tau_W, WORK_W_dor, -1, info_dor)
779 0 : LWORK_dor = INT(WORK_W_dor(1))
780 0 : DEALLOCATE (WORK_W_dor)
781 0 : ALLOCATE (WORK_W_dor(LWORK_dor))
782 0 : WORK_W_dor = 0.0_dp
783 : CALL DORGQR(homo*virtual, num_Z_vectors, MIN(homo*virtual, num_Z_vectors), W_vectors_reshaped, homo*virtual, &
784 0 : Tau_W, WORK_W_dor, LWORK_dor, info_dor)
785 0 : IF (info_orth /= 0) THEN
786 0 : CPABORT("QR Decomp Step 2 doesnt work")
787 : END IF
788 :
789 0 : W_vectors(:, :, :) = RESHAPE(W_vectors_reshaped, (/homo, virtual, num_Z_vectors/))
790 :
791 0 : DEALLOCATE (WORK_W, WORK_W_dor, Tau_W, W_vectors_reshaped)
792 :
793 0 : CALL timestop(handle)
794 :
795 0 : END SUBROUTINE
796 :
797 : ! **************************************************************************************************
798 : !> \brief ...
799 : !> \param homo ...
800 : !> \param virtual ...
801 : !> \param Subspace_residuals_reshaped ...
802 : !> \param Subspace_new_eigenval ...
803 : !> \param Eigenval ...
804 : !> \param num_new_t ...
805 : !> \param Subspace_add_dir ...
806 : ! **************************************************************************************************
807 0 : SUBROUTINE compute_new_directions(homo, virtual, Subspace_residuals_reshaped, Subspace_new_eigenval, Eigenval, &
808 : num_new_t, Subspace_add_dir)
809 :
810 : INTEGER :: homo, virtual
811 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_residuals_reshaped
812 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_new_eigenval
813 : REAL(KIND=dp), DIMENSION(:) :: Eigenval
814 : INTEGER :: num_new_t
815 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Subspace_add_dir
816 :
817 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_new_directions'
818 :
819 : INTEGER :: a_virt, handle, i_occ, mu_subspace, &
820 : prec_neg
821 : REAL(KIND=dp) :: prec_scalar
822 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_add_dir_reshaped
823 :
824 0 : CALL timeset(routineN, handle)
825 :
826 0 : ALLOCATE (Subspace_add_dir_reshaped(homo*virtual, num_new_t))
827 :
828 0 : prec_neg = 0
829 0 : DO mu_subspace = 1, num_new_t
830 0 : DO i_occ = 1, homo
831 0 : DO a_virt = 1, virtual
832 : !MG to check: Indexorder and range of indices
833 0 : prec_scalar = -1/(Subspace_new_eigenval(mu_subspace) - (Eigenval(a_virt + homo) - Eigenval(i_occ)))
834 : IF (prec_scalar < 0) THEN
835 0 : prec_neg = prec_neg + 1
836 : !prec_scalar = - prec_scalar
837 : END IF
838 : Subspace_add_dir_reshaped((i_occ - 1)*virtual + a_virt, mu_subspace) = prec_scalar* &
839 0 : Subspace_residuals_reshaped((i_occ - 1)*virtual + a_virt, mu_subspace)
840 : END DO
841 : END DO
842 : END DO
843 :
844 0 : Subspace_add_dir(:, :, :) = RESHAPE(Subspace_add_dir_reshaped, (/homo, virtual, num_new_t/))
845 :
846 0 : DEALLOCATE (Subspace_add_dir_reshaped)
847 0 : CALL timestop(handle)
848 :
849 0 : END SUBROUTINE
850 :
851 : ! **************************************************************************************************
852 : !> \brief ...
853 : !> \param AZ_reshaped ...
854 : !> \param Z_vectors_reshaped ...
855 : !> \param Subspace_new_eigenval ...
856 : !> \param Subspace_new_eigenvec ...
857 : !> \param Subspace_residuals_reshaped ...
858 : !> \param homo ...
859 : !> \param virtual ...
860 : !> \param num_new_t ...
861 : !> \param num_Z_vectors ...
862 : !> \param Subspace_ritzvec ...
863 : ! **************************************************************************************************
864 0 : SUBROUTINE compute_residuals(AZ_reshaped, Z_vectors_reshaped, Subspace_new_eigenval, Subspace_new_eigenvec, &
865 : Subspace_residuals_reshaped, homo, virtual, num_new_t, num_Z_vectors, Subspace_ritzvec)
866 :
867 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: AZ_reshaped, Z_vectors_reshaped
868 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_new_eigenval
869 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_new_eigenvec, &
870 : Subspace_residuals_reshaped
871 : INTEGER :: homo, virtual, num_new_t, num_Z_vectors
872 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Subspace_ritzvec
873 :
874 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_residuals'
875 :
876 : INTEGER :: handle, mu_subspace
877 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_res_A, Subspace_res_ev
878 :
879 0 : CALL timeset(routineN, handle)
880 :
881 0 : ALLOCATE (Subspace_res_ev(homo*virtual, num_new_t))
882 0 : Subspace_res_ev = 0.0_dp
883 :
884 0 : ALLOCATE (Subspace_res_A(homo*virtual, num_new_t))
885 0 : Subspace_res_A = 0.0_dp
886 :
887 : !Compute all residuals in one loop, iterating over number of new/added t per iteration
888 0 : DO mu_subspace = 1, num_new_t
889 :
890 : CALL DGEMM("N", "N", homo*virtual, 1, num_Z_vectors, 1.0_dp, Z_vectors_reshaped, homo*virtual, &
891 0 : Subspace_new_eigenvec(:, mu_subspace), num_Z_vectors, 0.0_dp, Subspace_res_ev(:, mu_subspace), homo*virtual)
892 :
893 : CALL DGEMM("N", "N", homo*virtual, 1, num_Z_vectors, 1.0_dp, AZ_reshaped, homo*virtual, &
894 0 : Subspace_new_eigenvec(:, mu_subspace), num_Z_vectors, 0.0_dp, Subspace_res_A(:, mu_subspace), homo*virtual)
895 :
896 : Subspace_residuals_reshaped(:, mu_subspace) = Subspace_new_eigenval(mu_subspace)*Subspace_res_ev(:, mu_subspace) &
897 0 : - Subspace_res_A(:, mu_subspace)
898 :
899 : END DO
900 0 : Subspace_ritzvec(:, :, :) = RESHAPE(Subspace_res_ev, (/homo, virtual, num_new_t/))
901 0 : DEALLOCATE (Subspace_res_ev, Subspace_res_A)
902 :
903 0 : CALL timestop(handle)
904 :
905 0 : END SUBROUTINE
906 :
907 : ! **************************************************************************************************
908 : !> \brief ...
909 : !> \param AZ_reshaped ...
910 : !> \param Z_vectors_reshaped ...
911 : !> \param num_Z_vectors ...
912 : !> \param Subspace_new_eigenval ...
913 : !> \param Subspace_new_eigenvec ...
914 : !> \param num_new_t ...
915 : !> \param Subspace_full_eigenval ...
916 : !> \param para_env ...
917 : !> \param unit_nr ...
918 : ! **************************************************************************************************
919 0 : SUBROUTINE compute_diagonalize_ZAZ(AZ_reshaped, Z_vectors_reshaped, num_Z_vectors, Subspace_new_eigenval, &
920 : Subspace_new_eigenvec, num_new_t, Subspace_full_eigenval, para_env, unit_nr)
921 :
922 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: AZ_reshaped, Z_vectors_reshaped
923 : INTEGER, INTENT(in) :: num_Z_vectors
924 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_new_eigenval
925 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_new_eigenvec
926 : INTEGER, INTENT(in) :: num_new_t
927 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_full_eigenval
928 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
929 : INTEGER, INTENT(in) :: unit_nr
930 :
931 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_diagonalize_ZAZ'
932 :
933 : INTEGER :: handle, i_Z_vector, j_Z_vector, LWORK, &
934 : ZAZ_diag_info
935 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: WORK
936 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ZAZ
937 :
938 0 : CALL timeset(routineN, handle)
939 :
940 0 : ALLOCATE (ZAZ(num_Z_vectors, num_Z_vectors))
941 0 : ZAZ(:, :) = 0.0_dp
942 :
943 : !Flatten AZ and Z matrices of a certain j_Z_vector w.r.t. occ and virt indices
944 : !Multiply for each j_Z_vec and write into matrix of dim (num_Z_vec, num_Z_vec)
945 0 : DO i_Z_vector = 1, num_Z_vectors
946 0 : DO j_Z_vector = 1, num_Z_vectors
947 0 : ZAZ(j_Z_vector, i_Z_vector) = DOT_PRODUCT(Z_vectors_reshaped(:, j_Z_vector), AZ_reshaped(:, i_Z_vector))
948 : END DO
949 : END DO
950 0 : IF (unit_nr > 0) THEN
951 0 : WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "Before Diag"
952 : END IF
953 :
954 : !MG to do: Check for symmetry of ZAZ!
955 0 : ALLOCATE (WORK(1))
956 0 : WORK = 0.0_dp
957 0 : CALL DSYEV("V", "U", num_Z_vectors, ZAZ, num_Z_vectors, Subspace_full_eigenval, WORK, -1, ZAZ_diag_info)
958 0 : LWORK = INT(WORK(1))
959 0 : DEALLOCATE (WORK)
960 0 : ALLOCATE (WORK(LWORK))
961 0 : WORK = 0.0_dp
962 : !MG to check: Usage of symmetric routine okay? (Correct LWORK?)
963 0 : CALL DSYEV("V", "U", num_Z_vectors, ZAZ, num_Z_vectors, Subspace_full_eigenval, WORK, LWORK, ZAZ_diag_info)
964 :
965 0 : IF (ZAZ_diag_info /= 0) THEN
966 0 : CPABORT("ZAZ could not be diagonalized successfully.")
967 : END IF
968 :
969 0 : IF (unit_nr > 0) THEN
970 0 : WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "After Diag"
971 : END IF
972 :
973 0 : Subspace_new_eigenval(1:num_new_t) = Subspace_full_eigenval(1:num_new_t)
974 0 : Subspace_new_eigenvec(:, 1:num_new_t) = ZAZ(:, 1:num_new_t)
975 0 : DEALLOCATE (WORK)
976 0 : DEALLOCATE (ZAZ)
977 :
978 0 : CALL timestop(handle)
979 :
980 0 : END SUBROUTINE
981 :
982 : ! **************************************************************************************************
983 : !> \brief ...
984 : !> \param BZ ...
985 : !> \param Z_vectors ...
986 : !> \param B_iaQ_bse_local ...
987 : !> \param B_bar_iaQ_bse_local ...
988 : !> \param M_ji_tmp ...
989 : !> \param homo ...
990 : !> \param virtual ...
991 : !> \param num_Z_vectors ...
992 : !> \param local_RI_size ...
993 : !> \param para_env ...
994 : ! **************************************************************************************************
995 0 : SUBROUTINE compute_BZ(BZ, Z_vectors, B_iaQ_bse_local, B_bar_iaQ_bse_local, &
996 0 : M_ji_tmp, homo, virtual, num_Z_vectors, local_RI_size, para_env)
997 :
998 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: BZ, Z_vectors, B_iaQ_bse_local, &
999 : B_bar_iaQ_bse_local
1000 : REAL(KIND=dp), DIMENSION(:, :) :: M_ji_tmp
1001 : INTEGER :: homo, virtual, num_Z_vectors, &
1002 : local_RI_size
1003 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1004 :
1005 : INTEGER :: i_Z_vector, LLL
1006 :
1007 0 : BZ(:, :, :) = 0.0_dp
1008 :
1009 : !CALL compute_v_ia_jb_part(BZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, &
1010 : ! num_Z_vectors, homo, virtual)
1011 :
1012 0 : DO i_Z_vector = 1, num_Z_vectors
1013 :
1014 0 : DO LLL = 1, local_RI_size
1015 :
1016 : ! M_ji^P = sum_b Z_jb*B_bi^P
1017 : CALL DGEMM("N", "T", homo, homo, virtual, 1.0_dp, Z_vectors(:, :, i_Z_vector), homo, &
1018 0 : B_iaQ_bse_local(:, :, LLL), homo, 0.0_dp, M_ji_tmp, homo)
1019 : ! (BZ)_ia = sum_jP M_ij^P*B^bar_ja^P
1020 : CALL DGEMM("T", "N", homo, virtual, homo, 1.0_dp, M_ji_tmp, homo, &
1021 0 : B_bar_iaQ_bse_local, homo, 1.0_dp, BZ(:, :, i_Z_vector), homo)
1022 :
1023 : END DO
1024 :
1025 : END DO
1026 :
1027 : ! we make the sum to sum over all RI basis functions
1028 0 : CALL para_env%sum(BZ)
1029 :
1030 0 : END SUBROUTINE
1031 :
1032 : ! **************************************************************************************************
1033 : !> \brief ...
1034 : !> \param AZ ...
1035 : !> \param Z_vectors ...
1036 : !> \param B_iaQ_bse_local ...
1037 : !> \param B_bar_ijQ_bse_local ...
1038 : !> \param B_abQ_bse_local ...
1039 : !> \param M_ia_tmp ...
1040 : !> \param RI_vector ...
1041 : !> \param Eigenval ...
1042 : !> \param homo ...
1043 : !> \param virtual ...
1044 : !> \param num_Z_vectors ...
1045 : !> \param local_RI_size ...
1046 : !> \param para_env ...
1047 : !> \param bse_spin_config ...
1048 : !> \param z_space_energy_cutoff ...
1049 : !> \param i_iter ...
1050 : !> \param bse_full_diag_debug ...
1051 : !> \param Full_exc_spectrum ...
1052 : !> \param unit_nr ...
1053 : ! **************************************************************************************************
1054 0 : SUBROUTINE compute_AZ(AZ, Z_vectors, B_iaQ_bse_local, B_bar_ijQ_bse_local, B_abQ_bse_local, M_ia_tmp, &
1055 0 : RI_vector, Eigenval, homo, virtual, num_Z_vectors, local_RI_size, &
1056 : para_env, bse_spin_config, z_space_energy_cutoff, i_iter, bse_full_diag_debug, &
1057 : Full_exc_spectrum, unit_nr)
1058 :
1059 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: AZ, Z_vectors, B_iaQ_bse_local, &
1060 : B_bar_ijQ_bse_local, B_abQ_bse_local
1061 : REAL(KIND=dp), DIMENSION(:, :) :: M_ia_tmp, RI_vector
1062 : REAL(KIND=dp), DIMENSION(:) :: Eigenval
1063 : INTEGER :: homo, virtual, num_Z_vectors, &
1064 : local_RI_size
1065 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1066 : INTEGER :: bse_spin_config
1067 : REAL(KIND=dp) :: z_space_energy_cutoff
1068 : INTEGER :: i_iter
1069 : LOGICAL :: bse_full_diag_debug
1070 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Full_exc_spectrum
1071 : INTEGER :: unit_nr
1072 :
1073 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_AZ'
1074 :
1075 : INTEGER :: a, a_virt, b, diag_info, handle, i, &
1076 : i_occ, i_Z_vector, j, LLL, LWORK, m, n
1077 : REAL(KIND=dp) :: eigen_diff
1078 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: WORK
1079 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: A_full_reshaped
1080 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: A_full, v_iajb, W_ijab
1081 :
1082 0 : CALL timeset(routineN, handle)
1083 0 : AZ(:, :, :) = 0.0_dp
1084 :
1085 0 : IF (i_iter == 1 .AND. bse_full_diag_debug) THEN
1086 0 : ALLOCATE (W_ijab(homo, homo, virtual, virtual))
1087 0 : ALLOCATE (A_full(homo, virtual, homo, virtual))
1088 0 : ALLOCATE (A_full_reshaped(homo*virtual, homo*virtual))
1089 0 : ALLOCATE (Full_exc_spectrum(homo*virtual))
1090 0 : W_ijab = 0.0_dp
1091 0 : A_full = 0.0_dp
1092 0 : A_full_reshaped = 0.0_dp
1093 0 : Full_exc_spectrum = 0.0_dp
1094 : END IF
1095 :
1096 : CALL compute_v_ia_jb_part(AZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, &
1097 : num_Z_vectors, homo, virtual, bse_spin_config, v_iajb, bse_full_diag_debug, i_iter, &
1098 0 : para_env)
1099 :
1100 0 : DO i_Z_vector = 1, num_Z_vectors
1101 :
1102 : ! JW TO DO: OMP PARALLELIZATION
1103 0 : DO LLL = 1, local_RI_size
1104 :
1105 : ! M_ja^P = sum_b Z_jb*B_ba^P
1106 : CALL DGEMM("N", "N", homo, virtual, virtual, 1.0_dp, Z_vectors(:, :, i_Z_vector), homo, &
1107 0 : B_abQ_bse_local(:, :, LLL), virtual, 0.0_dp, M_ia_tmp, homo)
1108 :
1109 : ! (AZ)_ia = sum_jP B_bar_ij^P*M_ja^P
1110 : CALL DGEMM("N", "N", homo, virtual, homo, -1.0_dp, B_bar_ijQ_bse_local(:, :, LLL), homo, &
1111 0 : M_ia_tmp, homo, 1.0_dp, AZ(:, :, i_Z_vector), homo)
1112 :
1113 : END DO
1114 : END DO
1115 :
1116 0 : IF (i_iter == 1 .AND. bse_full_diag_debug) THEN
1117 0 : W_ijab = 0.0_dp
1118 : !Create screened 4c integrals for check
1119 0 : DO LLL = 1, local_RI_size
1120 0 : DO i = 1, homo
1121 0 : DO j = 1, homo
1122 0 : DO a = 1, virtual
1123 0 : DO b = 1, virtual
1124 0 : W_ijab(i, j, a, b) = W_ijab(i, j, a, b) + B_bar_ijQ_bse_local(i, j, LLL)*B_abQ_bse_local(a, b, LLL)
1125 : END DO
1126 : END DO
1127 : END DO
1128 : END DO
1129 : END DO
1130 : ! we make the mp_sum to sum over all RI basis functions
1131 0 : CALL para_env%sum(W_ijab)
1132 : END IF
1133 :
1134 : ! we make the mp_sum to sum over all RI basis functions
1135 0 : CALL para_env%sum(AZ)
1136 :
1137 : ! add (e_a-e_i)*Z_ia
1138 0 : DO i_occ = 1, homo
1139 0 : DO a_virt = 1, virtual
1140 :
1141 0 : eigen_diff = Eigenval(a_virt + homo) - Eigenval(i_occ)
1142 0 : IF (unit_nr > 0 .AND. i_iter == 1) THEN
1143 0 : WRITE (unit_nr, *) "Ediff at (i_occ,a_virt)=", i_occ, a_virt, " is: ", eigen_diff
1144 : END IF
1145 :
1146 0 : AZ(i_occ, a_virt, :) = AZ(i_occ, a_virt, :) + Z_vectors(i_occ, a_virt, :)*eigen_diff
1147 :
1148 : END DO
1149 : END DO
1150 :
1151 : !cut off contributions, which are too high in the excitation spectrum
1152 0 : IF (z_space_energy_cutoff > 0) THEN
1153 0 : DO i_occ = 1, homo
1154 0 : DO a_virt = 1, virtual
1155 :
1156 0 : IF (Eigenval(a_virt + homo) > z_space_energy_cutoff .OR. -Eigenval(i_occ) > z_space_energy_cutoff) THEN
1157 0 : AZ(i_occ, a_virt, :) = 0
1158 : END IF
1159 :
1160 : END DO
1161 : END DO
1162 : END IF
1163 :
1164 : !Debugging purposes: full diagonalization of A
1165 0 : IF (i_iter == 1 .AND. bse_full_diag_debug) THEN
1166 0 : n = 0
1167 0 : DO i = 1, homo
1168 0 : DO a = 1, virtual
1169 0 : n = n + 1
1170 0 : m = 0
1171 0 : DO j = 1, homo
1172 0 : DO b = 1, virtual
1173 0 : m = m + 1
1174 0 : IF (a == b .AND. i == j) THEN
1175 0 : eigen_diff = Eigenval(a + homo) - Eigenval(i)
1176 : ELSE
1177 0 : eigen_diff = 0
1178 : END IF
1179 0 : A_full_reshaped(n, m) = eigen_diff + 2*v_iajb(i, a, j, b) - W_ijab(i, j, a, b)
1180 0 : A_full(i, a, j, b) = eigen_diff + 2*v_iajb(i, a, j, b) - W_ijab(i, j, a, b)
1181 : END DO
1182 : END DO
1183 : END DO
1184 : END DO
1185 :
1186 : !MG to do: Check for symmetry of ZAZ!
1187 0 : ALLOCATE (WORK(1))
1188 0 : WORK = 0.0_dp
1189 0 : CALL DSYEV("N", "U", homo*virtual, A_full_reshaped, homo*virtual, Full_exc_spectrum, WORK, -1, diag_info)
1190 0 : LWORK = INT(WORK(1))
1191 0 : DEALLOCATE (WORK)
1192 0 : ALLOCATE (WORK(LWORK))
1193 0 : WORK = 0.0_dp
1194 : !MG to check: Usage of symmetric routine okay? (Correct LWORK?)
1195 0 : CALL DSYEV("N", "U", homo*virtual, A_full_reshaped, homo*virtual, Full_exc_spectrum, WORK, LWORK, diag_info)
1196 :
1197 0 : DEALLOCATE (WORK)
1198 :
1199 0 : DEALLOCATE (W_ijab, v_iajb, A_full, A_full_reshaped)
1200 : END IF
1201 :
1202 0 : CALL timestop(handle)
1203 :
1204 0 : END SUBROUTINE
1205 :
1206 : ! **************************************************************************************************
1207 : !> \brief ...
1208 : !> \param AZ ...
1209 : !> \param Z_vectors ...
1210 : !> \param B_iaQ_bse_local ...
1211 : !> \param RI_vector ...
1212 : !> \param local_RI_size ...
1213 : !> \param num_Z_vectors ...
1214 : !> \param homo ...
1215 : !> \param virtual ...
1216 : !> \param bse_spin_config ...
1217 : !> \param v_iajb ...
1218 : !> \param bse_full_diag_debug ...
1219 : !> \param i_iter ...
1220 : !> \param para_env ...
1221 : ! **************************************************************************************************
1222 0 : SUBROUTINE compute_v_ia_jb_part(AZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, &
1223 : num_Z_vectors, homo, virtual, bse_spin_config, v_iajb, bse_full_diag_debug, i_iter, &
1224 : para_env)
1225 :
1226 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: AZ, Z_vectors, B_iaQ_bse_local
1227 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: RI_vector
1228 : INTEGER, INTENT(IN) :: local_RI_size, num_Z_vectors, homo, &
1229 : virtual, bse_spin_config
1230 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: v_iajb
1231 : LOGICAL :: bse_full_diag_debug
1232 : INTEGER, INTENT(IN) :: i_iter
1233 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1234 :
1235 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_v_ia_jb_part'
1236 :
1237 : INTEGER :: a, a_virt, b, handle, i, i_occ, &
1238 : i_Z_vector, j, LLL
1239 : REAL(KIND=dp) :: alpha
1240 :
1241 : !debugging:
1242 :
1243 0 : CALL timeset(routineN, handle)
1244 :
1245 : !Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
1246 0 : SELECT CASE (bse_spin_config)
1247 : CASE (bse_singlet)
1248 0 : alpha = 2.0_dp
1249 : CASE (bse_triplet)
1250 0 : alpha = 0.0_dp
1251 : END SELECT
1252 :
1253 0 : RI_vector = 0.0_dp
1254 :
1255 : ! v_P = sum_jb B_jb^P Z_jb
1256 0 : DO LLL = 1, local_RI_size
1257 0 : DO i_Z_vector = 1, num_Z_vectors
1258 0 : DO i_occ = 1, homo
1259 0 : DO a_virt = 1, virtual
1260 :
1261 : RI_vector(LLL, i_Z_vector) = RI_vector(LLL, i_Z_vector) + &
1262 : Z_vectors(i_occ, a_virt, i_Z_vector)* &
1263 0 : B_iaQ_bse_local(i_occ, a_virt, LLL)
1264 :
1265 : END DO
1266 : END DO
1267 : END DO
1268 : END DO
1269 :
1270 : ! AZ = sum_P B_ia^P*v_P + ...
1271 0 : DO LLL = 1, local_RI_size
1272 0 : DO i_Z_vector = 1, num_Z_vectors
1273 0 : DO i_occ = 1, homo
1274 0 : DO a_virt = 1, virtual
1275 : !MG to check: Minus sign at v oder W? Factor for triplet/singlet
1276 : AZ(i_occ, a_virt, i_Z_vector) = AZ(i_occ, a_virt, i_Z_vector) + &
1277 : alpha*RI_vector(LLL, i_Z_vector)* &
1278 0 : B_iaQ_bse_local(i_occ, a_virt, LLL)
1279 :
1280 : END DO
1281 : END DO
1282 : END DO
1283 : END DO
1284 0 : IF (i_iter == 1 .AND. bse_full_diag_debug) THEN
1285 0 : ALLOCATE (v_iajb(homo, virtual, homo, virtual))
1286 0 : v_iajb = 0.0_dp
1287 : !Create unscreened 4c integrals for check
1288 0 : DO LLL = 1, local_RI_size
1289 0 : DO i = 1, homo
1290 0 : DO j = 1, homo
1291 0 : DO a = 1, virtual
1292 0 : DO b = 1, virtual
1293 0 : v_iajb(i, a, j, b) = v_iajb(i, a, j, b) + B_iaQ_bse_local(i, a, LLL)*B_iaQ_bse_local(j, b, LLL)
1294 : END DO
1295 : END DO
1296 : END DO
1297 : END DO
1298 : END DO
1299 : ! we make the mp_sum to sum over all RI basis functions
1300 0 : CALL para_env%sum(v_iajb)
1301 : END IF
1302 :
1303 0 : CALL timestop(handle)
1304 :
1305 0 : END SUBROUTINE
1306 :
1307 : ! **************************************************************************************************
1308 : !> \brief ...Eigenval
1309 : !> \param Z_vectors ...
1310 : !> \param Eigenval ...
1311 : !> \param num_Z_vectors ...
1312 : !> \param homo ...
1313 : !> \param virtual ...
1314 : ! **************************************************************************************************
1315 0 : SUBROUTINE initial_guess_Z_vectors(Z_vectors, Eigenval, num_Z_vectors, homo, virtual)
1316 :
1317 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: Z_vectors
1318 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
1319 : INTEGER, INTENT(IN) :: num_Z_vectors, homo, virtual
1320 :
1321 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initial_guess_Z_vectors'
1322 :
1323 : INTEGER :: a_virt, handle, i_occ, i_Z_vector, &
1324 : min_loc(2)
1325 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigen_diff_ia
1326 :
1327 0 : CALL timeset(routineN, handle)
1328 :
1329 0 : ALLOCATE (eigen_diff_ia(homo, virtual))
1330 :
1331 0 : DO i_occ = 1, homo
1332 0 : DO a_virt = 1, virtual
1333 0 : eigen_diff_ia(i_occ, a_virt) = Eigenval(a_virt + homo) - Eigenval(i_occ)
1334 : END DO
1335 : END DO
1336 :
1337 0 : DO i_Z_vector = 1, num_Z_vectors
1338 :
1339 0 : min_loc = MINLOC(eigen_diff_ia)
1340 :
1341 0 : Z_vectors(min_loc(1), min_loc(2), i_Z_vector) = 1.0_dp
1342 :
1343 0 : eigen_diff_ia(min_loc(1), min_loc(2)) = 1.0E20_dp
1344 :
1345 : END DO
1346 :
1347 0 : DEALLOCATE (eigen_diff_ia)
1348 :
1349 0 : CALL timestop(handle)
1350 :
1351 0 : END SUBROUTINE
1352 :
1353 : ! **************************************************************************************************
1354 : !> \brief ...
1355 : !> \param fm_mat_S_ab_bse ...
1356 : !> \param fm_mat_S ...
1357 : !> \param fm_mat_S_bar_ia_bse ...
1358 : !> \param fm_mat_S_bar_ij_bse ...
1359 : !> \param B_bar_ijQ_bse_local ...
1360 : !> \param B_abQ_bse_local ...
1361 : !> \param B_bar_iaQ_bse_local ...
1362 : !> \param B_iaQ_bse_local ...
1363 : !> \param dimen_RI ...
1364 : !> \param homo ...
1365 : !> \param virtual ...
1366 : !> \param gd_array ...
1367 : !> \param color_sub ...
1368 : !> \param para_env ...
1369 : ! **************************************************************************************************
1370 0 : SUBROUTINE fill_local_3c_arrays(fm_mat_S_ab_bse, fm_mat_S, &
1371 : fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
1372 : B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
1373 : B_iaQ_bse_local, dimen_RI, homo, virtual, &
1374 : gd_array, color_sub, para_env)
1375 :
1376 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ab_bse, fm_mat_S, &
1377 : fm_mat_S_bar_ia_bse, &
1378 : fm_mat_S_bar_ij_bse
1379 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1380 : INTENT(OUT) :: B_bar_ijQ_bse_local, B_abQ_bse_local, &
1381 : B_bar_iaQ_bse_local, B_iaQ_bse_local
1382 : INTEGER, INTENT(IN) :: dimen_RI, homo, virtual
1383 : TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
1384 : INTEGER, INTENT(IN) :: color_sub
1385 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1386 :
1387 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_local_3c_arrays'
1388 :
1389 : INTEGER :: handle
1390 :
1391 0 : CALL timeset(routineN, handle)
1392 :
1393 0 : CALL allocate_and_fill_local_array(B_iaQ_bse_local, fm_mat_S, gd_array, color_sub, homo, virtual, dimen_RI, para_env)
1394 :
1395 : CALL allocate_and_fill_local_array(B_bar_iaQ_bse_local, fm_mat_S_bar_ia_bse, gd_array, color_sub, homo, virtual, &
1396 0 : dimen_RI, para_env)
1397 :
1398 : CALL allocate_and_fill_local_array(B_bar_ijQ_bse_local, fm_mat_S_bar_ij_bse, gd_array, color_sub, homo, homo, &
1399 0 : dimen_RI, para_env)
1400 :
1401 : CALL allocate_and_fill_local_array(B_abQ_bse_local, fm_mat_S_ab_bse, gd_array, color_sub, virtual, virtual, &
1402 0 : dimen_RI, para_env)
1403 :
1404 0 : CALL timestop(handle)
1405 :
1406 0 : END SUBROUTINE
1407 :
1408 : ! **************************************************************************************************
1409 : !> \brief ...
1410 : !> \param B_local ...
1411 : !> \param fm_mat_S ...
1412 : !> \param gd_array ...
1413 : !> \param color_sub ...
1414 : !> \param small_size ...
1415 : !> \param big_size ...
1416 : !> \param dimen_RI ...
1417 : !> \param para_env ...
1418 : ! **************************************************************************************************
1419 0 : SUBROUTINE allocate_and_fill_local_array(B_local, fm_mat_S, gd_array, &
1420 : color_sub, small_size, big_size, dimen_RI, para_env)
1421 :
1422 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1423 : INTENT(OUT) :: B_local
1424 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
1425 : TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
1426 : INTEGER, INTENT(IN) :: color_sub, small_size, big_size, dimen_RI
1427 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1428 :
1429 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_and_fill_local_array'
1430 :
1431 : INTEGER :: combi_index, end_RI, handle, handle1, i_comm, i_entry, iiB, imepos, jjB, &
1432 : level_big_size, level_small_size, ncol_local, nrow_local, num_comm_cycles, RI_index, &
1433 : size_RI, start_RI
1434 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, mepos_from_RI_index, &
1435 0 : num_entries_rec, num_entries_send
1436 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1437 : REAL(KIND=dp) :: matrix_el
1438 : TYPE(integ_mat_buffer_type), ALLOCATABLE, &
1439 0 : DIMENSION(:) :: buffer_rec, buffer_send
1440 0 : TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
1441 :
1442 0 : CALL timeset(routineN, handle)
1443 :
1444 0 : ALLOCATE (mepos_from_RI_index(dimen_RI))
1445 0 : mepos_from_RI_index = 0
1446 :
1447 0 : DO imepos = 0, para_env%num_pe - 1
1448 :
1449 0 : CALL get_group_dist(gd_array, pos=imepos, starts=start_RI, ends=end_RI)
1450 :
1451 0 : mepos_from_RI_index(start_RI:end_RI) = imepos
1452 :
1453 : END DO
1454 :
1455 : ! color_sub is automatically the number of the process since every subgroup has only one MPI rank
1456 0 : CALL get_group_dist(gd_array, color_sub, start_RI, end_RI, size_RI)
1457 :
1458 0 : ALLOCATE (B_local(small_size, big_size, 1:size_RI))
1459 :
1460 0 : ALLOCATE (num_entries_send(0:para_env%num_pe - 1))
1461 0 : ALLOCATE (num_entries_rec(0:para_env%num_pe - 1))
1462 :
1463 0 : ALLOCATE (req_array(1:para_env%num_pe, 4))
1464 :
1465 0 : ALLOCATE (entry_counter(0:para_env%num_pe - 1))
1466 :
1467 : CALL cp_fm_get_info(matrix=fm_mat_S, &
1468 : nrow_local=nrow_local, &
1469 : ncol_local=ncol_local, &
1470 : row_indices=row_indices, &
1471 0 : col_indices=col_indices)
1472 :
1473 0 : num_comm_cycles = 10
1474 :
1475 : ! communicate not all due to huge memory overhead, since for every number in fm_mat_S, we store
1476 : ! three additional ones (RI index, first MO index, second MO index!!)
1477 0 : DO i_comm = 0, num_comm_cycles - 1
1478 :
1479 0 : num_entries_send = 0
1480 0 : num_entries_rec = 0
1481 :
1482 : ! loop over RI index to get the number of sent entries
1483 0 : DO jjB = 1, nrow_local
1484 :
1485 0 : RI_index = row_indices(jjB)
1486 :
1487 0 : IF (MODULO(RI_index, num_comm_cycles) /= i_comm) CYCLE
1488 :
1489 0 : imepos = mepos_from_RI_index(RI_index)
1490 :
1491 0 : num_entries_send(imepos) = num_entries_send(imepos) + ncol_local
1492 :
1493 : END DO
1494 :
1495 0 : CALL para_env%alltoall(num_entries_send, num_entries_rec, 1)
1496 :
1497 0 : ALLOCATE (buffer_rec(0:para_env%num_pe - 1))
1498 0 : ALLOCATE (buffer_send(0:para_env%num_pe - 1))
1499 :
1500 : ! allocate data message and corresponding indices
1501 0 : DO imepos = 0, para_env%num_pe - 1
1502 :
1503 0 : ALLOCATE (buffer_rec(imepos)%msg(num_entries_rec(imepos)))
1504 0 : buffer_rec(imepos)%msg = 0.0_dp
1505 :
1506 0 : ALLOCATE (buffer_send(imepos)%msg(num_entries_send(imepos)))
1507 0 : buffer_send(imepos)%msg = 0.0_dp
1508 :
1509 0 : ALLOCATE (buffer_rec(imepos)%indx(num_entries_rec(imepos), 3))
1510 0 : buffer_rec(imepos)%indx = 0
1511 :
1512 0 : ALLOCATE (buffer_send(imepos)%indx(num_entries_send(imepos), 3))
1513 0 : buffer_send(imepos)%indx = 0
1514 :
1515 : END DO
1516 :
1517 0 : entry_counter(:) = 0
1518 :
1519 : ! loop over RI index for filling the send-buffer
1520 0 : DO jjB = 1, nrow_local
1521 :
1522 0 : RI_index = row_indices(jjB)
1523 :
1524 0 : IF (MODULO(RI_index, num_comm_cycles) /= i_comm) CYCLE
1525 :
1526 0 : imepos = mepos_from_RI_index(RI_index)
1527 :
1528 0 : DO iiB = 1, ncol_local
1529 :
1530 0 : combi_index = col_indices(iiB)
1531 0 : level_small_size = MAX(1, combi_index - 1)/MAX(big_size, 2) + 1
1532 0 : level_big_size = combi_index - (level_small_size - 1)*big_size
1533 :
1534 0 : entry_counter(imepos) = entry_counter(imepos) + 1
1535 :
1536 0 : buffer_send(imepos)%msg(entry_counter(imepos)) = fm_mat_S%local_data(jjB, iiB)
1537 :
1538 0 : buffer_send(imepos)%indx(entry_counter(imepos), 1) = RI_index
1539 0 : buffer_send(imepos)%indx(entry_counter(imepos), 2) = level_small_size
1540 0 : buffer_send(imepos)%indx(entry_counter(imepos), 3) = level_big_size
1541 :
1542 : END DO
1543 :
1544 : END DO
1545 :
1546 0 : CALL timeset("BSE_comm_data", handle1)
1547 :
1548 0 : CALL communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, buffer_send, req_array)
1549 :
1550 0 : CALL timestop(handle1)
1551 :
1552 : ! fill B_local
1553 0 : DO imepos = 0, para_env%num_pe - 1
1554 :
1555 0 : DO i_entry = 1, num_entries_rec(imepos)
1556 :
1557 0 : RI_index = buffer_rec(imepos)%indx(i_entry, 1) - start_RI + 1
1558 0 : level_small_size = buffer_rec(imepos)%indx(i_entry, 2)
1559 0 : level_big_size = buffer_rec(imepos)%indx(i_entry, 3)
1560 :
1561 0 : matrix_el = buffer_rec(imepos)%msg(i_entry)
1562 :
1563 0 : B_local(level_small_size, level_big_size, RI_index) = matrix_el
1564 :
1565 : END DO
1566 :
1567 : END DO
1568 :
1569 0 : DO imepos = 0, para_env%num_pe - 1
1570 0 : DEALLOCATE (buffer_send(imepos)%msg)
1571 0 : DEALLOCATE (buffer_send(imepos)%indx)
1572 0 : DEALLOCATE (buffer_rec(imepos)%msg)
1573 0 : DEALLOCATE (buffer_rec(imepos)%indx)
1574 : END DO
1575 :
1576 0 : DEALLOCATE (buffer_rec, buffer_send)
1577 :
1578 : END DO
1579 :
1580 0 : DEALLOCATE (num_entries_send, num_entries_rec)
1581 :
1582 0 : DEALLOCATE (mepos_from_RI_index)
1583 :
1584 0 : DEALLOCATE (entry_counter, req_array)
1585 :
1586 0 : CALL timestop(handle)
1587 :
1588 0 : END SUBROUTINE
1589 :
1590 : END MODULE bse_iterative
|
