(git:ed6f26b)
Loading...
Searching...
No Matches
bse_iterative.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 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! **************************************************************************************************
14MODULE 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
38CONTAINS
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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: en_diffs, en_diffs_exact, full_exc_spectrum, &
77 res_norms, subspace_full_eigenval, subspace_new_eigenval, subspace_prev_eigenval
78 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: az_reshaped, m_ia_tmp, m_ji_tmp, ri_vector, &
79 subspace_new_eigenvec, subspace_residuals_reshaped, z_vectors_reshaped
80 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: az, bz, subspace_add_dir, &
81 subspace_ritzvec, w_vectors, z_vectors
82
83 CALL timeset(routinen, handle)
84
85 !MG to del
86 !Debug flag for exact diagonalization (only using lapack!!!)
87 bse_full_diag_debug = .true.
88 num_en_unconverged = -1
89 num_res_unconverged = -1
90 num_exact_en_unconverged = -1
91
92 bse_davidson_abort_cond = mp2_env%bse%davidson_abort_cond
93 num_exc_en = mp2_env%bse%num_exc_en
94 num_add_start_z_space = mp2_env%bse%num_add_start_z_space
95 fac_max_z_space = mp2_env%bse%fac_max_z_space
96 num_new_t = mp2_env%bse%num_new_t
97 num_davidson_iter = mp2_env%bse%num_davidson_iter
98 eps_res = mp2_env%bse%eps_res
99 eps_exc_en = mp2_env%bse%eps_exc_en
100 z_space_energy_cutoff = mp2_env%bse%z_space_energy_cutoff
101
102 num_z_vectors_init = num_exc_en + num_add_start_z_space
103
104 IF (unit_nr > 0) THEN
105 WRITE (unit_nr, *) "bse_spin_config", bse_spin_config
106 WRITE (unit_nr, *) "num_exc_en", num_exc_en
107 WRITE (unit_nr, *) "num_add_start_z_space", num_add_start_z_space
108 WRITE (unit_nr, *) "num_Z_vectors_init", num_z_vectors_init
109 WRITE (unit_nr, *) "fac_max_z_space", fac_max_z_space
110 WRITE (unit_nr, *) "num_new_t", num_new_t
111 WRITE (unit_nr, *) "eps_res", eps_res
112 WRITE (unit_nr, *) "num_davidson_iter", num_davidson_iter
113 WRITE (unit_nr, *) "eps_exc_en", eps_exc_en
114 WRITE (unit_nr, *) "bse_davidson_abort_cond", bse_davidson_abort_cond
115 WRITE (unit_nr, *) "z_space_energy_cutoff", z_space_energy_cutoff
116 WRITE (unit_nr, *) "Printing B_bar_iaQ_bse_local of shape", shape(b_bar_iaq_bse_local)
117 END IF
118
119 local_ri_size = SIZE(b_iaq_bse_local, 3)
120
121 num_z_vectors = num_z_vectors_init
122 num_max_z_space = num_z_vectors_init*fac_max_z_space
123
124 !Check input parameters and correct them if necessary
125 IF (num_new_t > num_z_vectors_init) THEN
126 num_new_t = num_z_vectors_init
127 IF (unit_nr > 0) THEN
128 CALL cp_warn(__location__, "Number of added directions has to be smaller/equals than "// &
129 "initial dimension. Corrected num_new_t accordingly.")
130 END IF
131 END IF
132 IF (unit_nr > 0) THEN
133 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 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 "Choose another num_exc_en and num_add_start_z_space or adapt your basis set.")
139 END IF
140 IF (num_max_z_space .GE. homo*virtual) THEN
141 fac_max_z_space = homo*virtual/num_z_vectors_init
142 num_max_z_space = num_z_vectors_init*fac_max_z_space
143
144 IF (fac_max_z_space == 0) THEN
145 CALL cp_abort(__location__, "Maximal dimension was too large and could not be corrected. "// &
146 "Choose another fac_max_z_space and num_Z_vectors_init or adapt your basis set.")
147 ELSE
148 IF (unit_nr > 0) THEN
149 CALL cp_warn(__location__, "Maximal dimension of Z space has to be smaller than homo*virtual. "// &
150 "Corrected fac_max_z_space accordingly.")
151 END IF
152 END IF
153 END IF
154
155 DO i_iter = 1, num_davidson_iter
156 IF (unit_nr > 0) THEN
157 WRITE (unit_nr, *) "Allocating Z_vec,AZ,BZ with dimensions (homo,virt,num_Z)", homo, virtual, num_z_vectors
158 WRITE (unit_nr, *) 'ProcNr', para_env%mepos, 'you really enter here for i_iter', i_iter
159 END IF
160 ALLOCATE (z_vectors(homo, virtual, num_z_vectors))
161 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 IF (i_iter == 1) THEN
166 CALL initial_guess_z_vectors(z_vectors, eigenval, num_z_vectors, homo, virtual)
167 ALLOCATE (subspace_prev_eigenval(num_exc_en))
168 subspace_prev_eigenval = 0.0_dp
169 ELSE
170 z_vectors(:, :, :) = w_vectors(:, :, :)
171 DEALLOCATE (w_vectors)
172 END IF
173 IF (unit_nr > 0) THEN
174 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 homo, virtual, num_z_vectors, local_ri_size, num_new_t)
181 IF (unit_nr > 0) THEN
182 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 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 IF (unit_nr > 0) THEN
196 WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "Computed AZ"
197 END IF
198
199 !MG to check: Reshaping correct?
200 az_reshaped(:, :) = reshape(az, (/homo*virtual, num_z_vectors/))
201 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 subspace_new_eigenvec, num_new_t, subspace_full_eigenval, para_env, unit_nr)
206 IF (unit_nr > 0) THEN
207 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 num_exc_en, max_en_diff, en_diffs)
213 IF (unit_nr > 0) THEN
214 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 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 i_iter, max_res_norm, unit_nr, res_norms)
223
224 davidson_converged = -1
225 IF (num_res_unconverged == 0 .AND. bse_davidson_abort_cond /= 0) THEN
226 davidson_converged = 1
227 success_abort_string = "RESIDUALS"
228 ELSE IF (num_en_unconverged == 0 .AND. (bse_davidson_abort_cond /= 1)) THEN
229 davidson_converged = 1
230 success_abort_string = "ENERGIES"
231 ELSE IF (i_iter == num_davidson_iter) THEN
232 davidson_converged = -100
233 success_abort_string = "-----"
234 ELSE
235 davidson_converged = -1
236 END IF
237
238 IF (bse_davidson_abort_cond == 0) THEN
239 bse_davidson_abort_cond_string = "ENERGY"
240 ELSE IF (bse_davidson_abort_cond == 1) THEN
241 bse_davidson_abort_cond_string = "RESIDUAL"
242 ELSE
243 bse_davidson_abort_cond_string = "EITHER"
244 END IF
245
246 IF (davidson_converged == 1) THEN
247 CALL postprocess_bse(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 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 DEALLOCATE (az, bz, &
255 z_vectors, m_ia_tmp, m_ji_tmp, ri_vector, subspace_prev_eigenval, &
256 subspace_new_eigenval, subspace_new_eigenvec, subspace_residuals_reshaped, &
257 subspace_add_dir, az_reshaped, z_vectors_reshaped, subspace_ritzvec, subspace_full_eigenval)
258
259 EXIT
260 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 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 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 num_new_t, subspace_add_dir)
273
274 !If exact-diag: compute difference to exact eigenvalues
275 IF (bse_full_diag_debug) THEN
276 ALLOCATE (en_diffs_exact(num_exc_en))
277 num_exact_en_unconverged = 0
278 DO j_print = 1, num_exc_en
279 en_diffs_exact(j_print) = abs(subspace_full_eigenval(j_print) - full_exc_spectrum(j_print))
280 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 num_z_vectors, num_new_t, num_max_z_space, homo, virtual, i_iter, unit_nr)
287
288 !Copy eigenvalues for threshold
289 subspace_prev_eigenval(:) = subspace_full_eigenval(:num_exc_en)
290
291 DEALLOCATE (az, & !BZ,
292 z_vectors, m_ia_tmp, m_ji_tmp, ri_vector, &
293 subspace_new_eigenval, subspace_new_eigenvec, subspace_residuals_reshaped, &
294 subspace_add_dir, az_reshaped, z_vectors_reshaped, subspace_ritzvec, subspace_full_eigenval, &
295 res_norms, en_diffs)
296
297 IF (bse_full_diag_debug) THEN
298 DEALLOCATE (en_diffs_exact)
299 END IF
300
301 !Orthonorm:
302 CALL orthonormalize_w(w_vectors, num_z_vectors, homo, virtual)
303
304 END DO
305
306 CALL timestop(handle)
307
308 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 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 CALL timeset(routinen, handle)
338
339 IF (num_z_vectors + num_new_t .LE. num_max_z_space) THEN
340 w_vectors(:, :, :num_z_vectors) = z_vectors(:, :, :)
341 w_vectors(:, :, num_z_vectors + 1:) = subspace_add_dir
342 num_z_vectors = num_z_vectors + num_new_t
343 ELSE
344 IF (unit_nr > 0) THEN
345 WRITE (unit_nr, *) "Resetting dimension in i_iter=", i_iter
346 END IF
347 DEALLOCATE (w_vectors)
348 ALLOCATE (w_vectors(homo, virtual, 2*num_new_t))
349 w_vectors(:, :, :num_new_t) = subspace_ritzvec(:, :, :)
350 w_vectors(:, :, num_new_t + 1:) = subspace_add_dir
351 num_z_vectors = 2*num_new_t
352 END IF
353
354 CALL timestop(handle)
355
356 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 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 CALL timeset(routinen, handle)
403
404 ALLOCATE (az(homo, virtual, num_z_vectors))
405 az = 0.0_dp
406
407 ALLOCATE (z_vectors_reshaped(homo*virtual, num_z_vectors))
408 z_vectors_reshaped = 0.0_dp
409
410 ALLOCATE (az_reshaped(homo*virtual, num_z_vectors))
411 az_reshaped = 0.0_dp
412
413 ALLOCATE (bz(homo, virtual, num_z_vectors))
414 bz = 0.0_dp
415
416 ALLOCATE (m_ia_tmp(homo, virtual))
417 m_ia_tmp = 0.0_dp
418
419 ALLOCATE (m_ji_tmp(homo, homo))
420 m_ji_tmp = 0.0_dp
421
422 ALLOCATE (ri_vector(local_ri_size, num_z_vectors))
423 ri_vector = 0.0_dp
424
425 ALLOCATE (subspace_new_eigenval(num_new_t))
426 subspace_new_eigenval = 0.0_dp
427
428 ALLOCATE (subspace_full_eigenval(num_z_vectors))
429 subspace_full_eigenval = 0.0_dp
430
431 ALLOCATE (subspace_new_eigenvec(num_z_vectors, num_new_t))
432 subspace_new_eigenvec = 0.0_dp
433
434 ALLOCATE (subspace_residuals_reshaped(homo*virtual, num_new_t))
435 subspace_residuals_reshaped = 0.0_dp
436
437 ALLOCATE (subspace_ritzvec(homo, virtual, num_new_t))
438 subspace_ritzvec = 0.0_dp
439
440 ALLOCATE (subspace_add_dir(homo, virtual, num_new_t))
441 subspace_add_dir = 0.0_dp
442
443 ALLOCATE (w_vectors(homo, virtual, num_z_vectors + num_new_t))
444 w_vectors = 0.0_dp
445
446 CALL timestop(handle)
447
448 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 SUBROUTINE postprocess_bse(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 = 'postprocess_bse'
491
492 CHARACTER(LEN=10) :: multiplet
493 INTEGER :: handle, i
494 REAL(kind=dp) :: alpha
495
496 CALL timeset(routinen, handle)
497
498 !Prepare variables for printing
499 SELECT CASE (bse_spin_config)
500 CASE (bse_singlet)
501 alpha = 2.0_dp
502 multiplet = "Singlet"
503 CASE (bse_triplet)
504 alpha = 0.0_dp
505 multiplet = "Triplet"
506 END SELECT
507
508 IF (unit_nr > 0) THEN
509 WRITE (unit_nr, *) ' '
510 WRITE (unit_nr, '(T3,A)') '******************************************************************************'
511 WRITE (unit_nr, '(T3,A)') '** **'
512 WRITE (unit_nr, '(T3,A)') '** BSE-TDA EXCITONIC ENERGIES **'
513 WRITE (unit_nr, '(T3,A)') '** **'
514 WRITE (unit_nr, '(T3,A)') '******************************************************************************'
515 WRITE (unit_nr, '(T3,A)') ' '
516 WRITE (unit_nr, '(T3,A)') ' '
517 WRITE (unit_nr, '(T3,A)') ' The excitation energies are calculated by iteratively diagonalizing: '
518 WRITE (unit_nr, '(T3,A)') ' '
519 WRITE (unit_nr, '(T3,A)') ' A_iajb = (E_a-E_i) delta_ij delta_ab + alpha * v_iajb - W_ijab '
520 WRITE (unit_nr, '(T3,A)') ' '
521 WRITE (unit_nr, '(T3,A48,A7,A12,F3.1)') &
522 ' The spin-dependent factor for the requested ', multiplet, " is alpha = ", alpha
523 WRITE (unit_nr, '(T3,A)') ' '
524 WRITE (unit_nr, '(T3,A16,T50,A22)') &
525 ' Excitonic level', 'Excitation energy (eV)'
526 !prints actual energies values
527 DO i = 1, num_exc_en
528 WRITE (unit_nr, '(T3,I16,T50,F22.3)') i, subspace_full_eigenval(i)*evolt
529 END DO
530
531 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 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 IF (num_en_unconverged > 0) THEN
541 WRITE (unit_nr, '(T3,A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
542 WRITE (unit_nr, '(T3,A2,T79,A2)') '!!', "!!"
543 WRITE (unit_nr, '(T3,A2,T8,A65,T79,A2)') '!!', "THERE ARE UNCONVERGED ENERGIES PRINTED OUT, SOMETHING WENT WRONG!", "!!"
544 WRITE (unit_nr, '(T3,A2,T79,A2)') '!!', "!!"
545 WRITE (unit_nr, '(T3,A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
546 END IF
547 END IF
548
549 CALL timestop(handle)
550
551 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 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 CALL timeset(routinen, handle)
594
595 WRITE (unit_nr, '(T3,A)') '******************************************************************************'
596 WRITE (unit_nr, '(T3,A2,T15,A49,T79,A2)') &
597 '**', "Parameters of the BSE-Davidson solver:", "**"
598 WRITE (unit_nr, '(T3,A2,T79,A2)') &
599 '**', "**"
600 WRITE (unit_nr, '(T3,A2,T79,A2)') &
601 '**', "**"
602 WRITE (unit_nr, '(T3,A2,T10,A16,I5,A12,I5,A8,T79,A2)') &
603 '**', "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 '**', "because of ", success_abort_string, " using ", &
606 bse_davidson_abort_cond_string, " threshold condition", "**"
607 WRITE (unit_nr, '(T3,A2,T79,A2)') &
608 '**', "**"
609 WRITE (unit_nr, '(T3,A2,T10,A32,T65,I11,T79,A2)') &
610 '**', "The Z space has at the end dim. ", num_z_vectors, "**"
611 WRITE (unit_nr, '(T3,A2,T10,A45,T65,I11,T79,A2)') &
612 '**', "Number of unconverged residuals in subspace: ", num_res_unconverged, "**"
613 WRITE (unit_nr, '(T3,A2,T10,A35,T65,E11.4,T79,A2)') &
614 '**', "largest unconverged residual (eV): ", max_res_norm*evolt, "**"
615 WRITE (unit_nr, '(T3,A2,T10,A45,T65,E11.4,T79,A2)') &
616 '**', "threshold for convergence of residuals (eV): ", eps_res*evolt, "**"
617 WRITE (unit_nr, '(T3,A2,T10,A45,T65,I11,T79,A2)') &
618 '**', "Number of desired, but unconverged energies: ", num_en_unconverged, "**"
619 WRITE (unit_nr, '(T3,A2,T10,A44,T65,E11.4,T79,A2)') &
620 '**', "largest unconverged energy difference (eV): ", max_en_diff*evolt, "**"
621 WRITE (unit_nr, '(T3,A2,T10,A44,T65,E11.4,T79,A2)') &
622 '**', "threshold for convergence of energies (eV): ", eps_exc_en*evolt, "**"
623 WRITE (unit_nr, '(T3,A2,T10,A40,T65,I11,T79,A2)') &
624 '**', "number of computed excitation energies: ", num_exc_en, "**"
625
626 IF (z_space_energy_cutoff > 0) THEN
627 WRITE (unit_nr, '(T3,A2,T10,A37,T65,E11.4,T79,A2)') &
628 '**', "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 '**', "number of Z space at the beginning: ", num_z_vectors_init, "**"
633 WRITE (unit_nr, '(T3,A2,T10,A30,T65,I11,T79,A2)') &
634 '**', "maximal dimension of Z space: ", num_max_z_space, "**"
635 WRITE (unit_nr, '(T3,A2,T10,A31,T65,I11,T79,A2)') &
636 '**', "added directions per iteration: ", num_new_t, "**"
637 WRITE (unit_nr, '(T3,A2,T79,A2)') &
638 '**', "**"
639 WRITE (unit_nr, '(T3,A2,T79,A2)') &
640 '**', "**"
641 WRITE (unit_nr, '(T3,A)') '******************************************************************************'
642 WRITE (unit_nr, '(T3,A)') ' '
643
644 CALL timestop(handle)
645
646 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 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 CALL timeset(routinen, handle)
673
674 num_en_unconverged = 0
675 ALLOCATE (en_diffs(num_exc_en))
676 DO mu_l = 1, num_exc_en
677 en_diffs(mu_l) = abs(subspace_full_eigenval(mu_l) - subspace_prev_eigenval(mu_l))
678 IF (en_diffs(mu_l) > eps_exc_en) num_en_unconverged = num_en_unconverged + 1
679 END DO
680 max_en_diff = maxval(en_diffs)
681
682 CALL timestop(handle)
683
684 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 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 CALL timeset(routinen, handle)
713
714 num_res_unconverged = 0
715 ALLOCATE (res_norms(num_new_t))
716 DO mu_l = 1, num_new_t
717 res_norms(mu_l) = norm2(subspace_residuals_reshaped(:, mu_l))
718 IF (res_norms(mu_l) > eps_res) THEN
719 num_res_unconverged = num_res_unconverged + 1
720 IF (unit_nr > 0) THEN
721 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 max_res_norm = maxval(res_norms)
726 IF (unit_nr > 0) THEN
727 WRITE (unit_nr, *) "Maximal unconverged res (of ", num_res_unconverged, &
728 " unconverged res in this step) in i_iter=", i_iter, "is:", max_res_norm
729 END IF
730
731 CALL timestop(handle)
732
733 END SUBROUTINE
734
735! **************************************************************************************************
736!> \brief ...
737!> \param W_vectors ...
738!> \param num_Z_vectors ...
739!> \param homo ...
740!> \param virtual ...
741! **************************************************************************************************
742 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tau_w, work_w, work_w_dor
752 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: w_vectors_reshaped
753
754 CALL timeset(routinen, handle)
755
756 ALLOCATE (w_vectors_reshaped(homo*virtual, num_z_vectors))
757 w_vectors_reshaped(:, :) = reshape(w_vectors, (/homo*virtual, num_z_vectors/))
758
759 ALLOCATE (tau_w(min(homo*virtual, num_z_vectors)))
760 tau_w = 0.0_dp
761
762 ALLOCATE (work_w(1))
763 work_w = 0.0_dp
764
765 ALLOCATE (work_w_dor(1))
766 work_w_dor = 0.0_dp
767
768 CALL dgeqrf(homo*virtual, num_z_vectors, w_vectors_reshaped, homo*virtual, tau_w, work_w, -1, info_orth)
769 lwork_w = int(work_w(1))
770 DEALLOCATE (work_w)
771 ALLOCATE (work_w(lwork_w))
772 work_w = 0.0_dp
773 CALL dgeqrf(homo*virtual, num_z_vectors, w_vectors_reshaped, homo*virtual, tau_w, work_w, lwork_w, info_orth)
774 IF (info_orth /= 0) THEN
775 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 tau_w, work_w_dor, -1, info_dor)
779 lwork_dor = int(work_w_dor(1))
780 DEALLOCATE (work_w_dor)
781 ALLOCATE (work_w_dor(lwork_dor))
782 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 tau_w, work_w_dor, lwork_dor, info_dor)
785 IF (info_orth /= 0) THEN
786 cpabort("QR Decomp Step 2 doesnt work")
787 END IF
788
789 w_vectors(:, :, :) = reshape(w_vectors_reshaped, (/homo, virtual, num_z_vectors/))
790
791 DEALLOCATE (work_w, work_w_dor, tau_w, w_vectors_reshaped)
792
793 CALL timestop(handle)
794
795 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 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 CALL timeset(routinen, handle)
825
826 ALLOCATE (subspace_add_dir_reshaped(homo*virtual, num_new_t))
827
828 prec_neg = 0
829 DO mu_subspace = 1, num_new_t
830 DO i_occ = 1, homo
831 DO a_virt = 1, virtual
832 !MG to check: Indexorder and range of indices
833 prec_scalar = -1/(subspace_new_eigenval(mu_subspace) - (eigenval(a_virt + homo) - eigenval(i_occ)))
834 IF (prec_scalar < 0) THEN
835 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 subspace_residuals_reshaped((i_occ - 1)*virtual + a_virt, mu_subspace)
840 END DO
841 END DO
842 END DO
843
844 subspace_add_dir(:, :, :) = reshape(subspace_add_dir_reshaped, (/homo, virtual, num_new_t/))
845
846 DEALLOCATE (subspace_add_dir_reshaped)
847 CALL timestop(handle)
848
849 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: subspace_res_a, subspace_res_ev
878
879 CALL timeset(routinen, handle)
880
881 ALLOCATE (subspace_res_ev(homo*virtual, num_new_t))
882 subspace_res_ev = 0.0_dp
883
884 ALLOCATE (subspace_res_a(homo*virtual, num_new_t))
885 subspace_res_a = 0.0_dp
886
887 !Compute all residuals in one loop, iterating over number of new/added t per iteration
888 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 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 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 - subspace_res_a(:, mu_subspace)
898
899 END DO
900 subspace_ritzvec(:, :, :) = reshape(subspace_res_ev, (/homo, virtual, num_new_t/))
901 DEALLOCATE (subspace_res_ev, subspace_res_a)
902
903 CALL timestop(handle)
904
905 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
936 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: zaz
937
938 CALL timeset(routinen, handle)
939
940 ALLOCATE (zaz(num_z_vectors, num_z_vectors))
941 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 DO i_z_vector = 1, num_z_vectors
946 DO j_z_vector = 1, num_z_vectors
947 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 IF (unit_nr > 0) THEN
951 WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "Before Diag"
952 END IF
953
954 !MG to do: Check for symmetry of ZAZ!
955 ALLOCATE (work(1))
956 work = 0.0_dp
957 CALL dsyev("V", "U", num_z_vectors, zaz, num_z_vectors, subspace_full_eigenval, work, -1, zaz_diag_info)
958 lwork = int(work(1))
959 DEALLOCATE (work)
960 ALLOCATE (work(lwork))
961 work = 0.0_dp
962 !MG to check: Usage of symmetric routine okay? (Correct LWORK?)
963 CALL dsyev("V", "U", num_z_vectors, zaz, num_z_vectors, subspace_full_eigenval, work, lwork, zaz_diag_info)
964
965 IF (zaz_diag_info /= 0) THEN
966 cpabort("ZAZ could not be diagonalized successfully.")
967 END IF
968
969 IF (unit_nr > 0) THEN
970 WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "After Diag"
971 END IF
972
973 subspace_new_eigenval(1:num_new_t) = subspace_full_eigenval(1:num_new_t)
974 subspace_new_eigenvec(:, 1:num_new_t) = zaz(:, 1:num_new_t)
975 DEALLOCATE (work)
976 DEALLOCATE (zaz)
977
978 CALL timestop(handle)
979
980 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 SUBROUTINE compute_bz(BZ, Z_vectors, B_iaQ_bse_local, B_bar_iaQ_bse_local, &
996 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 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 DO i_z_vector = 1, num_z_vectors
1013
1014 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 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 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 CALL para_env%sum(bz)
1029
1030 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 SUBROUTINE compute_az(AZ, Z_vectors, B_iaQ_bse_local, B_bar_ijQ_bse_local, B_abQ_bse_local, M_ia_tmp, &
1055 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
1079 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a_full_reshaped
1080 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: a_full, v_iajb, w_ijab
1081
1082 CALL timeset(routinen, handle)
1083 az(:, :, :) = 0.0_dp
1084
1085 IF (i_iter == 1 .AND. bse_full_diag_debug) THEN
1086 ALLOCATE (w_ijab(homo, homo, virtual, virtual))
1087 ALLOCATE (a_full(homo, virtual, homo, virtual))
1088 ALLOCATE (a_full_reshaped(homo*virtual, homo*virtual))
1089 ALLOCATE (full_exc_spectrum(homo*virtual))
1090 w_ijab = 0.0_dp
1091 a_full = 0.0_dp
1092 a_full_reshaped = 0.0_dp
1093 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 para_env)
1099
1100 DO i_z_vector = 1, num_z_vectors
1101
1102 DO lll = 1, local_ri_size
1103
1104 ! M_ja^P = sum_b Z_jb*B_ba^P
1105 CALL dgemm("N", "N", homo, virtual, virtual, 1.0_dp, z_vectors(:, :, i_z_vector), homo, &
1106 b_abq_bse_local(:, :, lll), virtual, 0.0_dp, m_ia_tmp, homo)
1107
1108 ! (AZ)_ia = sum_jP B_bar_ij^P*M_ja^P
1109 CALL dgemm("N", "N", homo, virtual, homo, -1.0_dp, b_bar_ijq_bse_local(:, :, lll), homo, &
1110 m_ia_tmp, homo, 1.0_dp, az(:, :, i_z_vector), homo)
1111
1112 END DO
1113 END DO
1114
1115 IF (i_iter == 1 .AND. bse_full_diag_debug) THEN
1116 w_ijab = 0.0_dp
1117 !Create screened 4c integrals for check
1118 DO lll = 1, local_ri_size
1119 DO i = 1, homo
1120 DO j = 1, homo
1121 DO a = 1, virtual
1122 DO b = 1, virtual
1123 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)
1124 END DO
1125 END DO
1126 END DO
1127 END DO
1128 END DO
1129 ! we make the mp_sum to sum over all RI basis functions
1130 CALL para_env%sum(w_ijab)
1131 END IF
1132
1133 ! we make the mp_sum to sum over all RI basis functions
1134 CALL para_env%sum(az)
1135
1136 ! add (e_a-e_i)*Z_ia
1137 DO i_occ = 1, homo
1138 DO a_virt = 1, virtual
1139
1140 eigen_diff = eigenval(a_virt + homo) - eigenval(i_occ)
1141 IF (unit_nr > 0 .AND. i_iter == 1) THEN
1142 WRITE (unit_nr, *) "Ediff at (i_occ,a_virt)=", i_occ, a_virt, " is: ", eigen_diff
1143 END IF
1144
1145 az(i_occ, a_virt, :) = az(i_occ, a_virt, :) + z_vectors(i_occ, a_virt, :)*eigen_diff
1146
1147 END DO
1148 END DO
1149
1150 !cut off contributions, which are too high in the excitation spectrum
1151 IF (z_space_energy_cutoff > 0) THEN
1152 DO i_occ = 1, homo
1153 DO a_virt = 1, virtual
1154
1155 IF (eigenval(a_virt + homo) > z_space_energy_cutoff .OR. -eigenval(i_occ) > z_space_energy_cutoff) THEN
1156 az(i_occ, a_virt, :) = 0
1157 END IF
1158
1159 END DO
1160 END DO
1161 END IF
1162
1163 !Debugging purposes: full diagonalization of A
1164 IF (i_iter == 1 .AND. bse_full_diag_debug) THEN
1165 n = 0
1166 DO i = 1, homo
1167 DO a = 1, virtual
1168 n = n + 1
1169 m = 0
1170 DO j = 1, homo
1171 DO b = 1, virtual
1172 m = m + 1
1173 IF (a == b .AND. i == j) THEN
1174 eigen_diff = eigenval(a + homo) - eigenval(i)
1175 ELSE
1176 eigen_diff = 0
1177 END IF
1178 a_full_reshaped(n, m) = eigen_diff + 2*v_iajb(i, a, j, b) - w_ijab(i, j, a, b)
1179 a_full(i, a, j, b) = eigen_diff + 2*v_iajb(i, a, j, b) - w_ijab(i, j, a, b)
1180 END DO
1181 END DO
1182 END DO
1183 END DO
1184
1185 !MG to do: Check for symmetry of ZAZ!
1186 ALLOCATE (work(1))
1187 work = 0.0_dp
1188 CALL dsyev("N", "U", homo*virtual, a_full_reshaped, homo*virtual, full_exc_spectrum, work, -1, diag_info)
1189 lwork = int(work(1))
1190 DEALLOCATE (work)
1191 ALLOCATE (work(lwork))
1192 work = 0.0_dp
1193 !MG to check: Usage of symmetric routine okay? (Correct LWORK?)
1194 CALL dsyev("N", "U", homo*virtual, a_full_reshaped, homo*virtual, full_exc_spectrum, work, lwork, diag_info)
1195
1196 DEALLOCATE (work)
1197
1198 DEALLOCATE (w_ijab, v_iajb, a_full, a_full_reshaped)
1199 END IF
1200
1201 CALL timestop(handle)
1202
1203 END SUBROUTINE
1204
1205! **************************************************************************************************
1206!> \brief ...
1207!> \param AZ ...
1208!> \param Z_vectors ...
1209!> \param B_iaQ_bse_local ...
1210!> \param RI_vector ...
1211!> \param local_RI_size ...
1212!> \param num_Z_vectors ...
1213!> \param homo ...
1214!> \param virtual ...
1215!> \param bse_spin_config ...
1216!> \param v_iajb ...
1217!> \param bse_full_diag_debug ...
1218!> \param i_iter ...
1219!> \param para_env ...
1220! **************************************************************************************************
1221 SUBROUTINE compute_v_ia_jb_part(AZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, &
1222 num_Z_vectors, homo, virtual, bse_spin_config, v_iajb, bse_full_diag_debug, i_iter, &
1223 para_env)
1224
1225 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: az, z_vectors, b_iaq_bse_local
1226 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: ri_vector
1227 INTEGER, INTENT(IN) :: local_ri_size, num_z_vectors, homo, &
1228 virtual, bse_spin_config
1229 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: v_iajb
1230 LOGICAL :: bse_full_diag_debug
1231 INTEGER, INTENT(IN) :: i_iter
1232 TYPE(mp_para_env_type), INTENT(IN) :: para_env
1233
1234 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_v_ia_jb_part'
1235
1236 INTEGER :: a, a_virt, b, handle, i, i_occ, &
1237 i_z_vector, j, lll
1238 REAL(kind=dp) :: alpha
1239
1240!debugging:
1241
1242 CALL timeset(routinen, handle)
1243
1244 !Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
1245 SELECT CASE (bse_spin_config)
1246 CASE (bse_singlet)
1247 alpha = 2.0_dp
1248 CASE (bse_triplet)
1249 alpha = 0.0_dp
1250 END SELECT
1251
1252 ri_vector = 0.0_dp
1253
1254 ! v_P = sum_jb B_jb^P Z_jb
1255 DO lll = 1, local_ri_size
1256 DO i_z_vector = 1, num_z_vectors
1257 DO i_occ = 1, homo
1258 DO a_virt = 1, virtual
1259
1260 ri_vector(lll, i_z_vector) = ri_vector(lll, i_z_vector) + &
1261 z_vectors(i_occ, a_virt, i_z_vector)* &
1262 b_iaq_bse_local(i_occ, a_virt, lll)
1263
1264 END DO
1265 END DO
1266 END DO
1267 END DO
1268
1269 ! AZ = sum_P B_ia^P*v_P + ...
1270 DO lll = 1, local_ri_size
1271 DO i_z_vector = 1, num_z_vectors
1272 DO i_occ = 1, homo
1273 DO a_virt = 1, virtual
1274 !MG to check: Minus sign at v oder W? Factor for triplet/singlet
1275 az(i_occ, a_virt, i_z_vector) = az(i_occ, a_virt, i_z_vector) + &
1276 alpha*ri_vector(lll, i_z_vector)* &
1277 b_iaq_bse_local(i_occ, a_virt, lll)
1278
1279 END DO
1280 END DO
1281 END DO
1282 END DO
1283 IF (i_iter == 1 .AND. bse_full_diag_debug) THEN
1284 ALLOCATE (v_iajb(homo, virtual, homo, virtual))
1285 v_iajb = 0.0_dp
1286 !Create unscreened 4c integrals for check
1287 DO lll = 1, local_ri_size
1288 DO i = 1, homo
1289 DO j = 1, homo
1290 DO a = 1, virtual
1291 DO b = 1, virtual
1292 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)
1293 END DO
1294 END DO
1295 END DO
1296 END DO
1297 END DO
1298 ! we make the mp_sum to sum over all RI basis functions
1299 CALL para_env%sum(v_iajb)
1300 END IF
1301
1302 CALL timestop(handle)
1303
1304 END SUBROUTINE
1305
1306! **************************************************************************************************
1307!> \brief ...Eigenval
1308!> \param Z_vectors ...
1309!> \param Eigenval ...
1310!> \param num_Z_vectors ...
1311!> \param homo ...
1312!> \param virtual ...
1313! **************************************************************************************************
1314 SUBROUTINE initial_guess_z_vectors(Z_vectors, Eigenval, num_Z_vectors, homo, virtual)
1315
1316 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: z_vectors
1317 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenval
1318 INTEGER, INTENT(IN) :: num_z_vectors, homo, virtual
1319
1320 CHARACTER(LEN=*), PARAMETER :: routinen = 'initial_guess_Z_vectors'
1321
1322 INTEGER :: a_virt, handle, i_occ, i_z_vector, &
1323 min_loc(2)
1324 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigen_diff_ia
1325
1326 CALL timeset(routinen, handle)
1327
1328 ALLOCATE (eigen_diff_ia(homo, virtual))
1329
1330 DO i_occ = 1, homo
1331 DO a_virt = 1, virtual
1332 eigen_diff_ia(i_occ, a_virt) = eigenval(a_virt + homo) - eigenval(i_occ)
1333 END DO
1334 END DO
1335
1336 DO i_z_vector = 1, num_z_vectors
1337
1338 min_loc = minloc(eigen_diff_ia)
1339
1340 z_vectors(min_loc(1), min_loc(2), i_z_vector) = 1.0_dp
1341
1342 eigen_diff_ia(min_loc(1), min_loc(2)) = 1.0e20_dp
1343
1344 END DO
1345
1346 DEALLOCATE (eigen_diff_ia)
1347
1348 CALL timestop(handle)
1349
1350 END SUBROUTINE
1351
1352 ! **************************************************************************************************
1353!> \brief ...
1354!> \param fm_mat_S_ab_bse ...
1355!> \param fm_mat_S ...
1356!> \param fm_mat_S_bar_ia_bse ...
1357!> \param fm_mat_S_bar_ij_bse ...
1358!> \param B_bar_ijQ_bse_local ...
1359!> \param B_abQ_bse_local ...
1360!> \param B_bar_iaQ_bse_local ...
1361!> \param B_iaQ_bse_local ...
1362!> \param dimen_RI ...
1363!> \param homo ...
1364!> \param virtual ...
1365!> \param gd_array ...
1366!> \param color_sub ...
1367!> \param para_env ...
1368! **************************************************************************************************
1369 SUBROUTINE fill_local_3c_arrays(fm_mat_S_ab_bse, fm_mat_S, &
1370 fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
1371 B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
1372 B_iaQ_bse_local, dimen_RI, homo, virtual, &
1373 gd_array, color_sub, para_env)
1374
1375 TYPE(cp_fm_type), INTENT(IN) :: fm_mat_s_ab_bse, fm_mat_s, &
1376 fm_mat_s_bar_ia_bse, &
1377 fm_mat_s_bar_ij_bse
1378 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1379 INTENT(OUT) :: b_bar_ijq_bse_local, b_abq_bse_local, &
1380 b_bar_iaq_bse_local, b_iaq_bse_local
1381 INTEGER, INTENT(IN) :: dimen_ri, homo, virtual
1382 TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
1383 INTEGER, INTENT(IN) :: color_sub
1384 TYPE(mp_para_env_type), INTENT(IN) :: para_env
1385
1386 CHARACTER(LEN=*), PARAMETER :: routinen = 'fill_local_3c_arrays'
1387
1388 INTEGER :: handle
1389
1390 CALL timeset(routinen, handle)
1391
1392 CALL allocate_and_fill_local_array(b_iaq_bse_local, fm_mat_s, gd_array, color_sub, homo, virtual, dimen_ri, para_env)
1393
1394 CALL allocate_and_fill_local_array(b_bar_iaq_bse_local, fm_mat_s_bar_ia_bse, gd_array, color_sub, homo, virtual, &
1395 dimen_ri, para_env)
1396
1397 CALL allocate_and_fill_local_array(b_bar_ijq_bse_local, fm_mat_s_bar_ij_bse, gd_array, color_sub, homo, homo, &
1398 dimen_ri, para_env)
1399
1400 CALL allocate_and_fill_local_array(b_abq_bse_local, fm_mat_s_ab_bse, gd_array, color_sub, virtual, virtual, &
1401 dimen_ri, para_env)
1402
1403 CALL timestop(handle)
1404
1405 END SUBROUTINE
1406
1407! **************************************************************************************************
1408!> \brief ...
1409!> \param B_local ...
1410!> \param fm_mat_S ...
1411!> \param gd_array ...
1412!> \param color_sub ...
1413!> \param small_size ...
1414!> \param big_size ...
1415!> \param dimen_RI ...
1416!> \param para_env ...
1417! **************************************************************************************************
1418 SUBROUTINE allocate_and_fill_local_array(B_local, fm_mat_S, gd_array, &
1419 color_sub, small_size, big_size, dimen_RI, para_env)
1420
1421 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1422 INTENT(OUT) :: b_local
1423 TYPE(cp_fm_type), INTENT(IN) :: fm_mat_s
1424 TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
1425 INTEGER, INTENT(IN) :: color_sub, small_size, big_size, dimen_ri
1426 TYPE(mp_para_env_type), INTENT(IN) :: para_env
1427
1428 CHARACTER(LEN=*), PARAMETER :: routinen = 'allocate_and_fill_local_array'
1429
1430 INTEGER :: combi_index, end_ri, handle, handle1, i_comm, i_entry, iib, imepos, jjb, &
1431 level_big_size, level_small_size, ncol_local, nrow_local, num_comm_cycles, ri_index, &
1432 size_ri, start_ri
1433 INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, mepos_from_ri_index, &
1434 num_entries_rec, num_entries_send
1435 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1436 REAL(kind=dp) :: matrix_el
1437 TYPE(integ_mat_buffer_type), ALLOCATABLE, &
1438 DIMENSION(:) :: buffer_rec, buffer_send
1439 TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
1440
1441 CALL timeset(routinen, handle)
1442
1443 ALLOCATE (mepos_from_ri_index(dimen_ri))
1444 mepos_from_ri_index = 0
1445
1446 DO imepos = 0, para_env%num_pe - 1
1447
1448 CALL get_group_dist(gd_array, pos=imepos, starts=start_ri, ends=end_ri)
1449
1450 mepos_from_ri_index(start_ri:end_ri) = imepos
1451
1452 END DO
1453
1454 ! color_sub is automatically the number of the process since every subgroup has only one MPI rank
1455 CALL get_group_dist(gd_array, color_sub, start_ri, end_ri, size_ri)
1456
1457 ALLOCATE (b_local(small_size, big_size, 1:size_ri))
1458
1459 ALLOCATE (num_entries_send(0:para_env%num_pe - 1))
1460 ALLOCATE (num_entries_rec(0:para_env%num_pe - 1))
1461
1462 ALLOCATE (req_array(1:para_env%num_pe, 4))
1463
1464 ALLOCATE (entry_counter(0:para_env%num_pe - 1))
1465
1466 CALL cp_fm_get_info(matrix=fm_mat_s, &
1467 nrow_local=nrow_local, &
1468 ncol_local=ncol_local, &
1469 row_indices=row_indices, &
1470 col_indices=col_indices)
1471
1472 num_comm_cycles = 10
1473
1474 ! communicate not all due to huge memory overhead, since for every number in fm_mat_S, we store
1475 ! three additional ones (RI index, first MO index, second MO index!!)
1476 DO i_comm = 0, num_comm_cycles - 1
1477
1478 num_entries_send = 0
1479 num_entries_rec = 0
1480
1481 ! loop over RI index to get the number of sent entries
1482 DO jjb = 1, nrow_local
1483
1484 ri_index = row_indices(jjb)
1485
1486 IF (modulo(ri_index, num_comm_cycles) /= i_comm) cycle
1487
1488 imepos = mepos_from_ri_index(ri_index)
1489
1490 num_entries_send(imepos) = num_entries_send(imepos) + ncol_local
1491
1492 END DO
1493
1494 CALL para_env%alltoall(num_entries_send, num_entries_rec, 1)
1495
1496 ALLOCATE (buffer_rec(0:para_env%num_pe - 1))
1497 ALLOCATE (buffer_send(0:para_env%num_pe - 1))
1498
1499 ! allocate data message and corresponding indices
1500 DO imepos = 0, para_env%num_pe - 1
1501
1502 ALLOCATE (buffer_rec(imepos)%msg(num_entries_rec(imepos)))
1503 buffer_rec(imepos)%msg = 0.0_dp
1504
1505 ALLOCATE (buffer_send(imepos)%msg(num_entries_send(imepos)))
1506 buffer_send(imepos)%msg = 0.0_dp
1507
1508 ALLOCATE (buffer_rec(imepos)%indx(num_entries_rec(imepos), 3))
1509 buffer_rec(imepos)%indx = 0
1510
1511 ALLOCATE (buffer_send(imepos)%indx(num_entries_send(imepos), 3))
1512 buffer_send(imepos)%indx = 0
1513
1514 END DO
1515
1516 entry_counter(:) = 0
1517
1518 ! loop over RI index for filling the send-buffer
1519 DO jjb = 1, nrow_local
1520
1521 ri_index = row_indices(jjb)
1522
1523 IF (modulo(ri_index, num_comm_cycles) /= i_comm) cycle
1524
1525 imepos = mepos_from_ri_index(ri_index)
1526
1527 DO iib = 1, ncol_local
1528
1529 combi_index = col_indices(iib)
1530 level_small_size = max(1, combi_index - 1)/max(big_size, 2) + 1
1531 level_big_size = combi_index - (level_small_size - 1)*big_size
1532
1533 entry_counter(imepos) = entry_counter(imepos) + 1
1534
1535 buffer_send(imepos)%msg(entry_counter(imepos)) = fm_mat_s%local_data(jjb, iib)
1536
1537 buffer_send(imepos)%indx(entry_counter(imepos), 1) = ri_index
1538 buffer_send(imepos)%indx(entry_counter(imepos), 2) = level_small_size
1539 buffer_send(imepos)%indx(entry_counter(imepos), 3) = level_big_size
1540
1541 END DO
1542
1543 END DO
1544
1545 CALL timeset("BSE_comm_data", handle1)
1546
1547 CALL communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, buffer_send, req_array)
1548
1549 CALL timestop(handle1)
1550
1551 ! fill B_local
1552 DO imepos = 0, para_env%num_pe - 1
1553
1554 DO i_entry = 1, num_entries_rec(imepos)
1555
1556 ri_index = buffer_rec(imepos)%indx(i_entry, 1) - start_ri + 1
1557 level_small_size = buffer_rec(imepos)%indx(i_entry, 2)
1558 level_big_size = buffer_rec(imepos)%indx(i_entry, 3)
1559
1560 matrix_el = buffer_rec(imepos)%msg(i_entry)
1561
1562 b_local(level_small_size, level_big_size, ri_index) = matrix_el
1563
1564 END DO
1565
1566 END DO
1567
1568 DO imepos = 0, para_env%num_pe - 1
1569 DEALLOCATE (buffer_send(imepos)%msg)
1570 DEALLOCATE (buffer_send(imepos)%indx)
1571 DEALLOCATE (buffer_rec(imepos)%msg)
1572 DEALLOCATE (buffer_rec(imepos)%indx)
1573 END DO
1574
1575 DEALLOCATE (buffer_rec, buffer_send)
1576
1577 END DO
1578
1579 DEALLOCATE (num_entries_send, num_entries_rec)
1580
1581 DEALLOCATE (mepos_from_ri_index)
1582
1583 DEALLOCATE (entry_counter, req_array)
1584
1585 CALL timestop(handle)
1586
1587 END SUBROUTINE
1588
1589END MODULE bse_iterative
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:206
group_dist_types::get_group_dist
Definition group_dist_types.F:47
bse_iterative
Iterative routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_iterative.F:14
bse_iterative::do_subspace_iterations
subroutine, public do_subspace_iterations(b_bar_ijq_bse_local, b_abq_bse_local, b_bar_iaq_bse_local, b_iaq_bse_local, homo, virtual, bse_spin_config, unit_nr, eigenval, para_env, mp2_env)
...
Definition bse_iterative.F:57
bse_iterative::fill_local_3c_arrays
subroutine, public fill_local_3c_arrays(fm_mat_s_ab_bse, fm_mat_s, fm_mat_s_bar_ia_bse, fm_mat_s_bar_ij_bse, b_bar_ijq_bse_local, b_abq_bse_local, b_bar_iaq_bse_local, b_iaq_bse_local, dimen_ri, homo, virtual, gd_array, color_sub, para_env)
...
Definition bse_iterative.F:1374
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1009
group_dist_types
Types to describe group distributions.
Definition group_dist_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::bse_singlet
integer, parameter, public bse_singlet
Definition input_constants.F:1311
input_constants::bse_triplet
integer, parameter, public bse_triplet
Definition input_constants.F:1311
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mp2_types
Types needed for MP2 calculations.
Definition mp2_types.F:14
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
rpa_communication
Auxiliary routines necessary to redistribute an fm_matrix from a given blacs_env to another.
Definition rpa_communication.F:14
rpa_communication::communicate_buffer
subroutine, public communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, buffer_send, req_array, do_indx, do_msg)
...
Definition rpa_communication.F:434
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
group_dist_types::group_dist_d1_type
Definition group_dist_types.F:35
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
message_passing::mp_request_type
Definition message_passing.F:623
mp2_types::integ_mat_buffer_type
Definition mp2_types.F:409
mp2_types::mp2_type
Definition mp2_types.F:345