(git:5e7fe52)
Loading...
Searching...
No Matches
bse_full_diag.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines for the full diagonalization of GW + Bethe-Salpeter for computing
10!> electronic excitations
11!> \par History
12!> 10.2023 created [Maximilian Graml]
13! **************************************************************************************************
15
42 USE cp_fm_types, ONLY: cp_fm_create,&
55 USE kinds, ONLY: dp
57 USE mp2_types, ONLY: mp2_type
61#include "./base/base_uses.f90"
62
63 IMPLICIT NONE
64
65 PRIVATE
66
67 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_full_diag'
68
71
72CONTAINS
73
74! **************************************************************************************************
75!> \brief Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W)
76!> A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
77!> ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
78!> α is a spin-dependent factor
79!> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
80!> W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
81!> \param fm_mat_S_ia_bse ...
82!> \param fm_mat_S_bar_ij_bse ...
83!> \param fm_mat_S_ab_bse ...
84!> \param fm_A ...
85!> \param Eigenval ...
86!> \param unit_nr ...
87!> \param homo ...
88!> \param virtual ...
89!> \param dimen_RI ...
90!> \param mp2_env ...
91!> \param para_env ...
92!> \param qs_env ...
93! **************************************************************************************************
94 SUBROUTINE create_a(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, &
95 fm_A, Eigenval, unit_nr, &
96 homo, virtual, dimen_RI, mp2_env, &
97 para_env, qs_env)
98
99 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_s_ia_bse, fm_mat_s_bar_ij_bse, &
100 fm_mat_s_ab_bse
101 TYPE(cp_fm_type), INTENT(INOUT) :: fm_a
102 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenval
103 INTEGER, INTENT(IN) :: unit_nr
104 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
105 INTEGER, INTENT(IN) :: dimen_ri
106 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
107 TYPE(mp_para_env_type), INTENT(INOUT) :: para_env
108 TYPE(qs_environment_type), POINTER :: qs_env
109
110 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_A'
111
112 INTEGER :: a_virt_row, handle, i_occ_row, i_row_global, ii, isp, j_col_global, jj, k_isp, &
113 k_ov, n_ov_joint, ncol_local_a, nrow_local_a, nspins, sizeeigen
114 INTEGER, ALLOCATABLE, DIMENSION(:) :: eig_offsets, n_ov, offsets
115 INTEGER, DIMENSION(4) :: reordering
116 INTEGER, DIMENSION(:), POINTER :: col_indices_a, row_indices_a
117 REAL(kind=dp) :: alpha, alpha_screening, eigen_diff
118 TYPE(cp_blacs_env_type), POINTER :: blacs_env
119 TYPE(cp_fm_struct_type), POINTER :: fm_struct_a, fm_struct_s_joint, &
120 fm_struct_w
121 TYPE(cp_fm_type) :: fm_a_copy, fm_s_joint, fm_w
122 TYPE(dft_control_type), POINTER :: dft_control
123 TYPE(excited_energy_type), POINTER :: ex_env
124 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
125
126 CALL timeset(routinen, handle)
127
128 nspins = SIZE(homo)
129 ALLOCATE (n_ov(nspins), offsets(nspins), eig_offsets(nspins))
130 CALL get_bse_spin_block_layout(homo, virtual, n_ov, offsets, n_ov_joint)
131 ! Flat Eigenval layout: sigma-block isp at eig_offsets(isp)+1 .. eig_offsets(isp)+homo(isp)+virtual(isp)
132 eig_offsets(1) = 0
133 DO isp = 2, nspins
134 eig_offsets(isp) = eig_offsets(isp - 1) + homo(isp - 1) + virtual(isp - 1)
135 END DO
136
137 NULLIFY (dft_control, tddfpt_control)
138 CALL get_qs_env(qs_env, dft_control=dft_control)
139 tddfpt_control => dft_control%tddfpt2_control
140
141 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
142 WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating A'
143 END IF
144
145 !Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
146 SELECT CASE (mp2_env%bse%bse_spin_config)
147 CASE (bse_singlet)
148 alpha = 2.0_dp
149 CASE (bse_triplet)
150 alpha = 0.0_dp
151 END SELECT
152 ! For open-shell (nspins>1): each spin block contributes once; SPIN_CONFIG is ignored.
153 IF (nspins > 1) THEN
154 CALL cp_warn(__location__, &
155 "BSE: SPIN_CONFIG ignored for open-shell reference; using alpha=1.")
156 alpha = 1.0_dp
157 END IF
158
159 IF (mp2_env%bse%screening_method == bse_screening_alpha) THEN
160 alpha_screening = mp2_env%bse%screening_factor
161 ELSE
162 alpha_screening = 1.0_dp
163 END IF
164
165 ! create the blacs env for ij matrices (NOT fm_mat_S_ia_bse%matrix_struct related parallel_gemms!)
166 NULLIFY (blacs_env)
167 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
168
169 ! We have to use the same blacs_env for A as for the matrices fm_mat_S_ia_bse from RPA
170 ! Logic: A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
171 ! We create v_ia,jb and W_ij,ab, then we communicate entries from local W_ij,ab
172 ! to the full matrix v_ia,jb. By adding these and the energy diffenences: v_ia,jb -> A_ia,jb
173 ! We use the A matrix already from the start instead of v
174 CALL cp_fm_struct_create(fm_struct_a, context=fm_mat_s_ia_bse(1)%matrix_struct%context, &
175 nrow_global=n_ov_joint, ncol_global=n_ov_joint, &
176 para_env=fm_mat_s_ia_bse(1)%matrix_struct%para_env)
177 CALL cp_fm_create(fm_a, fm_struct_a, name="fm_A_iajb")
178 CALL cp_fm_set_all(fm_a, 0.0_dp)
179 ! fm_A_copy only used in the TDDFPT do_bse_w_only path (closed-shell only)
180 IF (tddfpt_control%do_bse_w_only .AND. nspins == 1) THEN
181 CALL cp_fm_create(fm_a_copy, fm_struct_a, name="fm_A_iajb")
182 CALL cp_fm_set_all(fm_a_copy, 0.0_dp)
183 END IF
184
185 ! Create A matrix from GW Energies, v_ia,jb and W_ij,ab
186 ! v_ia,jb = \sum_P B^P_ia B^P_jb (Coulomb)
187 IF ((.NOT. tddfpt_control%do_bse) .AND. (.NOT. tddfpt_control%do_bse_w_only)) THEN
188 IF (nspins > 1) THEN
189 ! Assemble joint ia-slab for a single Coulomb gemm across all spin blocks
190 CALL cp_fm_struct_create(fm_struct_s_joint, &
191 context=fm_mat_s_ia_bse(1)%matrix_struct%context, &
192 nrow_global=dimen_ri, ncol_global=n_ov_joint, &
193 para_env=fm_mat_s_ia_bse(1)%matrix_struct%para_env)
194 CALL cp_fm_create(fm_s_joint, fm_struct_s_joint, name="fm_S_ia_joint")
195 CALL cp_fm_set_all(fm_s_joint, 0.0_dp)
196 CALL assemble_joint_ov_slab(fm_mat_s_ia_bse, offsets, n_ov, dimen_ri, fm_s_joint)
197 CALL parallel_gemm(transa="T", transb="N", m=n_ov_joint, n=n_ov_joint, k=dimen_ri, &
198 alpha=alpha, matrix_a=fm_s_joint, matrix_b=fm_s_joint, beta=0.0_dp, &
199 matrix_c=fm_a)
200 CALL cp_fm_release(fm_s_joint)
201 CALL cp_fm_struct_release(fm_struct_s_joint)
202 ELSE
203 CALL parallel_gemm(transa="T", transb="N", m=homo(1)*virtual(1), n=homo(1)*virtual(1), &
204 k=dimen_ri, alpha=alpha, &
205 matrix_a=fm_mat_s_ia_bse(1), matrix_b=fm_mat_s_ia_bse(1), &
206 beta=0.0_dp, matrix_c=fm_a)
207 END IF
208 END IF
209
210 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
211 WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated A_iajb'
212 END IF
213
214 ! W term on sigma-diagonal blocks only: W^sigma_ij,ab = sum_P barB^P_ij B^P_ab
215 ! offsets(isp) places each block at the correct position in joint A.
216 ! For nspins=1: offsets(1)=0, equivalent to the original code.
217 DO isp = 1, nspins
218 IF (mp2_env%bse%screening_method /= bse_screening_rpa) THEN
219 CALL cp_fm_struct_create(fm_struct_w, context=fm_mat_s_ab_bse(isp)%matrix_struct%context, &
220 nrow_global=homo(isp)**2, ncol_global=virtual(isp)**2, &
221 para_env=fm_mat_s_ab_bse(isp)%matrix_struct%para_env)
222 CALL cp_fm_create(fm_w, fm_struct_w, name="fm_W_ijab")
223 CALL cp_fm_set_all(fm_w, 0.0_dp)
224 !W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab
225 CALL parallel_gemm(transa="T", transb="N", m=homo(isp)**2, n=virtual(isp)**2, &
226 k=dimen_ri, alpha=alpha_screening, &
227 matrix_a=fm_mat_s_bar_ij_bse(isp), matrix_b=fm_mat_s_ab_bse(isp), &
228 beta=0.0_dp, matrix_c=fm_w)
229 reordering = [1, 3, 2, 4]
230 CALL fm_general_add_bse(fm_a, fm_w, -1.0_dp, homo(isp), virtual(isp), &
231 virtual(isp), virtual(isp), unit_nr, reordering, mp2_env, &
232 row_offset=offsets(isp), col_offset=offsets(isp))
233 IF (nspins == 1 .AND. tddfpt_control%do_bse_w_only) THEN
234 CALL fm_general_add_bse(fm_a_copy, fm_w, -1.0_dp, homo(1), virtual(1), &
235 virtual(1), virtual(1), unit_nr, reordering, mp2_env)
236 END IF
237 ! W and A stash for TDDFPT path (closed-shell only; open-shell deferred)
238 IF (nspins == 1) THEN
239 IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
240 tddfpt_control%do_bse_gw_only) THEN
241 NULLIFY (ex_env)
242 CALL get_qs_env(qs_env, exstate_env=ex_env)
243 IF (.NOT. tddfpt_control%do_bse_gw_only) THEN
244 ALLOCATE (ex_env%bse_w_matrix_MO(1, 1))
245 ALLOCATE (ex_env%bse_a_matrix_MO(1, 1))
246 CALL cp_fm_create(ex_env%bse_w_matrix_MO(1, 1), fm_struct_w)
247 CALL cp_fm_create(ex_env%bse_a_matrix_MO(1, 1), fm_struct_a)
248 CALL cp_fm_to_fm(fm_w, ex_env%bse_w_matrix_MO(1, 1))
249 IF (tddfpt_control%do_bse_w_only) THEN
250 CALL cp_fm_to_fm(fm_a_copy, ex_env%bse_a_matrix_MO(1, 1))
251 ELSE
252 CALL cp_fm_to_fm(fm_a, ex_env%bse_a_matrix_MO(1, 1))
253 END IF
254 END IF
255 END IF
256 END IF
257 CALL cp_fm_release(fm_w)
258 CALL cp_fm_struct_release(fm_struct_w)
259 END IF
260 END DO
261
262 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
263 WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated W_ijab'
264 END IF
265 IF (nspins == 1 .AND. tddfpt_control%do_bse_w_only) CALL cp_fm_release(fm_a_copy)
266
267 ! Get local row/col indices for direct diagonal access
268 CALL cp_fm_get_info(matrix=fm_a, nrow_local=nrow_local_a, ncol_local=ncol_local_a, &
269 row_indices=row_indices_a, col_indices=col_indices_a)
270
271 !Add (ε_a-ε_i) on the diagonal of each sigma-block; cross-spin blocks have no ε contribution.
272 IF (.NOT. tddfpt_control%do_bse) THEN
273 DO ii = 1, nrow_local_a
274 i_row_global = row_indices_a(ii)
275 DO jj = 1, ncol_local_a
276 j_col_global = col_indices_a(jj)
277 IF (i_row_global == j_col_global) THEN
278 ! Decode spin: isp such that i_row_global in [offsets(isp)+1, offsets(isp)+n_ov(isp)]
279 isp = nspins
280 DO k_isp = 1, nspins - 1
281 IF (i_row_global <= offsets(k_isp) + n_ov(k_isp)) THEN
282 isp = k_isp
283 EXIT
284 END IF
285 END DO
286 k_ov = i_row_global - offsets(isp)
287 i_occ_row = (k_ov - 1)/virtual(isp) + 1
288 a_virt_row = mod(k_ov - 1, virtual(isp)) + 1
289 eigen_diff = eigenval(eig_offsets(isp) + a_virt_row + homo(isp)) - &
290 eigenval(eig_offsets(isp) + i_occ_row)
291 fm_a%local_data(ii, jj) = fm_a%local_data(ii, jj) + eigen_diff
292 END IF
293 END DO
294 END DO
295 END IF
296
297 ! GW eigenvalue stash for TDDFPT path (closed-shell only)
298 IF (nspins == 1) THEN
299 IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
300 tddfpt_control%do_bse_gw_only) THEN
301 sizeeigen = SIZE(eigenval)
302 ALLOCATE (ex_env%gw_eigen(sizeeigen))
303 ex_env%gw_eigen(:) = eigenval(:)
304 END IF
305 END IF
306
307 CALL cp_fm_struct_release(fm_struct_a)
308 DEALLOCATE (n_ov, offsets, eig_offsets)
309
310 CALL cp_blacs_env_release(blacs_env)
311
312 CALL timestop(handle)
313
314 END SUBROUTINE create_a
315
316! **************************************************************************************************
317!> \brief Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W)
318!> B_ia,jb = α * v_ia,jb - W_ib,aj
319!> α is a spin-dependent factor
320!> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
321!> W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
322!> \param fm_mat_S_ia_bse ...
323!> \param fm_mat_S_bar_ia_bse ...
324!> \param fm_B ...
325!> \param homo ...
326!> \param virtual ...
327!> \param dimen_RI ...
328!> \param unit_nr ...
329!> \param mp2_env ...
330! **************************************************************************************************
331 SUBROUTINE create_b(fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, fm_B, &
332 homo, virtual, dimen_RI, unit_nr, mp2_env)
333
334 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_s_ia_bse, fm_mat_s_bar_ia_bse
335 TYPE(cp_fm_type), INTENT(INOUT) :: fm_b
336 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
337 INTEGER, INTENT(IN) :: dimen_ri, unit_nr
338 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
339
340 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_B'
341
342 INTEGER :: handle, isp, n_ov_joint, nspins
343 INTEGER, ALLOCATABLE, DIMENSION(:) :: n_ov, offsets
344 INTEGER, DIMENSION(4) :: reordering
345 REAL(kind=dp) :: alpha, alpha_screening
346 TYPE(cp_fm_struct_type), POINTER :: fm_struct_b, fm_struct_s_joint, &
347 fm_struct_w
348 TYPE(cp_fm_type) :: fm_s_joint, fm_w
349
350 CALL timeset(routinen, handle)
351
352 nspins = SIZE(homo)
353 ALLOCATE (n_ov(nspins), offsets(nspins))
354 CALL get_bse_spin_block_layout(homo, virtual, n_ov, offsets, n_ov_joint)
355
356 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
357 WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating B'
358 END IF
359
360 ! Coulomb prefactor: SPIN_CONFIG sector for closed shell; for open shell each spin block
361 ! contributes once (alpha=1). create_A already emits the SPIN_CONFIG-ignored warning.
362 SELECT CASE (mp2_env%bse%bse_spin_config)
363 CASE (bse_singlet)
364 alpha = 2.0_dp
365 CASE (bse_triplet)
366 alpha = 0.0_dp
367 END SELECT
368 IF (nspins > 1) alpha = 1.0_dp
369
370 IF (mp2_env%bse%screening_method == bse_screening_alpha) THEN
371 alpha_screening = mp2_env%bse%screening_factor
372 ELSE
373 alpha_screening = 1.0_dp
374 END IF
375
376 ! Joint B over all spin blocks: B_ia,jb = alpha*(ia|bj) - W^sigma_ib,aj (W spin-diagonal)
377 NULLIFY (fm_struct_b)
378 CALL cp_fm_struct_create(fm_struct_b, context=fm_mat_s_ia_bse(1)%matrix_struct%context, &
379 nrow_global=n_ov_joint, ncol_global=n_ov_joint, &
380 para_env=fm_mat_s_ia_bse(1)%matrix_struct%para_env)
381 CALL cp_fm_create(fm_b, fm_struct_b, name="fm_B_iajb")
382 CALL cp_fm_set_all(fm_b, 0.0_dp)
383
384 ! Coulomb v_ia,jb = sum_P B^P_ia B^P_jb (= (ia|bj)); cross-spin blocks filled automatically.
385 IF (nspins > 1) THEN
386 NULLIFY (fm_struct_s_joint)
387 CALL cp_fm_struct_create(fm_struct_s_joint, &
388 context=fm_mat_s_ia_bse(1)%matrix_struct%context, &
389 nrow_global=dimen_ri, ncol_global=n_ov_joint, &
390 para_env=fm_mat_s_ia_bse(1)%matrix_struct%para_env)
391 CALL cp_fm_create(fm_s_joint, fm_struct_s_joint, name="fm_S_ia_joint")
392 CALL cp_fm_set_all(fm_s_joint, 0.0_dp)
393 CALL assemble_joint_ov_slab(fm_mat_s_ia_bse, offsets, n_ov, dimen_ri, fm_s_joint)
394 CALL parallel_gemm(transa="T", transb="N", m=n_ov_joint, n=n_ov_joint, k=dimen_ri, &
395 alpha=alpha, matrix_a=fm_s_joint, matrix_b=fm_s_joint, beta=0.0_dp, &
396 matrix_c=fm_b)
397 CALL cp_fm_release(fm_s_joint)
398 CALL cp_fm_struct_release(fm_struct_s_joint)
399 ELSE
400 CALL parallel_gemm(transa="T", transb="N", m=homo(1)*virtual(1), n=homo(1)*virtual(1), &
401 k=dimen_ri, alpha=alpha, &
402 matrix_a=fm_mat_s_ia_bse(1), matrix_b=fm_mat_s_ia_bse(1), &
403 beta=0.0_dp, matrix_c=fm_b)
404 END IF
405
406 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
407 WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated B_iajb'
408 END IF
409
410 ! W^sigma_ib,aj = sum_P barB^P_ib B^P_aj on sigma-diagonal blocks only (offsets place them).
411 ! reordering [1,4,3,2] maps W_ib,ja -> B_ia,jb. For nspins=1: offsets(1)=0 (original code).
412 IF (mp2_env%bse%screening_method /= bse_screening_rpa) THEN
413 DO isp = 1, nspins
414 NULLIFY (fm_struct_w)
415 CALL cp_fm_struct_create(fm_struct_w, &
416 context=fm_mat_s_ia_bse(isp)%matrix_struct%context, &
417 nrow_global=homo(isp)*virtual(isp), &
418 ncol_global=homo(isp)*virtual(isp), &
419 para_env=fm_mat_s_ia_bse(isp)%matrix_struct%para_env)
420 CALL cp_fm_create(fm_w, fm_struct_w, name="fm_W_ibaj")
421 CALL cp_fm_set_all(fm_w, 0.0_dp)
422 CALL parallel_gemm(transa="T", transb="N", m=homo(isp)*virtual(isp), &
423 n=homo(isp)*virtual(isp), k=dimen_ri, alpha=alpha_screening, &
424 matrix_a=fm_mat_s_bar_ia_bse(isp), matrix_b=fm_mat_s_ia_bse(isp), &
425 beta=0.0_dp, matrix_c=fm_w)
426 reordering = [1, 4, 3, 2]
427 CALL fm_general_add_bse(fm_b, fm_w, -1.0_dp, virtual(isp), virtual(isp), &
428 virtual(isp), virtual(isp), unit_nr, reordering, mp2_env, &
429 row_offset=offsets(isp), col_offset=offsets(isp))
430 CALL cp_fm_release(fm_w)
431 CALL cp_fm_struct_release(fm_struct_w)
432 END DO
433 END IF
434
435 CALL cp_fm_struct_release(fm_struct_b)
436 DEALLOCATE (n_ov, offsets)
437
438 CALL timestop(handle)
439
440 END SUBROUTINE create_b
441
442 ! **************************************************************************************************
443!> \brief Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
444!> (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
445!> We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
446!> of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
447!> \param fm_A ...
448!> \param fm_B ...
449!> \param fm_C ...
450!> \param fm_sqrt_A_minus_B ...
451!> \param fm_inv_sqrt_A_minus_B ...
452!> \param unit_nr ...
453!> \param mp2_env ...
454!> \param diag_est ...
455! **************************************************************************************************
456 SUBROUTINE create_hermitian_form_of_abba(fm_A, fm_B, fm_C, &
457 fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
458 unit_nr, mp2_env, diag_est)
459
460 TYPE(cp_fm_type), INTENT(IN) :: fm_a, fm_b
461 TYPE(cp_fm_type), INTENT(INOUT) :: fm_c, fm_sqrt_a_minus_b, &
462 fm_inv_sqrt_a_minus_b
463 INTEGER, INTENT(IN) :: unit_nr
464 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
465 REAL(kind=dp), INTENT(IN) :: diag_est
466
467 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_hermitian_form_of_ABBA'
468
469 INTEGER :: dim_mat, handle, n_dependent
470 REAL(kind=dp), DIMENSION(2) :: eigvals_ab_diff
471 TYPE(cp_fm_type) :: fm_a_minus_b, fm_a_plus_b, fm_dummy, &
472 fm_work_product
473
474 CALL timeset(routinen, handle)
475
476 IF (unit_nr > 0) THEN
477 WRITE (unit_nr, '(T2,A4,T7,A25,A39,ES6.0,A3)') 'BSE|', 'Diagonalizing aux. matrix', &
478 ' with size of A. This will take around ', diag_est, " s."
479 END IF
480
481 ! Create work matrices, which will hold A+B and A-B and their powers
482 ! C is created afterwards to save memory
483 ! Final result: C = (A-B)^0.5 (A+B) (A-B)^0.5 EQ.I
484 ! \_______/ \___/ \______/
485 ! fm_sqrt_A_minus_B fm_A_plus_B fm_sqrt_A_minus_B
486 ! (EQ.Ia) (EQ.Ib) (EQ.Ia)
487 ! Intermediate work matrices:
488 ! fm_inv_sqrt_A_minus_B: (A-B)^-0.5 EQ.II
489 ! fm_A_minus_B: (A-B) EQ.III
490 ! fm_work_product: (A-B)^0.5 (A+B) from (EQ.Ia) and (EQ.Ib) EQ.IV
491 CALL cp_fm_create(fm_a_plus_b, fm_a%matrix_struct)
492 CALL cp_fm_to_fm(fm_a, fm_a_plus_b)
493 CALL cp_fm_create(fm_a_minus_b, fm_a%matrix_struct)
494 CALL cp_fm_to_fm(fm_a, fm_a_minus_b)
495 CALL cp_fm_create(fm_sqrt_a_minus_b, fm_a%matrix_struct)
496 CALL cp_fm_set_all(fm_sqrt_a_minus_b, 0.0_dp)
497 CALL cp_fm_create(fm_inv_sqrt_a_minus_b, fm_a%matrix_struct)
498 CALL cp_fm_set_all(fm_inv_sqrt_a_minus_b, 0.0_dp)
499
500 CALL cp_fm_create(fm_work_product, fm_a%matrix_struct)
501
502 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
503 WRITE (unit_nr, '(T2,A10,T13,A19)') 'BSE|DEBUG|', 'Created work arrays'
504 END IF
505
506 ! Add/Substract B (cf. EQs. Ib and III)
507 CALL cp_fm_scale_and_add(1.0_dp, fm_a_plus_b, 1.0_dp, fm_b)
508 CALL cp_fm_scale_and_add(1.0_dp, fm_a_minus_b, -1.0_dp, fm_b)
509
510 ! cp_fm_power will overwrite matrix, therefore we create copies
511 CALL cp_fm_to_fm(fm_a_minus_b, fm_inv_sqrt_a_minus_b)
512
513 ! In order to avoid a second diagonalization (cp_fm_power), we create (A-B)^0.5 (EQ.Ia)
514 ! from (A-B)^-0.5 (EQ.II) by multiplication with (A-B) (EQ.III) afterwards.
515
516 ! Raise A-B to -0.5_dp, no quenching of eigenvectors, hence threshold=0.0_dp
517 CALL cp_fm_create(fm_dummy, fm_a%matrix_struct)
518 ! Create (A-B)^-0.5 (cf. EQ.II)
519 CALL cp_fm_power(fm_inv_sqrt_a_minus_b, fm_dummy, -0.5_dp, 0.0_dp, n_dependent, eigvals=eigvals_ab_diff)
520 CALL cp_fm_release(fm_dummy)
521 ! Raise an error in case the the matrix A-B is not positive definite (i.e. negative eigenvalues)
522 ! In this case, the procedure for hermitian form of ABBA is not applicable
523 IF (eigvals_ab_diff(1) < 0) THEN
524 CALL cp_abort(__location__, &
525 "Matrix (A-B) is not positive definite. "// &
526 "Hermitian diagonalization of full ABBA matrix is ill-defined.")
527 END IF
528
529 ! We keep fm_inv_sqrt_A_minus_B for print of singleparticle transitions of ABBA
530 ! We further create (A-B)^0.5 for the singleparticle transitions of ABBA
531 ! Create (A-B)^0.5= (A-B)^-0.5 * (A-B) (EQ.Ia)
532 CALL cp_fm_get_info(fm_a, nrow_global=dim_mat)
533 CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_inv_sqrt_a_minus_b, fm_a_minus_b, 0.0_dp, &
534 fm_sqrt_a_minus_b)
535
536 ! Compute and store LHS of C, i.e. (A-B)^0.5 (A+B) (EQ.IV)
537 CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_sqrt_a_minus_b, fm_a_plus_b, 0.0_dp, &
538 fm_work_product)
539
540 ! Release to save memory
541 CALL cp_fm_release(fm_a_plus_b)
542 CALL cp_fm_release(fm_a_minus_b)
543
544 ! Now create full
545 CALL cp_fm_create(fm_c, fm_a%matrix_struct)
546 CALL cp_fm_set_all(fm_c, 0.0_dp)
547 ! Compute C=(A-B)^0.5 (A+B) (A-B)^0.5 (EQ.I)
548 CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_work_product, fm_sqrt_a_minus_b, 0.0_dp, &
549 fm_c)
550 CALL cp_fm_release(fm_work_product)
551
552 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
553 WRITE (unit_nr, '(T2,A10,T13,A36)') 'BSE|DEBUG|', 'Filled C=(A-B)^0.5 (A+B) (A-B)^0.5'
554 END IF
555
556 CALL timestop(handle)
557 END SUBROUTINE create_hermitian_form_of_abba
558
559! **************************************************************************************************
560!> \brief Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
561!> Here, the eigenvectors Z^n relate to X^n via
562!> Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
563!> \param fm_C ...
564!> \param homo ...
565!> \param virtual ...
566!> \param homo_irred ...
567!> \param fm_sqrt_A_minus_B ...
568!> \param fm_inv_sqrt_A_minus_B ...
569!> \param unit_nr ...
570!> \param diag_est ...
571!> \param mp2_env ...
572!> \param qs_env ...
573!> \param mo_coeff ...
574! **************************************************************************************************
575 SUBROUTINE diagonalize_c(fm_C, homo, virtual, homo_irred, &
576 fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
577 unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
578
579 TYPE(cp_fm_type), INTENT(INOUT) :: fm_c
580 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual, homo_irred
581 TYPE(cp_fm_type), INTENT(INOUT) :: fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b
582 INTEGER, INTENT(IN) :: unit_nr
583 REAL(kind=dp), INTENT(IN) :: diag_est
584 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
585 TYPE(qs_environment_type), POINTER :: qs_env
586 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
587
588 CHARACTER(LEN=*), PARAMETER :: routinen = 'diagonalize_C'
589
590 INTEGER :: diag_info, handle, n_ov_joint, nspins
591 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exc_ens
592 TYPE(cp_fm_type) :: fm_eigvec_x, fm_eigvec_y, fm_eigvec_z, &
593 fm_mat_eigvec_transform_diff, &
594 fm_mat_eigvec_transform_sum
595
596 CALL timeset(routinen, handle)
597
598 nspins = SIZE(homo)
599 n_ov_joint = sum(homo*virtual)
600
601 IF (unit_nr > 0) THEN
602 WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing C. ', &
603 'This will take around ', diag_est, ' s.'
604 END IF
605
606 !We have now the full matrix C=(A-B)^0.5 (A+B) (A-B)^0.5
607 !Now: Diagonalize it
608 CALL cp_fm_create(fm_eigvec_z, fm_c%matrix_struct)
609
610 ALLOCATE (exc_ens(n_ov_joint))
611
612 CALL choose_eigv_solver(fm_c, fm_eigvec_z, exc_ens, diag_info)
613
614 IF (diag_info /= 0) THEN
615 CALL cp_abort(__location__, &
616 "Diagonalization of C=(A-B)^0.5 (A+B) (A-B)^0.5 failed in BSE")
617 END IF
618
619 ! C could have negative eigenvalues, since we do not explicitly check A+B
620 ! for positive definiteness (would make another O(N^6) Diagon. necessary)
621 ! Instead, we include a check here
622 IF (exc_ens(1) < 0) THEN
623 IF (unit_nr > 0) THEN
624 CALL cp_abort(__location__, &
625 "Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 has negative eigenvalues, i.e. "// &
626 "(A+B) is not positive definite.")
627 END IF
628 END IF
629 exc_ens = sqrt(exc_ens)
630
631 ! Prepare eigenvector for interpretation of singleparticle transitions
632 ! Compare: F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)
633 ! We aim for the upper part of the vector (X,Y) for a direct comparison with the TDA result
634
635 ! Following Furche, we basically use Eqs. (A10): First, we multiply
636 ! the (A-B)^+-0.5 with eigenvectors and then the eigenvalues
637 ! One has to be careful about the index structure, since the eigenvector matrix is not symmetric anymore!
638
639 ! First, Eq. I from (A10) from Furche: (X+Y)_n = (Ω_n)^-0.5 (A-B)^0.5 T_n
640 CALL cp_fm_create(fm_mat_eigvec_transform_sum, fm_c%matrix_struct)
641 CALL cp_fm_set_all(fm_mat_eigvec_transform_sum, 0.0_dp)
642 CALL parallel_gemm(transa="N", transb="N", m=n_ov_joint, n=n_ov_joint, k=n_ov_joint, alpha=1.0_dp, &
643 matrix_a=fm_sqrt_a_minus_b, matrix_b=fm_eigvec_z, beta=0.0_dp, &
644 matrix_c=fm_mat_eigvec_transform_sum)
645 CALL cp_fm_release(fm_sqrt_a_minus_b)
646 ! This normalizes the eigenvectors
647 CALL comp_eigvec_coeff_bse(fm_mat_eigvec_transform_sum, exc_ens, -0.5_dp, gamma=2.0_dp, do_transpose=.true.)
648
649 ! Second, Eq. II from (A10) from Furche: (X-Y)_n = (Ω_n)^0.5 (A-B)^-0.5 T_n
650 CALL cp_fm_create(fm_mat_eigvec_transform_diff, fm_c%matrix_struct)
651 CALL cp_fm_set_all(fm_mat_eigvec_transform_diff, 0.0_dp)
652 CALL parallel_gemm(transa="N", transb="N", m=n_ov_joint, n=n_ov_joint, k=n_ov_joint, alpha=1.0_dp, &
653 matrix_a=fm_inv_sqrt_a_minus_b, matrix_b=fm_eigvec_z, beta=0.0_dp, &
654 matrix_c=fm_mat_eigvec_transform_diff)
655 CALL cp_fm_release(fm_inv_sqrt_a_minus_b)
656 CALL cp_fm_release(fm_eigvec_z)
657
658 ! This normalizes the eigenvectors
659 CALL comp_eigvec_coeff_bse(fm_mat_eigvec_transform_diff, exc_ens, 0.5_dp, gamma=2.0_dp, do_transpose=.true.)
660
661 ! Now, we add the two equations to obtain X_n
662 ! Add overwrites the first argument, therefore we copy it beforehand
663 CALL cp_fm_create(fm_eigvec_x, fm_c%matrix_struct)
664 CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_x)
665 CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_x, 1.0_dp, fm_mat_eigvec_transform_diff)
666
667 ! Now, we subtract the two equations to obtain Y_n
668 ! Add overwrites the first argument, therefore we copy it beforehand
669 CALL cp_fm_create(fm_eigvec_y, fm_c%matrix_struct)
670 CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_y)
671 CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_y, -1.0_dp, fm_mat_eigvec_transform_diff)
672
673 !Cleanup
674 CALL cp_fm_release(fm_mat_eigvec_transform_diff)
675 CALL cp_fm_release(fm_mat_eigvec_transform_sum)
676
677 IF (nspins == 1) THEN
678 CALL postprocess_bse(exc_ens, fm_eigvec_x, mp2_env, qs_env, mo_coeff, &
679 homo(1), virtual(1), homo_irred(1), unit_nr, &
680 .false., fm_eigvec_y)
681 ELSE
682 ! Open-shell ABBA: helper forms X+Y internally and prints amplitudes (X and Y).
683 CALL bse_open_shell_optical(exc_ens, fm_eigvec_x, homo, virtual, homo_irred, &
684 .false., qs_env, mo_coeff, mp2_env, unit_nr, fm_eigvec_y)
685 END IF
686
687 DEALLOCATE (exc_ens)
688 CALL cp_fm_release(fm_eigvec_x)
689 CALL cp_fm_release(fm_eigvec_y)
690
691 CALL timestop(handle)
692
693 END SUBROUTINE diagonalize_c
694
695! **************************************************************************************************
696!> \brief Solving hermitian eigenvalue equation A X^n = Ω^n X^n
697!> \param fm_A ...
698!> \param homo ...
699!> \param virtual ...
700!> \param homo_irred ...
701!> \param unit_nr ...
702!> \param diag_est ...
703!> \param mp2_env ...
704!> \param qs_env ...
705!> \param mo_coeff ...
706! **************************************************************************************************
707 SUBROUTINE diagonalize_a(fm_A, homo, virtual, homo_irred, &
708 unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
709
710 TYPE(cp_fm_type), INTENT(INOUT) :: fm_a
711 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual, homo_irred
712 INTEGER, INTENT(IN) :: unit_nr
713 REAL(kind=dp), INTENT(IN) :: diag_est
714 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
715 TYPE(qs_environment_type), POINTER :: qs_env
716 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
717
718 CHARACTER(LEN=*), PARAMETER :: routinen = 'diagonalize_A'
719
720 INTEGER :: diag_info, handle, n_ov_joint, nspins
721 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exc_ens
722 TYPE(cp_fm_type) :: fm_eigvec
723
724 CALL timeset(routinen, handle)
725
726 nspins = SIZE(homo)
727 n_ov_joint = sum(homo*virtual)
728
729 IF (unit_nr > 0) THEN
730 WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing A. ', &
731 'This will take around ', diag_est, ' s.'
732 END IF
733
734 CALL cp_fm_create(fm_eigvec, fm_a%matrix_struct)
735
736 ALLOCATE (exc_ens(n_ov_joint))
737
738 CALL choose_eigv_solver(fm_a, fm_eigvec, exc_ens, diag_info)
739
740 IF (diag_info /= 0) THEN
741 CALL cp_abort(__location__, &
742 "Diagonalization of A failed in TDA-BSE")
743 END IF
744
745 IF (nspins == 1) THEN
746 CALL postprocess_bse(exc_ens, fm_eigvec, mp2_env, qs_env, mo_coeff, &
747 homo(1), virtual(1), homo_irred(1), unit_nr, .true.)
748 ELSE
749 CALL bse_open_shell_optical(exc_ens, fm_eigvec, homo, virtual, homo_irred, &
750 .true., qs_env, mo_coeff, mp2_env, unit_nr)
751 END IF
752
753 CALL cp_fm_release(fm_eigvec)
754 DEALLOCATE (exc_ens)
755
756 CALL timestop(handle)
757
758 END SUBROUTINE diagonalize_a
759
760! **************************************************************************************************
761!> \brief Open-shell (UKS) spin-summed post-processing for the joint spin-block space: joint
762!> excitation energies, per-spin transition amplitudes, and oscillator strengths. Mirrors
763!> postprocess_bse but spin-summed; exciton descriptors and NTOs are not yet implemented (CPWARN).
764!> \param Exc_ens joint excitation energies
765!> \param fm_eigvec_X joint X eigenvectors (excitations)
766!> \param homo per-spin reduced/active occupied counts
767!> \param virtual per-spin reduced/active virtual counts
768!> \param homo_irred per-spin full occupied counts (absolute-MO labels; N_e = sum)
769!> \param flag_tda .TRUE. -> TDA (coeff=X), .FALSE. -> ABBA (coeff=X+Y)
770!> \param qs_env ...
771!> \param mo_coeff per-spin MO coefficients
772!> \param mp2_env ...
773!> \param unit_nr ...
774!> \param fm_eigvec_Y joint Y eigenvectors (deexcitations; ABBA only)
775! **************************************************************************************************
776 SUBROUTINE bse_open_shell_optical(Exc_ens, fm_eigvec_X, homo, virtual, homo_irred, &
777 flag_tda, qs_env, mo_coeff, mp2_env, unit_nr, fm_eigvec_Y)
778
779 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exc_ens
780 TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_x
781 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual, homo_irred
782 LOGICAL, INTENT(IN) :: flag_tda
783 TYPE(qs_environment_type), POINTER :: qs_env
784 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
785 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
786 INTEGER, INTENT(IN) :: unit_nr
787 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_y
788
789 CHARACTER(LEN=*), PARAMETER :: routinen = 'bse_open_shell_optical'
790
791 CHARACTER(LEN=10) :: info_approximation, multiplet
792 INTEGER :: handle, idir, isp, jdir, n, n_ov_joint, &
793 nspins
794 INTEGER, ALLOCATABLE, DIMENSION(:) :: n_ov_sp, offsets_sp
795 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str_joint, ref_pt
796 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pol_res_joint, trans_mom_joint
797 TYPE(cp_fm_struct_type), POINTER :: fm_struct_dip_reord, fm_struct_sp, &
798 fm_struct_tmom
799 TYPE(cp_fm_type) :: fm_dip_reord_sp, fm_eigvec_sp, &
800 fm_trans_coeff
801 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dip_ab_sp, fm_dip_ai_sp, fm_dip_ij_sp
802 TYPE(cp_fm_type), DIMENSION(3) :: fm_trans_mom_joint
803
804 CALL timeset(routinen, handle)
805
806 nspins = SIZE(homo)
807 ALLOCATE (n_ov_sp(nspins), offsets_sp(nspins))
808 CALL get_bse_spin_block_layout(homo, virtual, n_ov_sp, offsets_sp, n_ov_joint)
809
810 ! LEN=10 locals auto-pad short literals with spaces (avoids the L-15 short-literal trap);
811 ! print_excitation_energies prints A6 of these, print_optical_properties prints them as-is.
812 multiplet = "UKS"
813 IF (flag_tda) THEN
814 info_approximation = " -TDA- "
815 ELSE
816 info_approximation = "-ABBA-"
817 END IF
818
819 IF (unit_nr > 0) THEN
820 WRITE (unit_nr, '(T2,A4,T7,A43)') 'BSE|', 'Joint open-shell BSE excitation energies:'
821 END IF
822 CALL print_excitation_energies(exc_ens, n_ov_joint, 1, flag_tda, multiplet, &
823 info_approximation, mp2_env, unit_nr)
824
825 ! Per-spin single-particle transition amplitudes (X via =>, Y via <=).
826 CALL print_transition_amplitudes(fm_eigvec_x, homo, virtual, homo_irred, &
827 info_approximation, mp2_env, unit_nr, fm_eigvec_y)
828
829 ! Transition coefficient for the spin-summed moment: X (TDA) or X+Y (ABBA).
830 CALL cp_fm_create(fm_trans_coeff, fm_eigvec_x%matrix_struct)
831 CALL cp_fm_to_fm(fm_eigvec_x, fm_trans_coeff)
832 IF (PRESENT(fm_eigvec_y)) CALL cp_fm_scale_and_add(1.0_dp, fm_trans_coeff, 1.0_dp, fm_eigvec_y)
833
834 ! Spin-summed transition moments: D^n_dir = sum_σ sum_{ia,σ} D^{dir,σ}_{ai} C_{ia,σ,n}
835 ! with C = X (TDA) or X+Y (ABBA); explicit spin sum, factor 1.0.
836 ALLOCATE (fm_dip_ai_sp(3), fm_dip_ij_sp(3), fm_dip_ab_sp(3), ref_pt(3))
837 ALLOCATE (oscill_str_joint(n_ov_joint), trans_mom_joint(3, 1, n_ov_joint))
838 ALLOCATE (pol_res_joint(3, 3, n_ov_joint))
839 trans_mom_joint(:, :, :) = 0.0_dp
840 NULLIFY (fm_struct_dip_reord, fm_struct_sp, fm_struct_tmom)
841 CALL cp_fm_struct_create(fm_struct_tmom, fm_trans_coeff%matrix_struct%para_env, &
842 fm_trans_coeff%matrix_struct%context, 1, n_ov_joint)
843 DO idir = 1, 3
844 CALL cp_fm_create(fm_trans_mom_joint(idir), fm_struct_tmom)
845 CALL cp_fm_set_all(fm_trans_mom_joint(idir), 0.0_dp)
846 END DO
847 DO isp = 1, nspins
848 CALL get_multipoles_mo(fm_dip_ai_sp, fm_dip_ij_sp, fm_dip_ab_sp, &
849 qs_env, mo_coeff(isp:isp), ref_pt, 1, &
850 homo(isp), virtual(isp), fm_trans_coeff%matrix_struct%context, &
851 ispin=isp)
852 NULLIFY (fm_struct_sp, fm_struct_dip_reord)
853 CALL cp_fm_struct_create(fm_struct_sp, fm_trans_coeff%matrix_struct%para_env, &
854 fm_trans_coeff%matrix_struct%context, n_ov_sp(isp), n_ov_joint)
855 CALL cp_fm_create(fm_eigvec_sp, fm_struct_sp)
856 CALL cp_fm_set_all(fm_eigvec_sp, 0.0_dp)
857 CALL cp_fm_to_fm_submat(fm_trans_coeff, fm_eigvec_sp, n_ov_sp(isp), n_ov_joint, &
858 offsets_sp(isp) + 1, 1, 1, 1)
859 CALL cp_fm_struct_create(fm_struct_dip_reord, fm_trans_coeff%matrix_struct%para_env, &
860 fm_trans_coeff%matrix_struct%context, 1, n_ov_sp(isp))
861 DO idir = 1, 3
862 CALL cp_fm_create(fm_dip_reord_sp, fm_struct_dip_reord, name="bse_dip_reord")
863 CALL cp_fm_set_all(fm_dip_reord_sp, 0.0_dp)
864 CALL fm_general_add_bse(fm_dip_reord_sp, fm_dip_ai_sp(idir), 1.0_dp, &
865 1, 1, 1, virtual(isp), unit_nr, [2, 4, 3, 1], mp2_env)
866 CALL parallel_gemm('N', 'N', 1, n_ov_joint, n_ov_sp(isp), 1.0_dp, &
867 fm_dip_reord_sp, fm_eigvec_sp, 1.0_dp, fm_trans_mom_joint(idir))
868 CALL cp_fm_release(fm_dip_reord_sp)
869 CALL cp_fm_release(fm_dip_ai_sp(idir))
870 CALL cp_fm_release(fm_dip_ij_sp(idir))
871 CALL cp_fm_release(fm_dip_ab_sp(idir))
872 END DO
873 CALL cp_fm_release(fm_eigvec_sp)
874 CALL cp_fm_struct_release(fm_struct_sp)
875 NULLIFY (fm_struct_sp)
876 CALL cp_fm_struct_release(fm_struct_dip_reord)
877 NULLIFY (fm_struct_dip_reord)
878 END DO
879 DO idir = 1, 3
880 CALL cp_fm_get_submatrix(fm_trans_mom_joint(idir), trans_mom_joint(idir, :, :))
881 CALL cp_fm_release(fm_trans_mom_joint(idir))
882 END DO
883 CALL cp_fm_struct_release(fm_struct_tmom)
884 DO n = 1, n_ov_joint
885 DO idir = 1, 3
886 DO jdir = 1, 3
887 pol_res_joint(idir, jdir, n) = 2.0_dp*exc_ens(n)*trans_mom_joint(idir, 1, n) &
888 *trans_mom_joint(jdir, 1, n)
889 END DO
890 END DO
891 oscill_str_joint(n) = 2.0_dp/3.0_dp*exc_ens(n)*sum(abs(trans_mom_joint(:, 1, n))**2)
892 END DO
893 CALL print_optical_properties(exc_ens, oscill_str_joint, trans_mom_joint, pol_res_joint, &
894 n_ov_joint, 1, sum(homo_irred), flag_tda, info_approximation, &
895 mp2_env, unit_nr, open_shell=.true.)
896 ! Open-shell post-processing is partial: energies, amplitudes, spin-summed oscillator strengths.
897 CALL cp_warn(__location__, &
898 "Open-shell (UKS) BSE: exciton descriptors and NTO analysis are not yet "// &
899 "implemented and have been skipped.")
900 CALL cp_fm_release(fm_trans_coeff)
901 DEALLOCATE (fm_dip_ai_sp, fm_dip_ij_sp, fm_dip_ab_sp, ref_pt)
902 DEALLOCATE (n_ov_sp, offsets_sp, oscill_str_joint, trans_mom_joint, pol_res_joint)
903
904 CALL timestop(handle)
905
906 END SUBROUTINE bse_open_shell_optical
907
908! **************************************************************************************************
909!> \brief Prints the success message (incl. energies) for full diag of BSE (TDA/full ABBA via flag)
910!> \param Exc_ens ...
911!> \param fm_eigvec_X ...
912!> \param mp2_env ...
913!> \param qs_env ...
914!> \param mo_coeff ...
915!> \param homo ...
916!> \param virtual ...
917!> \param homo_irred ...
918!> \param unit_nr ...
919!> \param flag_TDA ...
920!> \param fm_eigvec_Y ...
921! **************************************************************************************************
922 SUBROUTINE postprocess_bse(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &
923 homo, virtual, homo_irred, unit_nr, &
924 flag_TDA, fm_eigvec_Y)
925
926 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exc_ens
927 TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_x
928 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
929 TYPE(qs_environment_type), POINTER :: qs_env
930 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
931 INTEGER :: homo, virtual, homo_irred, unit_nr
932 LOGICAL, OPTIONAL :: flag_tda
933 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_y
934
935 CHARACTER(LEN=*), PARAMETER :: routinen = 'postprocess_bse'
936
937 CHARACTER(LEN=10) :: info_approximation, multiplet
938 INTEGER :: handle, i_exc, idir, n_moments_di, &
939 n_moments_quad
940 REAL(kind=dp) :: alpha
941 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str, ref_point_multipole
942 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: polarizability_residues, trans_mom_bse
943 TYPE(cp_fm_type) :: fm_x_ia, fm_y_ia
944 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dipole_ab_trunc, fm_dipole_ai_trunc, &
945 fm_dipole_ij_trunc, fm_quadpole_ab_trunc, fm_quadpole_ai_trunc, fm_quadpole_ij_trunc
946 TYPE(exciton_descr_type), ALLOCATABLE, &
947 DIMENSION(:) :: exc_descr
948
949 CALL timeset(routinen, handle)
950
951 !Prepare variables for printing
952 IF (mp2_env%bse%bse_spin_config == 0) THEN
953 multiplet = "Singlet"
954 alpha = 2.0_dp
955 ELSE
956 multiplet = "Triplet"
957 alpha = 0.0_dp
958 END IF
959 IF (.NOT. PRESENT(flag_tda)) THEN
960 flag_tda = .false.
961 END IF
962 IF (flag_tda) THEN
963 info_approximation = " -TDA- "
964 ELSE
965 info_approximation = "-ABBA-"
966 END IF
967
968 n_moments_di = 3
969 n_moments_quad = 9
970 ! Compute BSE dipoles and oscillator strengths - Keep in memory for later usage
971 ! Need dipoles also for spatial expectation values, which are well-defined also for triplets
972 ALLOCATE (fm_dipole_ij_trunc(n_moments_di))
973 ALLOCATE (fm_dipole_ab_trunc(n_moments_di))
974 ALLOCATE (fm_dipole_ai_trunc(n_moments_di))
975 ALLOCATE (ref_point_multipole(3))
976 ! Obtain dipoles in MO basis
977 CALL get_multipoles_mo(fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc, &
978 qs_env, mo_coeff, ref_point_multipole, 1, &
979 homo, virtual, fm_eigvec_x%matrix_struct%context)
980 ! Compute exciton descriptors from these multipoles
981 IF (mp2_env%bse%num_print_exc_descr > 0) THEN
982 ! Obtain quadrupoles in MO basis
983 ALLOCATE (fm_quadpole_ij_trunc(n_moments_quad))
984 ALLOCATE (fm_quadpole_ab_trunc(n_moments_quad))
985 ALLOCATE (fm_quadpole_ai_trunc(n_moments_quad))
986 CALL get_multipoles_mo(fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
987 qs_env, mo_coeff, ref_point_multipole, 2, &
988 homo, virtual, fm_eigvec_x%matrix_struct%context)
989 ! Iterate over excitation index outside of routine to make it compatible with tddft module
990 ALLOCATE (exc_descr(mp2_env%bse%num_print_exc_descr))
991 DO i_exc = 1, mp2_env%bse%num_print_exc_descr
992 CALL reshuffle_eigvec(fm_eigvec_x, fm_x_ia, homo, virtual, i_exc, &
993 .false., unit_nr, mp2_env)
994 IF (.NOT. flag_tda) THEN
995 CALL reshuffle_eigvec(fm_eigvec_y, fm_y_ia, homo, virtual, i_exc, &
996 .false., unit_nr, mp2_env)
997
998 CALL get_exciton_descriptors(exc_descr, fm_x_ia, &
999 fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
1000 fm_quadpole_ai_trunc, &
1001 i_exc, homo, virtual, &
1002 fm_y_ia)
1003 ELSE
1004 CALL get_exciton_descriptors(exc_descr, fm_x_ia, &
1005 fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
1006 fm_quadpole_ai_trunc, &
1007 i_exc, homo, virtual)
1008 END IF
1009 CALL cp_fm_release(fm_x_ia)
1010 IF (.NOT. flag_tda) THEN
1011 CALL cp_fm_release(fm_y_ia)
1012 END IF
1013 END DO
1014 END IF
1015
1016 IF (mp2_env%bse%bse_spin_config == 0) THEN
1017 CALL get_oscillator_strengths(fm_eigvec_x, exc_ens, fm_dipole_ai_trunc, &
1018 trans_mom_bse, oscill_str, polarizability_residues, &
1019 mp2_env, homo, virtual, unit_nr, &
1020 fm_eigvec_y)
1021 END IF
1022
1023 ! Prints basic definitions used in BSE calculation
1024 CALL print_output_header(homo, virtual, homo_irred, flag_tda, &
1025 multiplet, alpha, mp2_env, unit_nr)
1026
1027 ! Prints excitation energies up to user-specified number
1028 CALL print_excitation_energies(exc_ens, homo, virtual, flag_tda, multiplet, &
1029 info_approximation, mp2_env, unit_nr)
1030
1031 ! Print single particle transition amplitudes, i.e. components of eigenvectors X and Y
1032 CALL print_transition_amplitudes(fm_eigvec_x, [homo], [virtual], [homo_irred], &
1033 info_approximation, mp2_env, unit_nr, fm_eigvec_y)
1034
1035 ! Prints optical properties, if state is a singlet
1036 CALL print_optical_properties(exc_ens, oscill_str, trans_mom_bse, polarizability_residues, &
1037 homo, virtual, homo_irred, flag_tda, &
1038 info_approximation, mp2_env, unit_nr)
1039 ! Print exciton descriptors if keyword is invoked
1040 IF (mp2_env%bse%num_print_exc_descr > 0) THEN
1041 CALL print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, &
1042 mp2_env%bse%num_print_exc_descr, mp2_env%bse%bse_debug_print, &
1043 mp2_env%bse%print_directional_exc_descr, &
1044 'BSE|', qs_env)
1045 END IF
1046
1047 ! Compute and print excitation wavefunctions
1048 IF (mp2_env%bse%do_nto_analysis) THEN
1049 IF (unit_nr > 0) THEN
1050 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1051 WRITE (unit_nr, '(T2,A4,T7,A47)') &
1052 'BSE|', "Calculating Natural Transition Orbitals (NTOs)."
1053 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1054 END IF
1055 CALL calculate_ntos(fm_eigvec_x, fm_eigvec_y, &
1056 mo_coeff, homo, virtual, &
1057 info_approximation, &
1058 oscill_str, &
1059 qs_env, unit_nr, mp2_env)
1060 END IF
1061
1062 DO idir = 1, n_moments_di
1063 CALL cp_fm_release(fm_dipole_ai_trunc(idir))
1064 CALL cp_fm_release(fm_dipole_ij_trunc(idir))
1065 CALL cp_fm_release(fm_dipole_ab_trunc(idir))
1066 END DO
1067 IF (mp2_env%bse%num_print_exc_descr > 0) THEN
1068 DO idir = 1, n_moments_quad
1069 CALL cp_fm_release(fm_quadpole_ai_trunc(idir))
1070 CALL cp_fm_release(fm_quadpole_ij_trunc(idir))
1071 CALL cp_fm_release(fm_quadpole_ab_trunc(idir))
1072 END DO
1073 DEALLOCATE (fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc)
1074 DEALLOCATE (exc_descr)
1075 END IF
1076 DEALLOCATE (fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc)
1077 DEALLOCATE (ref_point_multipole)
1078 IF (mp2_env%bse%bse_spin_config == 0) THEN
1079 DEALLOCATE (oscill_str, trans_mom_bse, polarizability_residues)
1080 END IF
1081 IF (unit_nr > 0) THEN
1082 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1083 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1084 END IF
1085
1086 CALL timestop(handle)
1087
1088 END SUBROUTINE postprocess_bse
1089
1090END MODULE bse_full_diag
Routines for the full diagonalization of GW + Bethe-Salpeter for computing electronic excitations.
subroutine, public diagonalize_a(fm_a, homo, virtual, homo_irred, unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
Solving hermitian eigenvalue equation A X^n = Ω^n X^n.
subroutine, public diagonalize_c(fm_c, homo, virtual, homo_irred, fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, unit_nr, diag_est, mp2_env, qs_env, mo_coeff)
Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n . Here, the eigenvectors Z^n relate to X^n via Eq....
subroutine, public create_hermitian_form_of_abba(fm_a, fm_b, fm_c, fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, unit_nr, mp2_env, diag_est)
Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem (cf....
subroutine, public create_b(fm_mat_s_ia_bse, fm_mat_s_bar_ia_bse, fm_b, homo, virtual, dimen_ri, unit_nr, mp2_env)
Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W) B_ia,jb = α * v_ia,...
subroutine, public create_a(fm_mat_s_ia_bse, fm_mat_s_bar_ij_bse, fm_mat_s_ab_bse, fm_a, eigenval, unit_nr, homo, virtual, dimen_ri, mp2_env, para_env, qs_env)
Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W) A_ia,...
Routines for printing information in context of the BSE calculation.
Definition bse_print.F:13
subroutine, public print_output_header(homo, virtual, homo_irred, flag_tda, multiplet, alpha, mp2_env, unit_nr)
...
Definition bse_print.F:102
subroutine, public print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, num_print_exc_descr, print_checkvalue, print_directional_exc_descr, prefix_output, qs_env)
Prints exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018).
Definition bse_print.F:573
subroutine, public print_optical_properties(exc_ens, oscill_str, trans_mom_bse, polarizability_residues, homo, virtual, homo_irred, flag_tda, info_approximation, mp2_env, unit_nr, open_shell)
...
Definition bse_print.F:366
subroutine, public print_transition_amplitudes(fm_eigvec_x, homo, virtual, homo_irred, info_approximation, mp2_env, unit_nr, fm_eigvec_y)
...
Definition bse_print.F:271
subroutine, public print_excitation_energies(exc_ens, homo, virtual, flag_tda, multiplet, info_approximation, mp2_env, unit_nr)
...
Definition bse_print.F:222
Routines for computing excitonic properties, e.g. exciton diameter, from the BSE.
subroutine, public get_oscillator_strengths(fm_eigvec_x, exc_ens, fm_dipole_ai_trunc, trans_mom_bse, oscill_str, polarizability_residues, mp2_env, homo_red, virtual_red, unit_nr, fm_eigvec_y)
Compute and return BSE dipoles d_r^n = sqrt(2) sum_ia < ψ_i | r | ψ_a > ( X_ia^n + Y_ia^n ) and oscil...
subroutine, public get_exciton_descriptors(exc_descr, fm_x_ia, fm_multipole_ij_trunc, fm_multipole_ab_trunc, fm_multipole_ai_trunc, i_exc, homo, virtual, fm_y_ia)
...
subroutine, public calculate_ntos(fm_x, fm_y, mo_coeff, homo, virtual, info_approximation, oscill_str, qs_env, unit_nr, mp2_env)
...
Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_util.F:13
subroutine, public assemble_joint_ov_slab(fm_s_ia, offsets, n_ov, dimen_ri, fm_s_joint)
Column-concatenate per-spin ia-slabs into the joint dimen_RI x n_ov_joint slab. Sigma-block of spin i...
Definition bse_util.F:2112
subroutine, public fm_general_add_bse(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, nrow_secidx_out, ncol_secidx_out, unit_nr, reordering, mp2_env, row_offset, col_offset)
Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,...
Definition bse_util.F:203
subroutine, public get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, qs_env, mo_coeff, rpoint, n_moments, homo_red, virtual_red, context_bse, ispin)
...
Definition bse_util.F:1886
subroutine, public reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose, unit_nr, mp2_env)
...
Definition bse_util.F:1521
subroutine, public get_bse_spin_block_layout(homo_red, virt_red, n_ov, offsets, n_ov_joint)
Spin-block layout for the open-shell (joint) BSE matrix: per-spin OV-pair counts and the block offset...
Definition bse_util.F:1262
subroutine, public comp_eigvec_coeff_bse(fm_work, eig_vals, beta, gamma, do_transpose)
Routine for computing the coefficients of the eigenvectors of the BSE matrix from a multiplication wi...
Definition bse_util.F:800
methods related to the blacs parallel environment
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:262
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
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
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Types for excited states potential energies.
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public bse_singlet
integer, parameter, public bse_triplet
integer, parameter, public bse_screening_alpha
integer, parameter, public bse_screening_rpa
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Interface to the message passing library MPI.
Types needed for MP2 calculations.
Definition mp2_types.F:14
basic linear algebra operations for full matrixes
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
Contains information on the excited states energy.
stores all the informations relevant to an mpi environment