(git:1f285aa)
xas_tdp_utils.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Utilities for X-ray absorption spectroscopy using TDDFPT
10 !> \author AB (01.2018)
11 ! **************************************************************************************************
12 
13 MODULE xas_tdp_utils
14  USE cp_blacs_env, ONLY: cp_blacs_env_type
15  USE cp_cfm_diag, ONLY: cp_cfm_heevd
16  USE cp_cfm_types, ONLY: cp_cfm_create,&
17  cp_cfm_get_info,&
18  cp_cfm_get_submatrix,&
19  cp_cfm_release,&
20  cp_cfm_type,&
21  cp_fm_to_cfm
22  USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
23  cp_dbcsr_cholesky_invert
24  USE cp_dbcsr_diag, ONLY: cp_dbcsr_power
25  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
26  copy_fm_to_dbcsr,&
27  cp_dbcsr_sm_fm_multiply,&
28  dbcsr_allocate_matrix_set,&
29  dbcsr_deallocate_matrix_set
30  USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
31  cp_fm_scale,&
32  cp_fm_transpose,&
33  cp_fm_upper_to_full
34  USE cp_fm_diag, ONLY: choose_eigv_solver,&
35  cp_fm_geeig
36  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
37  cp_fm_struct_release,&
38  cp_fm_struct_type
39  USE cp_fm_types, ONLY: cp_fm_create,&
40  cp_fm_get_diag,&
41  cp_fm_get_info,&
42  cp_fm_get_submatrix,&
43  cp_fm_release,&
44  cp_fm_set_element,&
45  cp_fm_to_fm_submat,&
46  cp_fm_type
47  USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
48  USE dbcsr_api, ONLY: &
49  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, &
50  dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_block_p, &
51  dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
52  dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
53  dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_reserve_all_blocks, dbcsr_set, &
54  dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
55  USE input_constants, ONLY: ot_precond_full_single,&
56  tddfpt_singlet,&
57  tddfpt_spin_cons,&
58  tddfpt_spin_flip,&
59  tddfpt_triplet,&
60  xas_dip_len
61  USE kinds, ONLY: dp
62  USE mathlib, ONLY: get_diag
63  USE message_passing, ONLY: mp_para_env_type
64  USE parallel_gemm_api, ONLY: parallel_gemm
65  USE physcon, ONLY: a_fine
66  USE preconditioner_types, ONLY: destroy_preconditioner,&
67  init_preconditioner,&
68  preconditioner_type
69  USE qs_environment_types, ONLY: get_qs_env,&
70  qs_environment_type
71  USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
72  USE qs_mo_types, ONLY: get_mo_set,&
73  mo_set_type
74  USE qs_ot_eigensolver, ONLY: ot_eigensolver
75  USE xas_tdp_kernel, ONLY: kernel_coulomb_xc,&
76  kernel_exchange
77  USE xas_tdp_types, ONLY: donor_state_type,&
78  xas_tdp_control_type,&
79  xas_tdp_env_type
80 
81 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
82 #include "./base/base_uses.f90"
83 
84  IMPLICIT NONE
85  PRIVATE
86 
87  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_tdp_utils'
88 
89  PUBLIC :: setup_xas_tdp_prob, solve_xas_tdp_prob, include_rcs_soc, &
90  include_os_soc, rcs_amew_soc_elements
91 
92  !A helper type for SOC
93  TYPE dbcsr_soc_package_type
94  TYPE(dbcsr_type), POINTER :: dbcsr_sg => null()
95  TYPE(dbcsr_type), POINTER :: dbcsr_tp => null()
96  TYPE(dbcsr_type), POINTER :: dbcsr_sc => null()
97  TYPE(dbcsr_type), POINTER :: dbcsr_sf => null()
98  TYPE(dbcsr_type), POINTER :: dbcsr_prod => null()
99  TYPE(dbcsr_type), POINTER :: dbcsr_ovlp => null()
100  TYPE(dbcsr_type), POINTER :: dbcsr_tmp => null()
101  TYPE(dbcsr_type), POINTER :: dbcsr_work => null()
102  END TYPE dbcsr_soc_package_type
103 
104 CONTAINS
105 
106 ! **************************************************************************************************
107 !> \brief Builds the matrix that defines the XAS TDDFPT generalized eigenvalue problem to be solved
108 !> for excitation energies omega. The problem has the form omega*G*C = M*C, where C contains
109 !> the response orbitals coefficients. The matrix M and the metric G are stored in the given
110 !> donor_state
111 !> \param donor_state the donor_state for which the problem is restricted
112 !> \param qs_env ...
113 !> \param xas_tdp_env ...
114 !> \param xas_tdp_control ...
115 !> \note the matrix M is symmetric and has the form | M_d M_o |
116 !> | M_o M_d |,
117 !> -In the SPIN-RESTRICTED case:
118 !> depending on whether we consider singlet or triplet excitation, the diagonal (M_d) and
119 !> off-diagonal (M_o) parts of M differ:
120 !> - For singlet: M_d = A + 2B + C_aa + C_ab - D
121 !> M_o = 2B + C_aa + C_ab - E
122 !> - For triplet: M_d = A + C_aa - C_ab - D
123 !> M_o = C_aa - C_ab - E
124 !> where other subroutines computes the matrices A, B, E, D and G, which are:
125 !> - A: the ground-state contribution: F_pq*delta_IJ - epsilon_IJ*S_pq
126 !> - B: the Coulomb kernel ~(pI|Jq)
127 !> - C: the xc kernel c_aa (double derivatibe wrt to n_alpha) and C_ab (wrt n_alpha and n_beta)
128 !> - D: the on-digonal exact exchange kernel ~(pq|IJ)
129 !> - E: the off-diagonal exact exchange kernel ~(pJ|Iq)
130 !> - G: the metric S_pq*delta_IJ
131 !> For the xc functionals, C_aa + C_ab or C_aa - C_ab are stored in the same matrix
132 !> In the above definitions, I,J label the donnor MOs and p,q the sgfs of the basis
133 !>
134 !> -In the SPIN-UNRESTRICTED, spin-conserving case:
135 !> the on- and off-diagonal elements of M are:
136 !> M_d = A + B + C -D
137 !> M_o = B + C - E
138 !> where the submatrices A, B, C, D and E are:
139 !> - A: the groun-state contribution: (F_pq*delta_IJ - epsilon_IJ*S_pq) * delta_ab
140 !> - B: the Coulomb kernel: (pI_a|J_b q)
141 !> - C: the xc kernel: (pI_a|fxc_ab|J_b q)
142 !> - D: the on-diagonal exact-exchange kernel: (pq|I_a J_b) delta_ab
143 !> - E: the off-diagonal exact-exchange kernel: (pJ_b|I_a q) delta_ab
144 !> - G: the metric S_pq*delta_IJ*delta_ab
145 !> p,q label the sgfs, I,J the donro MOs and a,b the spins
146 !>
147 !> -In both above cases, the matrix M is always projected onto the unperturbed unoccupied
148 !> ground state: M <= Q * M * Q^T = (1 - SP) * M * (1 - PS)
149 !>
150 !> -In the SPIN-FLIP case:
151 !> Only the TDA is implemented, that is, there are only on-diagonal elements:
152 !> M_d = A + C - D
153 !> where the submatrices A, C and D are:
154 !> - A: the ground state-contribution: (F_pq*delta_IJ - epsilon_IJ*S_pq) * delta_ab, but here,
155 !> the alph-alpha quadrant has the beta Fock matrix and
156 !> the beta-beta quadrant has the alpha Fock matrix
157 !> - C: the SF xc kernel: (pI_a|fxc|J_bq), fxc = 1/m * (vxc_a -vxc_b)
158 !> - D: the on-diagonal exact-exchange kernel: (pq|I_a J_b) delta_ab
159 !> To ensure that all excitation start from a given spin to the opposite, we then multiply
160 !> by a Q projector where we swap the alpha-alpha and beta-beta spin-quadrants
161 !>
162 !> All possibilities: TDA or full-TDDFT, singlet or triplet, xc or hybrid, etc are treated
163 !> in the same routine to avoid recomputing stuff
164 !> Under TDA, only the on-diagonal elements of M are computed
165 !> In the case of non-TDA, one turns the problem Hermitian
166 ! **************************************************************************************************
167  SUBROUTINE setup_xas_tdp_prob(donor_state, qs_env, xas_tdp_env, xas_tdp_control)
168 
169  TYPE(donor_state_type), POINTER :: donor_state
170  TYPE(qs_environment_type), POINTER :: qs_env
171  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
172  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
173 
174  CHARACTER(len=*), PARAMETER :: routinen = 'setup_xas_tdp_prob'
175 
176  INTEGER :: handle
177  INTEGER, DIMENSION(:), POINTER :: submat_blk_size
178  LOGICAL :: do_coul, do_hfx, do_os, do_sc, do_sf, &
179  do_sg, do_tda, do_tp, do_xc
180  REAL(dp) :: eps_filter, sx
181  TYPE(dbcsr_distribution_type), POINTER :: submat_dist
182  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ex_ker, xc_ker
183  TYPE(dbcsr_type) :: matrix_a, matrix_a_sf, matrix_b, proj_q, &
184  proj_q_sf, work
185  TYPE(dbcsr_type), POINTER :: matrix_c_sc, matrix_c_sf, matrix_c_sg, matrix_c_tp, matrix_d, &
186  matrix_e_sc, sc_matrix_tdp, sf_matrix_tdp, sg_matrix_tdp, tp_matrix_tdp
187 
188  NULLIFY (sg_matrix_tdp, tp_matrix_tdp, submat_dist, submat_blk_size, matrix_c_sf)
189  NULLIFY (matrix_c_sg, matrix_c_tp, matrix_c_sc, matrix_d, matrix_e_sc)
190  NULLIFY (sc_matrix_tdp, sf_matrix_tdp, ex_ker, xc_ker)
191 
192  CALL timeset(routinen, handle)
193 
194 ! Initialization
195  do_os = xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks
196  do_sc = xas_tdp_control%do_spin_cons
197  do_sf = xas_tdp_control%do_spin_flip
198  do_sg = xas_tdp_control%do_singlet
199  do_tp = xas_tdp_control%do_triplet
200  do_xc = xas_tdp_control%do_xc
201  do_hfx = xas_tdp_control%do_hfx
202  do_coul = xas_tdp_control%do_coulomb
203  do_tda = xas_tdp_control%tamm_dancoff
204  sx = xas_tdp_control%sx
205  eps_filter = xas_tdp_control%eps_filter
206  IF (do_sc) THEN
207  ALLOCATE (donor_state%sc_matrix_tdp)
208  sc_matrix_tdp => donor_state%sc_matrix_tdp
209  END IF
210  IF (do_sf) THEN
211  ALLOCATE (donor_state%sf_matrix_tdp)
212  sf_matrix_tdp => donor_state%sf_matrix_tdp
213  END IF
214  IF (do_sg) THEN
215  ALLOCATE (donor_state%sg_matrix_tdp)
216  sg_matrix_tdp => donor_state%sg_matrix_tdp
217  END IF
218  IF (do_tp) THEN
219  ALLOCATE (donor_state%tp_matrix_tdp)
220  tp_matrix_tdp => donor_state%tp_matrix_tdp
221  END IF
222 
223 ! Get the dist and block size of all matrices A, B, C, etc
224  CALL compute_submat_dist_and_blk_size(donor_state, do_os, qs_env)
225  submat_dist => donor_state%dbcsr_dist
226  submat_blk_size => donor_state%blk_size
227 
228 ! Allocate and compute all the matrices A, B, C, etc we will need
229 
230  ! The projector(s) on the unoccupied unperturbed ground state 1-SP and associated work matrix
231  IF (do_sg .OR. do_tp .OR. do_sc) THEN !spin-conserving
232  CALL get_q_projector(proj_q, donor_state, do_os, xas_tdp_env)
233  END IF
234  IF (do_sf) THEN !spin-flip
235  CALL get_q_projector(proj_q_sf, donor_state, do_os, xas_tdp_env, do_sf=.true.)
236  END IF
237  CALL dbcsr_create(matrix=work, matrix_type=dbcsr_type_no_symmetry, dist=submat_dist, &
238  name="WORK", row_blk_size=submat_blk_size, col_blk_size=submat_blk_size)
239 
240  ! The ground state contribution(s)
241  IF (do_sg .OR. do_tp .OR. do_sc) THEN !spin-conserving
242  CALL build_gs_contribution(matrix_a, donor_state, do_os, qs_env)
243  END IF
244  IF (do_sf) THEN !spin-flip
245  CALL build_gs_contribution(matrix_a_sf, donor_state, do_os, qs_env, do_sf=.true.)
246  END IF
247 
248  ! The Coulomb and XC kernels. Internal analysis to know which matrix to compute
249  CALL dbcsr_allocate_matrix_set(xc_ker, 4)
250  ALLOCATE (xc_ker(1)%matrix, xc_ker(2)%matrix, xc_ker(3)%matrix, xc_ker(4)%matrix)
251  CALL kernel_coulomb_xc(matrix_b, xc_ker, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
252  matrix_c_sg => xc_ker(1)%matrix; matrix_c_tp => xc_ker(2)%matrix
253  matrix_c_sc => xc_ker(3)%matrix; matrix_c_sf => xc_ker(4)%matrix
254 
255  ! The exact exchange. Internal analysis to know which matrices to compute
256  CALL dbcsr_allocate_matrix_set(ex_ker, 2)
257  ALLOCATE (ex_ker(1)%matrix, ex_ker(2)%matrix)
258  CALL kernel_exchange(ex_ker, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
259  matrix_d => ex_ker(1)%matrix; matrix_e_sc => ex_ker(2)%matrix
260 
261  ! Build the metric G, also need its inverse in case of full-TDDFT
262  IF (do_tda) THEN
263  ALLOCATE (donor_state%metric(1))
264  CALL build_metric(donor_state%metric, donor_state, qs_env, do_os)
265  ELSE
266  ALLOCATE (donor_state%metric(2))
267  CALL build_metric(donor_state%metric, donor_state, qs_env, do_os, do_inv=.true.)
268  END IF
269 
270 ! Build the eigenvalue problem, depending on the case (TDA, singlet, triplet, hfx, etc ...)
271  IF (do_tda) THEN
272 
273  IF (do_sc) THEN ! open-shell spin-conserving under TDA
274 
275  ! The final matrix is M = A + B + C - D
276  CALL dbcsr_copy(sc_matrix_tdp, matrix_a, name="OS MATRIX TDP")
277  IF (do_coul) CALL dbcsr_add(sc_matrix_tdp, matrix_b, 1.0_dp, 1.0_dp)
278 
279  IF (do_xc) CALL dbcsr_add(sc_matrix_tdp, matrix_c_sc, 1.0_dp, 1.0_dp) !xc kernel
280  IF (do_hfx) CALL dbcsr_add(sc_matrix_tdp, matrix_d, 1.0_dp, -1.0_dp*sx) !scaled hfx
281 
282  ! The product with the Q projector
283  CALL dbcsr_multiply('N', 'N', 1.0_dp, proj_q, sc_matrix_tdp, 0.0_dp, work, filter_eps=eps_filter)
284  CALL dbcsr_multiply('N', 'T', 1.0_dp, work, proj_q, 0.0_dp, sc_matrix_tdp, filter_eps=eps_filter)
285 
286  END IF !do_sc
287 
288  IF (do_sf) THEN ! open-shell spin-flip under TDA
289 
290  ! The final matrix is M = A + C - D
291  CALL dbcsr_copy(sf_matrix_tdp, matrix_a_sf, name="OS MATRIX TDP")
292 
293  IF (do_xc) CALL dbcsr_add(sf_matrix_tdp, matrix_c_sf, 1.0_dp, 1.0_dp) !xc kernel
294  IF (do_hfx) CALL dbcsr_add(sf_matrix_tdp, matrix_d, 1.0_dp, -1.0_dp*sx) !scaled hfx
295 
296  ! Take the product with the (spin-flip) Q projector
297  CALL dbcsr_multiply('N', 'N', 1.0_dp, proj_q_sf, sf_matrix_tdp, 0.0_dp, work, filter_eps=eps_filter)
298  CALL dbcsr_multiply('N', 'T', 1.0_dp, work, proj_q_sf, 0.0_dp, sf_matrix_tdp, filter_eps=eps_filter)
299 
300  END IF !do_sf
301 
302  IF (do_sg) THEN ! singlets under TDA
303 
304  ! The final matrix is M = A + 2B + (C_aa + C_ab) - D
305  CALL dbcsr_copy(sg_matrix_tdp, matrix_a, name="SINGLET MATRIX TDP")
306  IF (do_coul) CALL dbcsr_add(sg_matrix_tdp, matrix_b, 1.0_dp, 2.0_dp)
307 
308  IF (do_xc) CALL dbcsr_add(sg_matrix_tdp, matrix_c_sg, 1.0_dp, 1.0_dp) ! xc kernel
309  IF (do_hfx) CALL dbcsr_add(sg_matrix_tdp, matrix_d, 1.0_dp, -1.0_dp*sx) ! scaled hfx
310 
311  ! Take the product with the Q projector:
312  CALL dbcsr_multiply('N', 'N', 1.0_dp, proj_q, sg_matrix_tdp, 0.0_dp, work, filter_eps=eps_filter)
313  CALL dbcsr_multiply('N', 'T', 1.0_dp, work, proj_q, 0.0_dp, sg_matrix_tdp, filter_eps=eps_filter)
314 
315  END IF !do_sg (TDA)
316 
317  IF (do_tp) THEN ! triplets under TDA
318 
319  ! The final matrix is M = A + (C_aa - C_ab) - D
320  CALL dbcsr_copy(tp_matrix_tdp, matrix_a, name="TRIPLET MATRIX TDP")
321 
322  IF (do_xc) CALL dbcsr_add(tp_matrix_tdp, matrix_c_tp, 1.0_dp, 1.0_dp) ! xc_kernel
323  IF (do_hfx) CALL dbcsr_add(tp_matrix_tdp, matrix_d, 1.0_dp, -1.0_dp*sx) ! scaled hfx
324 
325  ! Take the product with the Q projector:
326  CALL dbcsr_multiply('N', 'N', 1.0_dp, proj_q, tp_matrix_tdp, 0.0_dp, work, filter_eps=eps_filter)
327  CALL dbcsr_multiply('N', 'T', 1.0_dp, work, proj_q, 0.0_dp, tp_matrix_tdp, filter_eps=eps_filter)
328 
329  END IF !do_tp (TDA)
330 
331  ELSE ! not TDA
332 
333  ! In the case of full-TDDFT, the problem is turned Hermitian with the help of auxiliary
334  ! matrices AUX = (A-D+E)^(+-0.5) that are stored in donor_state
335  CALL build_aux_matrix(1.0e-8_dp, sx, matrix_a, matrix_d, matrix_e_sc, do_hfx, proj_q, &
336  work, donor_state, eps_filter, qs_env)
337 
338  IF (do_sc) THEN !full-TDDFT open-shell spin-conserving
339 
340  ! The final matrix is the sum of the on- and off-diagonal elements as in the description
341  ! M = A + 2B + 2C - D - E
342  CALL dbcsr_copy(sc_matrix_tdp, matrix_a, name="OS MATRIX TDP")
343  IF (do_coul) CALL dbcsr_add(sc_matrix_tdp, matrix_b, 1.0_dp, 2.0_dp)
344 
345  IF (do_hfx) THEN !scaled hfx
346  CALL dbcsr_add(sc_matrix_tdp, matrix_d, 1.0_dp, -1.0_dp*sx)
347  CALL dbcsr_add(sc_matrix_tdp, matrix_e_sc, 1.0_dp, -1.0_dp*sx)
348  END IF
349  IF (do_xc) THEN
350  CALL dbcsr_add(sc_matrix_tdp, matrix_c_sc, 1.0_dp, 2.0_dp)
351  END IF
352 
353  ! Take the product with the Q projector
354  CALL dbcsr_multiply('N', 'N', 1.0_dp, proj_q, sc_matrix_tdp, 0.0_dp, work, filter_eps=eps_filter)
355  CALL dbcsr_multiply('N', 'T', 1.0_dp, work, proj_q, 0.0_dp, sc_matrix_tdp, filter_eps=eps_filter)
356 
357  ! Take the product with the inverse metric
358  ! M <= G^-1 * M * G^-1
359  CALL dbcsr_multiply('N', 'N', 1.0_dp, donor_state%metric(2)%matrix, sc_matrix_tdp, &
360  0.0_dp, work, filter_eps=eps_filter)
361  CALL dbcsr_multiply('N', 'N', 1.0_dp, work, donor_state%metric(2)%matrix, 0.0_dp, &
362  sc_matrix_tdp, filter_eps=eps_filter)
363 
364  END IF
365 
366  IF (do_sg) THEN ! full-TDDFT singlets
367 
368  ! The final matrix is the sum of the on- and off-diagonal elements as in the description
369  ! M = A + 4B + 2(C_aa + C_ab) - D - E
370  CALL dbcsr_copy(sg_matrix_tdp, matrix_a, name="SINGLET MATRIX TDP")
371  IF (do_coul) CALL dbcsr_add(sg_matrix_tdp, matrix_b, 1.0_dp, 4.0_dp)
372 
373  IF (do_hfx) THEN !scaled hfx
374  CALL dbcsr_add(sg_matrix_tdp, matrix_d, 1.0_dp, -1.0_dp*sx)
375  CALL dbcsr_add(sg_matrix_tdp, matrix_e_sc, 1.0_dp, -1.0_dp*sx)
376  END IF
377  IF (do_xc) THEN !xc kernel
378  CALL dbcsr_add(sg_matrix_tdp, matrix_c_sg, 1.0_dp, 2.0_dp)
379  END IF
380 
381  ! Take the product with the Q projector
382  CALL dbcsr_multiply('N', 'N', 1.0_dp, proj_q, sg_matrix_tdp, 0.0_dp, work, filter_eps=eps_filter)
383  CALL dbcsr_multiply('N', 'T', 1.0_dp, work, proj_q, 0.0_dp, sg_matrix_tdp, filter_eps=eps_filter)
384 
385  ! Take the product with the inverse metric
386  ! M <= G^-1 * M * G^-1
387  CALL dbcsr_multiply('N', 'N', 1.0_dp, donor_state%metric(2)%matrix, sg_matrix_tdp, &
388  0.0_dp, work, filter_eps=eps_filter)
389  CALL dbcsr_multiply('N', 'N', 1.0_dp, work, donor_state%metric(2)%matrix, 0.0_dp, &
390  sg_matrix_tdp, filter_eps=eps_filter)
391 
392  END IF ! singlets
393 
394  IF (do_tp) THEN ! full-TDDFT triplets
395 
396  ! The final matrix is the sum of the on- and off-diagonal elements as in the description
397  ! M = A + 2(C_aa - C_ab) - D - E
398  CALL dbcsr_copy(tp_matrix_tdp, matrix_a, name="TRIPLET MATRIX TDP")
399 
400  IF (do_hfx) THEN !scaled hfx
401  CALL dbcsr_add(tp_matrix_tdp, matrix_d, 1.0_dp, -1.0_dp*sx)
402  CALL dbcsr_add(tp_matrix_tdp, matrix_e_sc, 1.0_dp, -1.0_dp*sx)
403  END IF
404  IF (do_xc) THEN
405  CALL dbcsr_add(tp_matrix_tdp, matrix_c_tp, 1.0_dp, 2.0_dp)
406  END IF
407 
408  ! Take the product with the Q projector
409  CALL dbcsr_multiply('N', 'N', 1.0_dp, proj_q, tp_matrix_tdp, 0.0_dp, work, filter_eps=eps_filter)
410  CALL dbcsr_multiply('N', 'T', 1.0_dp, work, proj_q, 0.0_dp, tp_matrix_tdp, filter_eps=eps_filter)
411 
412  ! Take the product with the inverse metric
413  ! M <= G^-1 * M * G^-1
414  CALL dbcsr_multiply('N', 'N', 1.0_dp, donor_state%metric(2)%matrix, tp_matrix_tdp, &
415  0.0_dp, work, filter_eps=eps_filter)
416  CALL dbcsr_multiply('N', 'N', 1.0_dp, work, donor_state%metric(2)%matrix, 0.0_dp, &
417  tp_matrix_tdp, filter_eps=eps_filter)
418 
419  END IF ! triplets
420 
421  END IF ! test on TDA
422 
423 ! Clean-up
424  CALL dbcsr_release(matrix_a)
425  CALL dbcsr_release(matrix_a_sf)
426  CALL dbcsr_release(matrix_b)
427  CALL dbcsr_release(proj_q)
428  CALL dbcsr_release(proj_q_sf)
429  CALL dbcsr_release(work)
430  CALL dbcsr_deallocate_matrix_set(ex_ker)
431  CALL dbcsr_deallocate_matrix_set(xc_ker)
432 
433  CALL timestop(handle)
434 
435  END SUBROUTINE setup_xas_tdp_prob
436 
437 ! **************************************************************************************************
438 !> \brief Solves the XAS TDP generalized eigenvalue problem omega*C = matrix_tdp*C using standard
439 !> full diagonalization methods. The problem is Hermitian (made that way even if not TDA)
440 !> \param donor_state ...
441 !> \param xas_tdp_control ...
442 !> \param xas_tdp_env ...
443 !> \param qs_env ...
444 !> \param ex_type whether we deal with singlets, triplets, spin-conserving open-shell or spin-flip
445 !> \note The computed eigenvalues and eigenvectors are stored in the donor_state
446 !> The eigenvectors are the LR-coefficients. In case of TDA, c^- is stored. In the general
447 !> case, the sum c^+ + c^- is stored.
448 !> - Spin-restricted:
449 !> In case both singlets and triplets are considered, this routine must be called twice. This
450 !> is the choice that was made because the body of the routine is exactly the same in both cases
451 !> Note that for singlet we solve for u = 1/sqrt(2)*(c_alpha + c_beta) = sqrt(2)*c
452 !> and that for triplets we solve for v = 1/sqrt(2)*(c_alpha - c_beta) = sqrt(2)*c
453 !> - Spin-unrestricted:
454 !> The problem is solved for the LR coefficients c_pIa as they are (not linear combination)
455 !> The routine might be called twice (once for spin-conservign, one for spin-flip)
456 ! **************************************************************************************************
457  SUBROUTINE solve_xas_tdp_prob(donor_state, xas_tdp_control, xas_tdp_env, qs_env, ex_type)
458 
459  TYPE(donor_state_type), POINTER :: donor_state
460  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
461  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
462  TYPE(qs_environment_type), POINTER :: qs_env
463  INTEGER, INTENT(IN) :: ex_type
464 
465  CHARACTER(len=*), PARAMETER :: routinen = 'solve_xas_tdp_prob'
466 
467  INTEGER :: first_ex, handle, i, imo, ispin, nao, &
468  ndo_mo, nelectron, nevals, nocc, nrow, &
469  nspins, ot_nevals
470  LOGICAL :: do_os, do_range, do_sf
471  REAL(dp) :: ot_elb
472  REAL(dp), ALLOCATABLE, DIMENSION(:) :: scaling, tmp_evals
473  REAL(dp), DIMENSION(:), POINTER :: lr_evals
474  TYPE(cp_blacs_env_type), POINTER :: blacs_env
475  TYPE(cp_fm_struct_type), POINTER :: ex_struct, fm_struct, ot_fm_struct
476  TYPE(cp_fm_type) :: c_diff, c_sum, lhs_matrix, rhs_matrix, &
477  work
478  TYPE(cp_fm_type), POINTER :: lr_coeffs
479  TYPE(dbcsr_type) :: tmp_mat, tmp_mat2
480  TYPE(dbcsr_type), POINTER :: matrix_tdp
481  TYPE(mp_para_env_type), POINTER :: para_env
482 
483  CALL timeset(routinen, handle)
484 
485  NULLIFY (para_env, blacs_env, fm_struct, matrix_tdp)
486  NULLIFY (ex_struct, lr_evals, lr_coeffs)
487  cpassert(ASSOCIATED(xas_tdp_env))
488 
489  do_os = .false.
490  do_sf = .false.
491  IF (ex_type == tddfpt_spin_cons) THEN
492  matrix_tdp => donor_state%sc_matrix_tdp
493  do_os = .true.
494  ELSE IF (ex_type == tddfpt_spin_flip) THEN
495  matrix_tdp => donor_state%sf_matrix_tdp
496  do_os = .true.
497  do_sf = .true.
498  ELSE IF (ex_type == tddfpt_singlet) THEN
499  matrix_tdp => donor_state%sg_matrix_tdp
500  ELSE IF (ex_type == tddfpt_triplet) THEN
501  matrix_tdp => donor_state%tp_matrix_tdp
502  END IF
503  CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env, nelectron_total=nelectron)
504 
505 ! Initialization
506  nspins = 1; IF (do_os) nspins = 2
507  CALL cp_fm_get_info(donor_state%gs_coeffs, nrow_global=nao)
508  CALL dbcsr_get_info(matrix_tdp, nfullrows_total=nrow)
509  ndo_mo = donor_state%ndo_mo
510  nocc = nelectron/2; IF (do_os) nocc = nelectron
511  nocc = ndo_mo*nocc
512 
513  !solve by energy_range or number of states ?
514  do_range = .false.
515  IF (xas_tdp_control%e_range > 0.0_dp) do_range = .true.
516 
517  ! create the fm infrastructure
518  CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nrow, &
519  para_env=para_env, ncol_global=nrow)
520  CALL cp_fm_create(rhs_matrix, fm_struct)
521  CALL cp_fm_create(work, fm_struct)
522 
523 ! Test on TDA
524  IF (xas_tdp_control%tamm_dancoff) THEN
525 
526  IF (xas_tdp_control%do_ot) THEN
527 
528  !need to precompute the number of evals for OT
529  IF (do_range) THEN
530 
531  !in case of energy range criterion, use LUMO eigenvalues as estimate
532  ot_elb = xas_tdp_env%lumo_evals(1)%array(1)
533  IF (do_os) ot_elb = min(ot_elb, xas_tdp_env%lumo_evals(2)%array(1))
534 
535  ot_nevals = count(xas_tdp_env%lumo_evals(1)%array - ot_elb .LE. xas_tdp_control%e_range)
536  IF (do_os) ot_nevals = ot_nevals + &
537  count(xas_tdp_env%lumo_evals(2)%array - ot_elb .LE. xas_tdp_control%e_range)
538 
539  ELSE
540 
541  ot_nevals = nspins*nao - nocc/ndo_mo
542  IF (xas_tdp_control%n_excited > 0 .AND. xas_tdp_control%n_excited < ot_nevals) THEN
543  ot_nevals = xas_tdp_control%n_excited
544  END IF
545  END IF
546  ot_nevals = ndo_mo*ot_nevals !as in input description, multiply by multiplicity of donor state
547 
548 ! Organize results data
549  first_ex = 1
550  ALLOCATE (tmp_evals(ot_nevals))
551  CALL cp_fm_struct_create(ot_fm_struct, context=blacs_env, para_env=para_env, &
552  nrow_global=nrow, ncol_global=ot_nevals)
553  CALL cp_fm_create(c_sum, ot_fm_struct)
554 
555  CALL xas_ot_solver(matrix_tdp, donor_state%metric(1)%matrix, c_sum, tmp_evals, ot_nevals, &
556  do_sf, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
557 
558  CALL cp_fm_struct_release(ot_fm_struct)
559 
560  ELSE
561 
562 ! Organize results data
563  first_ex = nocc + 1 !where to find the first proper eigenvalue
564  ALLOCATE (tmp_evals(nrow))
565  CALL cp_fm_create(c_sum, fm_struct)
566 
567 ! Get the main matrix_tdp as an fm
568  CALL copy_dbcsr_to_fm(matrix_tdp, rhs_matrix)
569 
570 ! Get the metric as a fm
571  CALL cp_fm_create(lhs_matrix, fm_struct)
572  CALL copy_dbcsr_to_fm(donor_state%metric(1)%matrix, lhs_matrix)
573 
574  !Diagonalisation (Cholesky decomposition). In TDA, c_sum = c^-
575  CALL cp_fm_geeig(rhs_matrix, lhs_matrix, c_sum, tmp_evals, work)
576 
577 ! TDA specific clean-up
578  CALL cp_fm_release(lhs_matrix)
579 
580  END IF
581 
582  ELSE ! not TDA
583 
584 ! Organize results data
585  first_ex = nocc + 1
586  ALLOCATE (tmp_evals(nrow))
587  CALL cp_fm_create(c_sum, fm_struct)
588 
589 ! Need to multiply the current matrix_tdp with the auxiliary matrix
590 ! tmp_mat = (A-D+E)^0.5 * M * (A-D+E)^0.5
591  CALL dbcsr_create(matrix=tmp_mat, template=matrix_tdp, matrix_type=dbcsr_type_no_symmetry)
592  CALL dbcsr_create(matrix=tmp_mat2, template=matrix_tdp, matrix_type=dbcsr_type_no_symmetry)
593  CALL dbcsr_multiply('N', 'N', 1.0_dp, donor_state%matrix_aux, matrix_tdp, &
594  0.0_dp, tmp_mat2, filter_eps=xas_tdp_control%eps_filter)
595  CALL dbcsr_multiply('N', 'N', 1.0_dp, tmp_mat2, donor_state%matrix_aux, &
596  0.0_dp, tmp_mat, filter_eps=xas_tdp_control%eps_filter)
597 
598 ! Get the matrix as a fm
599  CALL copy_dbcsr_to_fm(tmp_mat, rhs_matrix)
600 
601 ! Solve the "turned-Hermitian" eigenvalue problem
602  CALL choose_eigv_solver(rhs_matrix, work, tmp_evals)
603 
604 ! Currently, work = (A-D+E)^0.5 (c^+ - c^-) and tmp_evals = omega^2
605 ! Put tiny almost zero eigenvalues to zero (corresponding to occupied MOs)
606  WHERE (tmp_evals < 1.0e-4_dp) tmp_evals = 0.0_dp
607 
608 ! Retrieve c_diff = (c^+ - c^-) for normalization
609 ! (c^+ - c^-) = 1/omega^2 * M * (A-D+E)^0.5 * work
610  CALL cp_fm_create(c_diff, fm_struct)
611  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_tdp, donor_state%matrix_aux, &
612  0.0_dp, tmp_mat, filter_eps=xas_tdp_control%eps_filter)
613  CALL cp_dbcsr_sm_fm_multiply(tmp_mat, work, c_diff, ncol=nrow)
614 
615  ALLOCATE (scaling(nrow))
616  scaling = 0.0_dp
617  WHERE (abs(tmp_evals) > 1.0e-8_dp) scaling = 1.0_dp/tmp_evals
618  CALL cp_fm_column_scale(c_diff, scaling)
619 
620 ! Normalize with the metric: c_diff * G * c_diff = +- 1
621  scaling = 0.0_dp
622  CALL get_normal_scaling(scaling, c_diff, donor_state)
623  CALL cp_fm_column_scale(c_diff, scaling)
624 
625 ! Get the actual eigenvalues
626  tmp_evals = sqrt(tmp_evals)
627 
628 ! Get c_sum = (c^+ + c^-), which appears in all transition density related expressions
629 ! c_sum = -1/omega G^-1 * (A-D+E) * (c^+ - c^-)
630  CALL dbcsr_multiply('N', 'N', 1.0_dp, donor_state%matrix_aux, donor_state%matrix_aux, &
631  0.0_dp, tmp_mat2, filter_eps=xas_tdp_control%eps_filter)
632  CALL dbcsr_multiply('N', 'N', 1.0_dp, donor_state%metric(2)%matrix, tmp_mat2, &
633  0.0_dp, tmp_mat, filter_eps=xas_tdp_control%eps_filter)
634  CALL cp_dbcsr_sm_fm_multiply(tmp_mat, c_diff, c_sum, ncol=nrow)
635  WHERE (tmp_evals .NE. 0) scaling = -1.0_dp/tmp_evals
636  CALL cp_fm_column_scale(c_sum, scaling)
637 
638 ! Full TDDFT specific clean-up
639  CALL cp_fm_release(c_diff)
640  CALL dbcsr_release(tmp_mat)
641  CALL dbcsr_release(tmp_mat2)
642  DEALLOCATE (scaling)
643 
644  END IF ! TDA
645 
646 ! Full matrix clean-up
647  CALL cp_fm_release(rhs_matrix)
648  CALL cp_fm_release(work)
649 
650 ! Reorganize the eigenvalues, we want a lr_evals array with the proper dimension and where the
651 ! first element is the first eval. Need a case study on do_range/ot
652  IF (xas_tdp_control%do_ot) THEN
653 
654  nevals = ot_nevals
655 
656  ELSE IF (do_range) THEN
657 
658  WHERE (tmp_evals > tmp_evals(first_ex) + xas_tdp_control%e_range) tmp_evals = 0.0_dp
659  nevals = maxloc(tmp_evals, 1) - nocc
660 
661  ELSE
662 
663  !Determine the number of evals to keep base on N_EXCITED
664  nevals = nspins*nao - nocc/ndo_mo
665  IF (xas_tdp_control%n_excited > 0 .AND. xas_tdp_control%n_excited < nevals) THEN
666  nevals = xas_tdp_control%n_excited
667  END IF
668  nevals = ndo_mo*nevals !as in input description, multiply by # of donor MOs
669 
670  END IF
671 
672  ALLOCATE (lr_evals(nevals))
673  lr_evals(:) = tmp_evals(first_ex:first_ex + nevals - 1)
674 
675 ! Reorganize the eigenvectors in array of cp_fm so that each ndo_mo columns corresponds to an
676 ! excited state. Makes later calls to those easier and more efficient
677 ! In case of open-shell, we store the coeffs in the same logic as the matrix => first block where
678 ! the columns are the c_Ialpha and second block with columns as c_Ibeta
679  CALL cp_fm_struct_create(ex_struct, nrow_global=nao, ncol_global=ndo_mo*nspins*nevals, &
680  para_env=para_env, context=blacs_env)
681  ALLOCATE (lr_coeffs)
682  CALL cp_fm_create(lr_coeffs, ex_struct)
683 
684  DO i = 1, nevals
685  DO ispin = 1, nspins
686  DO imo = 1, ndo_mo
687 
688  CALL cp_fm_to_fm_submat(msource=c_sum, mtarget=lr_coeffs, &
689  nrow=nao, ncol=1, s_firstrow=((ispin - 1)*ndo_mo + imo - 1)*nao + 1, &
690  s_firstcol=first_ex + i - 1, t_firstrow=1, &
691  t_firstcol=(i - 1)*ndo_mo*nspins + (ispin - 1)*ndo_mo + imo)
692  END DO !imo
693  END DO !ispin
694  END DO !istate
695 
696  IF (ex_type == tddfpt_spin_cons) THEN
697  donor_state%sc_coeffs => lr_coeffs
698  donor_state%sc_evals => lr_evals
699  ELSE IF (ex_type == tddfpt_spin_flip) THEN
700  donor_state%sf_coeffs => lr_coeffs
701  donor_state%sf_evals => lr_evals
702  ELSE IF (ex_type == tddfpt_singlet) THEN
703  donor_state%sg_coeffs => lr_coeffs
704  donor_state%sg_evals => lr_evals
705  ELSE IF (ex_type == tddfpt_triplet) THEN
706  donor_state%tp_coeffs => lr_coeffs
707  donor_state%tp_evals => lr_evals
708  END IF
709 
710 ! Clean-up
711  CALL cp_fm_release(c_sum)
712  CALL cp_fm_struct_release(fm_struct)
713  CALL cp_fm_struct_release(ex_struct)
714 
715 ! Perform a partial clean-up of the donor_state
716  CALL dbcsr_release(matrix_tdp)
717 
718  CALL timestop(handle)
719 
720  END SUBROUTINE solve_xas_tdp_prob
721 
722 ! **************************************************************************************************
723 !> \brief An iterative solver based on OT for the TDA generalized eigV problem lambda Sx = Hx
724 !> \param matrix_tdp the RHS matrix (dbcsr)
725 !> \param metric the LHS matrix (dbcsr)
726 !> \param evecs the corresponding eigenvectors (fm)
727 !> \param evals the corresponding eigenvalues
728 !> \param neig the number of wanted eigenvalues
729 !> \param do_sf whther spin-flip TDDFT is on
730 !> \param donor_state ...
731 !> \param xas_tdp_env ...
732 !> \param xas_tdp_control ...
733 !> \param qs_env ...
734 ! **************************************************************************************************
735  SUBROUTINE xas_ot_solver(matrix_tdp, metric, evecs, evals, neig, do_sf, donor_state, xas_tdp_env, &
736  xas_tdp_control, qs_env)
737 
738  TYPE(dbcsr_type), POINTER :: matrix_tdp, metric
739  TYPE(cp_fm_type), INTENT(IN) :: evecs
740  REAL(dp), DIMENSION(:) :: evals
741  INTEGER, INTENT(IN) :: neig
742  LOGICAL :: do_sf
743  TYPE(donor_state_type), POINTER :: donor_state
744  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
745  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
746  TYPE(qs_environment_type), POINTER :: qs_env
747 
748  CHARACTER(len=*), PARAMETER :: routinen = 'xas_ot_solver'
749 
750  INTEGER :: handle, max_iter, ndo_mo, nelec_spin(2), &
751  nocc, nrow, output_unit
752  LOGICAL :: do_os
753  REAL(dp) :: eps_iter
754  TYPE(cp_blacs_env_type), POINTER :: blacs_env
755  TYPE(cp_fm_struct_type), POINTER :: ortho_struct
756  TYPE(cp_fm_type) :: ortho_space
757  TYPE(dbcsr_type), POINTER :: ot_prec
758  TYPE(mp_para_env_type), POINTER :: para_env
759  TYPE(preconditioner_type), POINTER :: precond
760 
761  NULLIFY (para_env, blacs_env, ortho_struct, ot_prec)
762 
763  CALL timeset(routinen, handle)
764 
765  output_unit = cp_logger_get_default_io_unit()
766  IF (output_unit > 0) THEN
767  WRITE (output_unit, "(/,T5,A)") &
768  "Using OT eigensolver for diagonalization: "
769  END IF
770 
771  do_os = xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks
772  ndo_mo = donor_state%ndo_mo
773  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env, nelectron_spin=nelec_spin)
774  CALL cp_fm_get_info(evecs, nrow_global=nrow)
775  max_iter = xas_tdp_control%ot_max_iter
776  eps_iter = xas_tdp_control%ot_eps_iter
777  nocc = nelec_spin(1)/2*ndo_mo
778  IF (do_os) nocc = sum(nelec_spin)*ndo_mo
779 
780 ! Initialize relevent matrices
781  ALLOCATE (ot_prec)
782  CALL dbcsr_create(ot_prec, template=matrix_tdp)
783  CALL cp_fm_struct_create(ortho_struct, context=blacs_env, para_env=para_env, &
784  nrow_global=nrow, ncol_global=nocc)
785  CALL cp_fm_create(ortho_space, ortho_struct)
786 
787  CALL prep_for_ot(evecs, ortho_space, ot_prec, neig, do_sf, donor_state, xas_tdp_env, &
788  xas_tdp_control, qs_env)
789 
790 ! Prepare the preconditioner
791  ALLOCATE (precond)
792  CALL init_preconditioner(precond, para_env, blacs_env)
793  precond%in_use = ot_precond_full_single ! because applying this conditioner is only a mm
794  precond%dbcsr_matrix => ot_prec
795 
796 ! Actually solving the eigenvalue problem
797  CALL ot_eigensolver(matrix_h=matrix_tdp, matrix_s=metric, matrix_c_fm=evecs, &
798  eps_gradient=eps_iter, iter_max=max_iter, silent=.false., &
799  ot_settings=xas_tdp_control%ot_settings, &
800  matrix_orthogonal_space_fm=ortho_space, &
801  preconditioner=precond)
802  CALL calculate_subspace_eigenvalues(evecs, matrix_tdp, evals_arg=evals)
803 
804 ! Clean-up
805  CALL cp_fm_struct_release(ortho_struct)
806  CALL cp_fm_release(ortho_space)
807  CALL dbcsr_release(ot_prec)
808  CALL destroy_preconditioner(precond)
809  DEALLOCATE (precond)
810 
811  CALL timestop(handle)
812 
813  END SUBROUTINE xas_ot_solver
814 
815 ! **************************************************************************************************
816 !> \brief Prepares all required matrices for the OT eigensolver (precond, ortho space and guesses)
817 !> \param guess the guess eigenvectors absed on LUMOs, in fm format
818 !> \param ortho the orthogonal space in fm format (occupied MOs)
819 !> \param precond the OT preconditioner in DBCSR format
820 !> \param neig ...
821 !> \param do_sf ...
822 !> \param donor_state ...
823 !> \param xas_tdp_env ...
824 !> \param xas_tdp_control ...
825 !> \param qs_env ...
826 !> \note Matrices are allocate before entry
827 ! **************************************************************************************************
828  SUBROUTINE prep_for_ot(guess, ortho, precond, neig, do_sf, donor_state, xas_tdp_env, &
829  xas_tdp_control, qs_env)
830 
831  TYPE(cp_fm_type), INTENT(IN) :: guess, ortho
832  TYPE(dbcsr_type) :: precond
833  INTEGER :: neig
834  LOGICAL :: do_sf
835  TYPE(donor_state_type), POINTER :: donor_state
836  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
837  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
838  TYPE(qs_environment_type), POINTER :: qs_env
839 
840  CHARACTER(len=*), PARAMETER :: routinen = 'prep_for_ot'
841 
842  INTEGER :: blk, handle, i, iblk, ido_mo, ispin, jblk, maxel, minel, nao, natom, ndo_mo, &
843  nelec_spin(2), nhomo(2), nlumo(2), nspins, start_block, start_col, start_row
844  LOGICAL :: do_os, found
845  REAL(dp), DIMENSION(:, :), POINTER :: pblock
846  TYPE(cp_fm_type), POINTER :: mo_coeff
847  TYPE(dbcsr_iterator_type) :: iter
848  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
849 
850  NULLIFY (mos, mo_coeff, pblock)
851 
852  !REMINDER on the organization of the xas_tdp matrix. It is DBCSR format, with a super bock
853  !structure. First block structure is spin quadrants: upper left is alpha-alpha spin and lower
854  !right is beta-beta spin. Then each quadrants is divided in a ndo_mo x ndo_mo grid (1x1 for 1s,
855  !2s, 3x3 for 2p). Each block in this grid has the normal DBCSR structure and dist, simply
856  !replicated. The resulting eigenvectors follow the same logic.
857 
858  CALL timeset(routinen, handle)
859 
860  do_os = xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks
861  nspins = 1; IF (do_os) nspins = 2
862  ndo_mo = donor_state%ndo_mo
863  CALL cp_fm_get_info(xas_tdp_env%lumo_evecs(1), nrow_global=nao)
864  CALL get_qs_env(qs_env, natom=natom, nelectron_spin=nelec_spin)
865 
866  !Compute the number of guesses for each spins
867  IF (do_os) THEN
868  minel = minloc(nelec_spin, 1)
869  maxel = 3 - minel
870  nlumo(minel) = (neig/ndo_mo + nelec_spin(maxel) - nelec_spin(minel))/2
871  nlumo(maxel) = neig/ndo_mo - nlumo(minel)
872  ELSE
873  nlumo(1) = neig/ndo_mo
874  END IF
875 
876  !Building the guess vectors based on the LUMOs. Copy LUMOs into approriate spin/do_mo
877  !quadrant/block. Order within a block does not matter
878  !Note: in spin-flip, the upper left quadrant is for beta-alpha transition, so guess are alpha LUMOs
879  start_row = 0
880  start_col = 0
881  DO ispin = 1, nspins
882  DO ido_mo = 1, ndo_mo
883 
884  CALL cp_fm_to_fm_submat(msource=xas_tdp_env%lumo_evecs(ispin), mtarget=guess, &
885  nrow=nao, ncol=nlumo(ispin), s_firstrow=1, s_firstcol=1, &
886  t_firstrow=start_row + 1, t_firstcol=start_col + 1)
887 
888  start_row = start_row + nao
889  start_col = start_col + nlumo(ispin)
890 
891  END DO
892  END DO
893 
894  !Build the orthogonal space according to the same principles, but based on occupied MOs
895  !Note: in spin-flip, the upper left quadrant is for beta-alpha transition, so ortho space is beta HOMOs
896  CALL get_qs_env(qs_env, mos=mos)
897  nhomo = 0
898  DO ispin = 1, nspins
899  CALL get_mo_set(mos(ispin), homo=nhomo(ispin))
900  END DO
901 
902  start_row = 0
903  start_col = 0
904  DO i = 1, nspins
905  ispin = i; IF (do_sf) ispin = 3 - i
906  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
907 
908  DO ido_mo = 1, ndo_mo
909 
910  CALL cp_fm_to_fm_submat(msource=mo_coeff, mtarget=ortho, nrow=nao, ncol=nhomo(ispin), &
911  s_firstrow=1, s_firstcol=1, &
912  t_firstrow=start_row + 1, t_firstcol=start_col + 1)
913 
914  start_row = start_row + nao
915  start_col = start_col + nhomo(ispin)
916 
917  END DO
918  END DO
919 
920  !Build the preconditioner. Copy the "canonical" pre-computed matrix into the proper spin/do_mo
921  !quadrants/blocks. The end matrix is purely block diagonal
922  DO ispin = 1, nspins
923 
924  CALL dbcsr_iterator_start(iter, xas_tdp_env%ot_prec(ispin)%matrix)
925  DO WHILE (dbcsr_iterator_blocks_left(iter))
926 
927  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
928 
929  CALL dbcsr_get_block_p(xas_tdp_env%ot_prec(ispin)%matrix, iblk, jblk, pblock, found)
930 
931  IF (found) THEN
932 
933  start_block = (ispin - 1)*ndo_mo*natom
934  DO ido_mo = 1, ndo_mo
935  CALL dbcsr_put_block(precond, start_block + iblk, start_block + jblk, pblock)
936 
937  start_block = start_block + natom
938 
939  END DO
940  END IF
941 
942  END DO !dbcsr iter
943  CALL dbcsr_iterator_stop(iter)
944  END DO
945 
946  CALL dbcsr_finalize(precond)
947 
948  CALL timestop(handle)
949 
950  END SUBROUTINE prep_for_ot
951 
952 ! **************************************************************************************************
953 !> \brief Returns the scaling to apply to normalize the LR eigenvectors.
954 !> \param scaling the scaling array to apply
955 !> \param lr_coeffs the linear response coefficients as a fm
956 !> \param donor_state ...
957 !> \note The LR coeffs are normalized when c^T G c = +- 1, G is the metric, c = c^- for TDA and
958 !> c = c^+ - c^- for the full problem
959 ! **************************************************************************************************
960  SUBROUTINE get_normal_scaling(scaling, lr_coeffs, donor_state)
961 
962  REAL(dp), ALLOCATABLE, DIMENSION(:) :: scaling
963  TYPE(cp_fm_type), INTENT(IN) :: lr_coeffs
964  TYPE(donor_state_type), POINTER :: donor_state
965 
966  INTEGER :: nrow, nscal, nvals
967  REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag
968  TYPE(cp_blacs_env_type), POINTER :: blacs_env
969  TYPE(cp_fm_struct_type), POINTER :: norm_struct, work_struct
970  TYPE(cp_fm_type) :: fm_norm, work
971  TYPE(mp_para_env_type), POINTER :: para_env
972 
973  NULLIFY (para_env, blacs_env, norm_struct, work_struct)
974 
975 ! Creating the matrix structures and initializing the work matrices
976  CALL cp_fm_get_info(lr_coeffs, context=blacs_env, para_env=para_env, &
977  matrix_struct=work_struct, ncol_global=nvals, nrow_global=nrow)
978  CALL cp_fm_struct_create(norm_struct, para_env=para_env, context=blacs_env, &
979  nrow_global=nvals, ncol_global=nvals)
980 
981  CALL cp_fm_create(work, work_struct)
982  CALL cp_fm_create(fm_norm, norm_struct)
983 
984 ! Taking c^T * G * C
985  CALL cp_dbcsr_sm_fm_multiply(donor_state%metric(1)%matrix, lr_coeffs, work, ncol=nvals)
986  CALL parallel_gemm('T', 'N', nvals, nvals, nrow, 1.0_dp, lr_coeffs, work, 0.0_dp, fm_norm)
987 
988 ! Computing the needed scaling
989  ALLOCATE (diag(nvals))
990  CALL cp_fm_get_diag(fm_norm, diag)
991  WHERE (abs(diag) > 1.0e-8_dp) diag = 1.0_dp/sqrt(abs(diag))
992 
993  nscal = SIZE(scaling)
994  scaling(1:nscal) = diag(1:nscal)
995 
996 ! Clean-up
997  CALL cp_fm_release(work)
998  CALL cp_fm_release(fm_norm)
999  CALL cp_fm_struct_release(norm_struct)
1000 
1001  END SUBROUTINE get_normal_scaling
1002 
1003 ! **************************************************************************************************
1004 !> \brief This subroutine computes the row/column block structure as well as the dbcsr ditrinution
1005 !> for the submatrices making up the generalized XAS TDP eigenvalue problem. They all share
1006 !> the same properties, which are based on the replication of the KS matrix. Stored in the
1007 !> donor_state_type
1008 !> \param donor_state ...
1009 !> \param do_os whether this is a open-shell calculation
1010 !> \param qs_env ...
1011 ! **************************************************************************************************
1012  SUBROUTINE compute_submat_dist_and_blk_size(donor_state, do_os, qs_env)
1013 
1014  TYPE(donor_state_type), POINTER :: donor_state
1015  LOGICAL, INTENT(IN) :: do_os
1016  TYPE(qs_environment_type), POINTER :: qs_env
1017 
1018  INTEGER :: group, i, nao, nblk_row, ndo_mo, nspins, &
1019  scol_dist, srow_dist
1020  INTEGER, DIMENSION(:), POINTER :: col_dist, col_dist_sub, row_blk_size, &
1021  row_dist, row_dist_sub, submat_blk_size
1022  INTEGER, DIMENSION(:, :), POINTER :: pgrid
1023  TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist, submat_dist
1024  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1025 
1026  NULLIFY (matrix_ks, dbcsr_dist, row_blk_size, row_dist, col_dist, pgrid, col_dist_sub)
1027  NULLIFY (row_dist_sub, submat_dist, submat_blk_size)
1028 
1029 ! The submatrices are indexed by M_{pi sig,qj tau}, where p,q label basis functions and i,j donor
1030 ! MOs and sig,tau the spins. For each spin and donor MOs combination, one has a submatrix of the
1031 ! size of the KS matrix (nao x nao) with the same dbcsr block structure
1032 
1033 ! Initialization
1034  ndo_mo = donor_state%ndo_mo
1035  CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks, dbcsr_dist=dbcsr_dist)
1036  CALL dbcsr_get_info(matrix_ks(1)%matrix, row_blk_size=row_blk_size)
1037  CALL dbcsr_distribution_get(dbcsr_dist, row_dist=row_dist, col_dist=col_dist, group=group, &
1038  pgrid=pgrid)
1039  nao = sum(row_blk_size)
1040  nblk_row = SIZE(row_blk_size)
1041  srow_dist = SIZE(row_dist)
1042  scol_dist = SIZE(col_dist)
1043  nspins = 1; IF (do_os) nspins = 2
1044 
1045 ! Creation if submatrix block size and col/row distribution
1046  ALLOCATE (submat_blk_size(ndo_mo*nspins*nblk_row))
1047  ALLOCATE (row_dist_sub(ndo_mo*nspins*srow_dist))
1048  ALLOCATE (col_dist_sub(ndo_mo*nspins*scol_dist))
1049 
1050  DO i = 1, ndo_mo*nspins
1051  submat_blk_size((i - 1)*nblk_row + 1:i*nblk_row) = row_blk_size
1052  row_dist_sub((i - 1)*srow_dist + 1:i*srow_dist) = row_dist
1053  col_dist_sub((i - 1)*scol_dist + 1:i*scol_dist) = col_dist
1054  END DO
1055 
1056 ! Create the submatrix dbcsr distribution
1057  ALLOCATE (submat_dist)
1058  CALL dbcsr_distribution_new(submat_dist, group=group, pgrid=pgrid, row_dist=row_dist_sub, &
1059  col_dist=col_dist_sub)
1060 
1061  donor_state%dbcsr_dist => submat_dist
1062  donor_state%blk_size => submat_blk_size
1063 
1064 ! Clean-up
1065  DEALLOCATE (col_dist_sub, row_dist_sub)
1066 
1067  END SUBROUTINE compute_submat_dist_and_blk_size
1068 
1069 ! **************************************************************************************************
1070 !> \brief Returns the projector on the unperturbed unoccupied ground state Q = 1 - SP on the block
1071 !> diagonal of a matrix with the standard size and distribution.
1072 !> \param proj_Q the matrix with the projector
1073 !> \param donor_state ...
1074 !> \param do_os whether it is open-shell calculation
1075 !> \param xas_tdp_env ...
1076 !> \param do_sf whether the projector should be prepared for spin-flip excitations
1077 !> \note In the spin-flip case, the alpha spins are sent to beta and vice-versa. The structure of
1078 !> the projector is changed by swapping the alpha-alpha with the beta-beta block, which
1079 !> naturally take the spin change into account. Only for open-shell.
1080 ! **************************************************************************************************
1081  SUBROUTINE get_q_projector(proj_Q, donor_state, do_os, xas_tdp_env, do_sf)
1082 
1083  TYPE(dbcsr_type), INTENT(INOUT) :: proj_q
1084  TYPE(donor_state_type), POINTER :: donor_state
1085  LOGICAL, INTENT(IN) :: do_os
1086  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1087  LOGICAL, INTENT(IN), OPTIONAL :: do_sf
1088 
1089  CHARACTER(len=*), PARAMETER :: routinen = 'get_q_projector'
1090 
1091  INTEGER :: blk, handle, iblk, imo, ispin, jblk, &
1092  nblk_row, ndo_mo, nspins
1093  INTEGER, DIMENSION(:), POINTER :: blk_size_q, row_blk_size
1094  LOGICAL :: found_block, my_dosf
1095  REAL(dp), DIMENSION(:), POINTER :: work_block
1096  TYPE(dbcsr_distribution_type), POINTER :: dist_q
1097  TYPE(dbcsr_iterator_type) :: iter
1098  TYPE(dbcsr_type), POINTER :: one_sp
1099 
1100  NULLIFY (work_block, one_sp, row_blk_size, dist_q, blk_size_q)
1101 
1102  CALL timeset(routinen, handle)
1103 
1104 ! Initialization
1105  nspins = 1; IF (do_os) nspins = 2
1106  ndo_mo = donor_state%ndo_mo
1107  one_sp => xas_tdp_env%q_projector(1)%matrix
1108  CALL dbcsr_get_info(one_sp, row_blk_size=row_blk_size)
1109  nblk_row = SIZE(row_blk_size)
1110  my_dosf = .false.
1111  IF (PRESENT(do_sf)) my_dosf = do_sf
1112  dist_q => donor_state%dbcsr_dist
1113  blk_size_q => donor_state%blk_size
1114 
1115  ! the projector is not symmetric
1116  CALL dbcsr_create(matrix=proj_q, name="PROJ Q", matrix_type=dbcsr_type_no_symmetry, dist=dist_q, &
1117  row_blk_size=blk_size_q, col_blk_size=blk_size_q)
1118 
1119 ! Fill the projector by looping over 1-SP and duplicating blocks. (all on the spin-MO block diagonal)
1120  DO ispin = 1, nspins
1121  one_sp => xas_tdp_env%q_projector(ispin)%matrix
1122 
1123  !if spin-flip, swap the alpha-alpha and beta-beta blocks
1124  IF (my_dosf) one_sp => xas_tdp_env%q_projector(3 - ispin)%matrix
1125 
1126  CALL dbcsr_iterator_start(iter, one_sp)
1127  DO WHILE (dbcsr_iterator_blocks_left(iter))
1128 
1129  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
1130 
1131  ! get the block
1132  CALL dbcsr_get_block_p(one_sp, iblk, jblk, work_block, found_block)
1133 
1134  IF (found_block) THEN
1135 
1136  DO imo = 1, ndo_mo
1137  CALL dbcsr_put_block(proj_q, ((ispin - 1)*ndo_mo + imo - 1)*nblk_row + iblk, &
1138  ((ispin - 1)*ndo_mo + imo - 1)*nblk_row + jblk, work_block)
1139  END DO
1140 
1141  END IF
1142  NULLIFY (work_block)
1143 
1144  END DO !iterator
1145  CALL dbcsr_iterator_stop(iter)
1146  END DO !ispin
1147 
1148  CALL dbcsr_finalize(proj_q)
1149 
1150  CALL timestop(handle)
1151 
1152  END SUBROUTINE get_q_projector
1153 
1154 ! **************************************************************************************************
1155 !> \brief Builds the matrix containing the ground state contribution to the matrix_tdp (aka matrix A)
1156 !> => A_{pis,qjt} = (F_pq*delta_ij - epsilon_ij*S_pq) delta_st, where:
1157 !> F is the KS matrix
1158 !> S is the overlap matrix
1159 !> epsilon_ij is the donor MO eigenvalue
1160 !> i,j labels the MOs, p,q the AOs and s,t the spins
1161 !> \param matrix_a pointer to a DBCSR matrix containing A
1162 !> \param donor_state ...
1163 !> \param do_os ...
1164 !> \param qs_env ...
1165 !> \param do_sf whether the ground state contribution should accommodate spin-flip
1166 !> \note Even localized non-canonical MOs are diagonalized in their subsapce => eps_ij = eps_ii*delta_ij
1167 !> Use GW2X corrected evals as eps_ii. If not GW2X correction, these are the default KS energies
1168 ! **************************************************************************************************
1169  SUBROUTINE build_gs_contribution(matrix_a, donor_state, do_os, qs_env, do_sf)
1170 
1171  TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a
1172  TYPE(donor_state_type), POINTER :: donor_state
1173  LOGICAL, INTENT(IN) :: do_os
1174  TYPE(qs_environment_type), POINTER :: qs_env
1175  LOGICAL, INTENT(IN), OPTIONAL :: do_sf
1176 
1177  CHARACTER(len=*), PARAMETER :: routinen = 'build_gs_contribution'
1178 
1179  INTEGER :: blk, handle, iblk, imo, ispin, jblk, &
1180  nblk_row, ndo_mo, nspins
1181  INTEGER, DIMENSION(:), POINTER :: blk_size_a, row_blk_size
1182  LOGICAL :: found_block, my_dosf
1183  REAL(dp), DIMENSION(:), POINTER :: work_block
1184  TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist, dist_a
1185  TYPE(dbcsr_iterator_type) :: iter
1186  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: m_ks, matrix_ks, matrix_s
1187  TYPE(dbcsr_type) :: work_matrix
1188 
1189  NULLIFY (matrix_ks, dbcsr_dist, row_blk_size, work_block, matrix_s, m_ks)
1190  NULLIFY (dist_a, blk_size_a)
1191 
1192  ! Note: matrix A is symmetric and block diagonal. If ADMM, the ks matrix is the total one,
1193  ! and it is corrected for eigenvalues (done at xas_tdp_init)
1194 
1195  CALL timeset(routinen, handle)
1196 
1197 ! Initialization
1198  nspins = 1; IF (do_os) nspins = 2
1199  ndo_mo = donor_state%ndo_mo
1200  CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks, matrix_s=matrix_s, dbcsr_dist=dbcsr_dist)
1201  CALL dbcsr_get_info(matrix_s(1)%matrix, row_blk_size=row_blk_size)
1202  nblk_row = SIZE(row_blk_size)
1203  dist_a => donor_state%dbcsr_dist
1204  blk_size_a => donor_state%blk_size
1205 
1206 ! Prepare the KS matrix pointer
1207  ALLOCATE (m_ks(nspins))
1208  m_ks(1)%matrix => matrix_ks(1)%matrix
1209  IF (do_os) m_ks(2)%matrix => matrix_ks(2)%matrix
1210 
1211 ! If spin-flip, swap the KS alpha-alpha and beta-beta quadrants.
1212  my_dosf = .false.
1213  IF (PRESENT(do_sf)) my_dosf = do_sf
1214  IF (my_dosf .AND. do_os) THEN
1215  m_ks(1)%matrix => matrix_ks(2)%matrix
1216  m_ks(2)%matrix => matrix_ks(1)%matrix
1217  END IF
1218 
1219 ! Creating the symmetric matrix A (and work)
1220  CALL dbcsr_create(matrix=matrix_a, name="MATRIX A", matrix_type=dbcsr_type_symmetric, &
1221  dist=dist_a, row_blk_size=blk_size_a, col_blk_size=blk_size_a)
1222  CALL dbcsr_create(matrix=work_matrix, name="WORK MAT", matrix_type=dbcsr_type_symmetric, &
1223  dist=dist_a, row_blk_size=blk_size_a, col_blk_size=blk_size_a)
1224 
1225  DO ispin = 1, nspins
1226 
1227 ! Loop over the blocks of KS and put them on the spin-MO block diagonal of matrix A
1228  CALL dbcsr_iterator_start(iter, m_ks(ispin)%matrix)
1229  DO WHILE (dbcsr_iterator_blocks_left(iter))
1230 
1231  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
1232 
1233 ! Get the block
1234  CALL dbcsr_get_block_p(m_ks(ispin)%matrix, iblk, jblk, work_block, found_block)
1235 
1236  IF (found_block) THEN
1237 
1238 ! The KS matrix only appears on diagonal of matrix A => loop over II donor MOs
1239  DO imo = 1, ndo_mo
1240 
1241 ! Put the block as it is
1242  CALL dbcsr_put_block(matrix_a, ((ispin - 1)*ndo_mo + imo - 1)*nblk_row + iblk, &
1243  ((ispin - 1)*ndo_mo + imo - 1)*nblk_row + jblk, work_block)
1244 
1245  END DO !imo
1246  END IF !found_block
1247  NULLIFY (work_block)
1248  END DO ! iteration on KS blocks
1249  CALL dbcsr_iterator_stop(iter)
1250 
1251 ! Loop over the blocks of S and put them on the block diagonal of work
1252 
1253  CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
1254  DO WHILE (dbcsr_iterator_blocks_left(iter))
1255 
1256  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
1257 
1258 ! Get the block
1259  CALL dbcsr_get_block_p(matrix_s(1)%matrix, iblk, jblk, work_block, found_block)
1260 
1261  IF (found_block) THEN
1262 
1263 ! Add S matrix on block diagonal as epsilon_ii*S_pq
1264  DO imo = 1, ndo_mo
1265 
1266  CALL dbcsr_put_block(work_matrix, ((ispin - 1)*ndo_mo + imo - 1)*nblk_row + iblk, &
1267  ((ispin - 1)*ndo_mo + imo - 1)*nblk_row + jblk, &
1268  donor_state%gw2x_evals(imo, ispin)*work_block)
1269  END DO !imo
1270  END IF !found block
1271  NULLIFY (work_block)
1272  END DO ! iteration on S blocks
1273  CALL dbcsr_iterator_stop(iter)
1274 
1275  END DO !ispin
1276  CALL dbcsr_finalize(matrix_a)
1277  CALL dbcsr_finalize(work_matrix)
1278 
1279 ! Take matrix_a = matrix_a - work
1280  CALL dbcsr_add(matrix_a, work_matrix, 1.0_dp, -1.0_dp)
1281  CALL dbcsr_finalize(matrix_a)
1282 
1283 ! Clean-up
1284  CALL dbcsr_release(work_matrix)
1285  DEALLOCATE (m_ks)
1286 
1287  CALL timestop(handle)
1288 
1289  END SUBROUTINE build_gs_contribution
1290 
1291 ! **************************************************************************************************
1292 !> \brief Creates the metric (aka matrix G) needed for the generalized eigenvalue problem and inverse
1293 !> => G_{pis,qjt} = S_pq*delta_ij*delta_st
1294 !> \param matrix_g dbcsr matrix containing G
1295 !> \param donor_state ...
1296 !> \param qs_env ...
1297 !> \param do_os if open-shell calculation
1298 !> \param do_inv if the inverse of G should be computed
1299 ! **************************************************************************************************
1300  SUBROUTINE build_metric(matrix_g, donor_state, qs_env, do_os, do_inv)
1301 
1302  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_g
1303  TYPE(donor_state_type), POINTER :: donor_state
1304  TYPE(qs_environment_type), POINTER :: qs_env
1305  LOGICAL, INTENT(IN) :: do_os
1306  LOGICAL, INTENT(IN), OPTIONAL :: do_inv
1307 
1308  CHARACTER(len=*), PARAMETER :: routinen = 'build_metric'
1309 
1310  INTEGER :: blk, handle, i, iblk, jblk, nao, &
1311  nblk_row, ndo_mo, nspins
1312  INTEGER, DIMENSION(:), POINTER :: blk_size_g, row_blk_size
1313  LOGICAL :: found_block, my_do_inv
1314  REAL(dp), DIMENSION(:), POINTER :: work_block
1315  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1316  TYPE(dbcsr_distribution_type), POINTER :: dist_g
1317  TYPE(dbcsr_iterator_type) :: iter
1318  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1319  TYPE(dbcsr_type) :: matrix_sinv
1320  TYPE(mp_para_env_type), POINTER :: para_env
1321 
1322  NULLIFY (matrix_s, row_blk_size, work_block, para_env, blacs_env, dist_g, blk_size_g)
1323 
1324  CALL timeset(routinen, handle)
1325 
1326 ! Initialization
1327  nspins = 1; IF (do_os) nspins = 2
1328  ndo_mo = donor_state%ndo_mo
1329  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
1330  CALL dbcsr_get_info(matrix_s(1)%matrix, row_blk_size=row_blk_size, nfullrows_total=nao)
1331  nblk_row = SIZE(row_blk_size)
1332  my_do_inv = .false.
1333  IF (PRESENT(do_inv)) my_do_inv = do_inv
1334  dist_g => donor_state%dbcsr_dist
1335  blk_size_g => donor_state%blk_size
1336 
1337 ! Creating the symmetric matrices G and G^-1 with the right size and distribution
1338  ALLOCATE (matrix_g(1)%matrix)
1339  CALL dbcsr_create(matrix=matrix_g(1)%matrix, name="MATRIX G", matrix_type=dbcsr_type_symmetric, &
1340  dist=dist_g, row_blk_size=blk_size_g, col_blk_size=blk_size_g)
1341 
1342 ! Fill the matrices G by looping over the block of S and putting them on the diagonal
1343  CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
1344  DO WHILE (dbcsr_iterator_blocks_left(iter))
1345 
1346  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
1347 
1348 ! Get the block
1349  CALL dbcsr_get_block_p(matrix_s(1)%matrix, iblk, jblk, work_block, found_block)
1350 
1351  IF (found_block) THEN
1352 
1353 ! Go over the diagonal of G => donor MOs ii, spin ss
1354  DO i = 1, ndo_mo*nspins
1355  CALL dbcsr_put_block(matrix_g(1)%matrix, (i - 1)*nblk_row + iblk, (i - 1)*nblk_row + jblk, work_block)
1356  END DO
1357 
1358  END IF
1359  NULLIFY (work_block)
1360 
1361  END DO ! dbcsr_iterator
1362  CALL dbcsr_iterator_stop(iter)
1363 
1364 ! Finalize
1365  CALL dbcsr_finalize(matrix_g(1)%matrix)
1366 
1367 ! If the inverse of G is required, do the same as above with the inverse
1368  IF (my_do_inv) THEN
1369 
1370  cpassert(SIZE(matrix_g) == 2)
1371 
1372  ! Create the matrix
1373  ALLOCATE (matrix_g(2)%matrix)
1374  CALL dbcsr_create(matrix=matrix_g(2)%matrix, name="MATRIX GINV", &
1375  matrix_type=dbcsr_type_symmetric, dist=dist_g, &
1376  row_blk_size=blk_size_g, col_blk_size=blk_size_g)
1377 
1378  ! Invert the overlap matrix
1379  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1380  CALL dbcsr_copy(matrix_sinv, matrix_s(1)%matrix)
1381  CALL cp_dbcsr_cholesky_decompose(matrix_sinv, para_env=para_env, blacs_env=blacs_env)
1382  CALL cp_dbcsr_cholesky_invert(matrix_sinv, para_env=para_env, blacs_env=blacs_env, upper_to_full=.true.)
1383 
1384 ! Fill the matrices G^-1 by looping over the block of S^-1 and putting them on the diagonal
1385  CALL dbcsr_iterator_start(iter, matrix_sinv)
1386  DO WHILE (dbcsr_iterator_blocks_left(iter))
1387 
1388  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
1389 
1390 ! Get the block
1391  CALL dbcsr_get_block_p(matrix_sinv, iblk, jblk, work_block, found_block)
1392 
1393  IF (found_block) THEN
1394 
1395 ! Go over the diagonal of G => donor MOs ii spin ss
1396  DO i = 1, ndo_mo*nspins
1397  CALL dbcsr_put_block(matrix_g(2)%matrix, (i - 1)*nblk_row + iblk, (i - 1)*nblk_row + jblk, work_block)
1398  END DO
1399 
1400  END IF
1401  NULLIFY (work_block)
1402 
1403  END DO ! dbcsr_iterator
1404  CALL dbcsr_iterator_stop(iter)
1405 
1406  ! Finalize
1407  CALL dbcsr_finalize(matrix_g(2)%matrix)
1408 
1409  ! Clean-up
1410  CALL dbcsr_release(matrix_sinv)
1411  END IF !do_inv
1412 
1413  CALL timestop(handle)
1414 
1415  END SUBROUTINE build_metric
1416 
1417 ! **************************************************************************************************
1418 !> \brief Builds the auxiliary matrix (A-D+E)^+0.5 needed for the transofrmation of the
1419 !> full-TDDFT problem into an Hermitian one
1420 !> \param threshold a threshold for allowed negative eigenvalues
1421 !> \param sx the amount of exact exchange
1422 !> \param matrix_a the ground state contribution matrix A
1423 !> \param matrix_d the on-diagonal exchange kernel matrix (ab|IJ)
1424 !> \param matrix_e the off-diagonal exchange kernel matrix (aJ|Ib)
1425 !> \param do_hfx if exact exchange is included
1426 !> \param proj_Q ...
1427 !> \param work ...
1428 !> \param donor_state ...
1429 !> \param eps_filter for the dbcsr multiplication
1430 !> \param qs_env ...
1431 ! **************************************************************************************************
1432  SUBROUTINE build_aux_matrix(threshold, sx, matrix_a, matrix_d, matrix_e, do_hfx, proj_Q, &
1433  work, donor_state, eps_filter, qs_env)
1434 
1435  REAL(dp), INTENT(IN) :: threshold, sx
1436  TYPE(dbcsr_type), INTENT(INOUT) :: matrix_a, matrix_d, matrix_e
1437  LOGICAL, INTENT(IN) :: do_hfx
1438  TYPE(dbcsr_type), INTENT(INOUT) :: proj_q, work
1439  TYPE(donor_state_type), POINTER :: donor_state
1440  REAL(dp), INTENT(IN) :: eps_filter
1441  TYPE(qs_environment_type), POINTER :: qs_env
1442 
1443  CHARACTER(len=*), PARAMETER :: routinen = 'build_aux_matrix'
1444 
1445  INTEGER :: handle, ndep
1446  REAL(dp) :: evals(2)
1447  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1448  TYPE(dbcsr_type) :: tmp_mat
1449  TYPE(mp_para_env_type), POINTER :: para_env
1450 
1451  NULLIFY (blacs_env, para_env)
1452 
1453  CALL timeset(routinen, handle)
1454 
1455  CALL dbcsr_copy(tmp_mat, matrix_a)
1456  IF (do_hfx) THEN
1457  CALL dbcsr_add(tmp_mat, matrix_d, 1.0_dp, -1.0_dp*sx) !scaled hfx
1458  CALL dbcsr_add(tmp_mat, matrix_e, 1.0_dp, 1.0_dp*sx)
1459  END IF
1460 
1461  ! Take the product with the Q projector:
1462  CALL dbcsr_multiply('N', 'N', 1.0_dp, proj_q, tmp_mat, 0.0_dp, work, filter_eps=eps_filter)
1463  CALL dbcsr_multiply('N', 'T', 1.0_dp, work, proj_q, 0.0_dp, tmp_mat, filter_eps=eps_filter)
1464 
1465  ! Actually computing and storing the auxiliary matrix
1466  ALLOCATE (donor_state%matrix_aux)
1467  CALL dbcsr_create(matrix=donor_state%matrix_aux, template=matrix_a, name="MAT AUX")
1468 
1469  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1470 
1471  ! good quality sqrt
1472  CALL cp_dbcsr_power(tmp_mat, 0.5_dp, threshold, ndep, para_env, blacs_env, eigenvalues=evals)
1473 
1474  CALL dbcsr_copy(donor_state%matrix_aux, tmp_mat)
1475 
1476  ! Warn the user if matrix not positive semi-definite
1477  IF (evals(1) < 0.0_dp .AND. abs(evals(1)) > threshold) THEN
1478  cpwarn("The full TDDFT problem might not have been soundly turned Hermitian. Try TDA.")
1479  END IF
1480 
1481  ! clean-up
1482  CALL dbcsr_release(tmp_mat)
1483 
1484  CALL timestop(handle)
1485 
1486  END SUBROUTINE build_aux_matrix
1487 
1488 ! **************************************************************************************************
1489 !> \brief Includes the SOC effects on the precomputed spin-conserving and spin-flip excitations
1490 !> from an open-shell calculation (UKS or ROKS). This is a perturbative treatment
1491 !> \param donor_state ...
1492 !> \param xas_tdp_env ...
1493 !> \param xas_tdp_control ...
1494 !> \param qs_env ...
1495 !> \note Using AMEWs, build an hermitian matrix with all excited states SOC coupling + the
1496 !> excitation energies on the diagonal. Then diagonalize it to get the new excitation
1497 !> energies and corresponding linear combinations of lr_coeffs.
1498 !> The AMEWs are normalized
1499 !> Only for open-shell calculations
1500 ! **************************************************************************************************
1501  SUBROUTINE include_os_soc(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
1502 
1503  TYPE(donor_state_type), POINTER :: donor_state
1504  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1505  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1506  TYPE(qs_environment_type), POINTER :: qs_env
1507 
1508  CHARACTER(len=*), PARAMETER :: routinen = 'include_os_soc'
1509 
1510  INTEGER :: group, handle, homo, iex, isc, isf, nao, &
1511  ndo_mo, ndo_so, nex, npcols, nprows, &
1512  nsc, nsf, ntot, tas(2), tbs(2)
1513  INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_dist, row_blk_size, &
1514  row_dist, row_dist_new
1515  INTEGER, DIMENSION(:, :), POINTER :: pgrid
1516  LOGICAL :: do_roks, do_uks
1517  REAL(dp) :: eps_filter, gs_sum, soc
1518  REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag, tmp_evals
1519  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: domo_soc_x, domo_soc_y, domo_soc_z, &
1520  gsex_block
1521  REAL(dp), DIMENSION(:), POINTER :: sc_evals, sf_evals
1522  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1523  TYPE(cp_cfm_type) :: evecs_cfm, pert_cfm
1524  TYPE(cp_fm_struct_type), POINTER :: full_struct, gsex_struct, prod_struct, &
1525  vec_struct, work_struct
1526  TYPE(cp_fm_type) :: gsex_fm, img_fm, prod_work, real_fm, &
1527  vec_soc_x, vec_soc_y, vec_soc_z, &
1528  vec_work, work_fm
1529  TYPE(cp_fm_type), POINTER :: gs_coeffs, mo_coeff, sc_coeffs, sf_coeffs
1530  TYPE(dbcsr_distribution_type), POINTER :: coeffs_dist, dbcsr_dist, prod_dist
1531  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1532  TYPE(dbcsr_soc_package_type) :: dbcsr_soc_package
1533  TYPE(dbcsr_type), POINTER :: dbcsr_ovlp, dbcsr_prod, dbcsr_sc, &
1534  dbcsr_sf, dbcsr_tmp, dbcsr_work, &
1535  orb_soc_x, orb_soc_y, orb_soc_z
1536  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1537  TYPE(mp_para_env_type), POINTER :: para_env
1538 
1539  NULLIFY (gs_coeffs, sc_coeffs, sf_coeffs, matrix_s, orb_soc_x, orb_soc_y, orb_soc_z, mos)
1540  NULLIFY (full_struct, para_env, blacs_env, mo_coeff, sc_evals, sf_evals, vec_struct, prod_struct)
1541  NULLIFY (work_struct, gsex_struct, col_dist, row_dist)
1542  NULLIFY (col_blk_size, row_blk_size, row_dist_new, pgrid, dbcsr_sc, dbcsr_sf, dbcsr_work)
1543  NULLIFY (dbcsr_tmp, dbcsr_ovlp, dbcsr_prod)
1544 
1545  CALL timeset(routinen, handle)
1546 
1547 ! Initialization
1548  sc_coeffs => donor_state%sc_coeffs
1549  sf_coeffs => donor_state%sf_coeffs
1550  sc_evals => donor_state%sc_evals
1551  sf_evals => donor_state%sf_evals
1552  nsc = SIZE(sc_evals)
1553  nsf = SIZE(sf_evals)
1554  ntot = 1 + nsc + nsf
1555  nex = nsc !by contrutciotn nsc == nsf, but keep 2 counts for clarity
1556  ndo_mo = donor_state%ndo_mo
1557  ndo_so = 2*ndo_mo
1558  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env, mos=mos, matrix_s=matrix_s)
1559  CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
1560  orb_soc_x => xas_tdp_env%orb_soc(1)%matrix
1561  orb_soc_y => xas_tdp_env%orb_soc(2)%matrix
1562  orb_soc_z => xas_tdp_env%orb_soc(3)%matrix
1563  do_roks = xas_tdp_control%do_roks
1564  do_uks = xas_tdp_control%do_uks
1565  eps_filter = xas_tdp_control%eps_filter
1566 
1567  ! For the GS coeffs, we use the same structure both for ROKS and UKS here => allows us to write
1568  ! general code later on, and not use IF (do_roks) statements every second line
1569  IF (do_uks) gs_coeffs => donor_state%gs_coeffs
1570  IF (do_roks) THEN
1571  CALL cp_fm_struct_create(vec_struct, context=blacs_env, para_env=para_env, &
1572  nrow_global=nao, ncol_global=ndo_so)
1573  ALLOCATE (gs_coeffs)
1574  CALL cp_fm_create(gs_coeffs, vec_struct)
1575 
1576  ! only alpha donor MOs are stored, need to copy them intoboth the alpha and the beta slot
1577  CALL cp_fm_to_fm_submat(msource=donor_state%gs_coeffs, mtarget=gs_coeffs, nrow=nao, &
1578  ncol=ndo_mo, s_firstrow=1, s_firstcol=1, t_firstrow=1, &
1579  t_firstcol=1)
1580  CALL cp_fm_to_fm_submat(msource=donor_state%gs_coeffs, mtarget=gs_coeffs, nrow=nao, &
1581  ncol=ndo_mo, s_firstrow=1, s_firstcol=1, t_firstrow=1, &
1582  t_firstcol=ndo_mo + 1)
1583 
1584  CALL cp_fm_struct_release(vec_struct)
1585  END IF
1586 
1587 ! Creating the real and the imaginary part of the SOC perturbation matrix
1588  CALL cp_fm_struct_create(full_struct, context=blacs_env, para_env=para_env, &
1589  nrow_global=ntot, ncol_global=ntot)
1590  CALL cp_fm_create(real_fm, full_struct)
1591  CALL cp_fm_create(img_fm, full_struct)
1592 
1593 ! Put the excitation energies on the diagonal of the real matrix. Element 1,1 is the ground state
1594  DO isc = 1, nsc
1595  CALL cp_fm_set_element(real_fm, 1 + isc, 1 + isc, sc_evals(isc))
1596  END DO
1597  DO isf = 1, nsf
1598  CALL cp_fm_set_element(real_fm, 1 + nsc + isf, 1 + nsc + isf, sf_evals(isf))
1599  END DO
1600 
1601 ! Create the bdcsr machinery
1602  CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist)
1603  CALL dbcsr_distribution_get(dbcsr_dist, group=group, row_dist=row_dist, pgrid=pgrid, &
1604  npcols=npcols, nprows=nprows)
1605  ALLOCATE (col_dist(nex), row_dist_new(nex))
1606  DO iex = 1, nex
1607  col_dist(iex) = modulo(npcols - iex, npcols)
1608  row_dist_new(iex) = modulo(nprows - iex, nprows)
1609  END DO
1610  ALLOCATE (coeffs_dist, prod_dist)
1611  CALL dbcsr_distribution_new(coeffs_dist, group=group, pgrid=pgrid, row_dist=row_dist, &
1612  col_dist=col_dist)
1613  CALL dbcsr_distribution_new(prod_dist, group=group, pgrid=pgrid, row_dist=row_dist_new, &
1614  col_dist=col_dist)
1615 
1616  !Create the matrices
1617  ALLOCATE (col_blk_size(nex))
1618  col_blk_size = ndo_so
1619  CALL dbcsr_get_info(matrix_s(1)%matrix, row_blk_size=row_blk_size)
1620 
1621  ALLOCATE (dbcsr_sc, dbcsr_sf, dbcsr_work, dbcsr_ovlp, dbcsr_tmp, dbcsr_prod)
1622  CALL dbcsr_create(matrix=dbcsr_sc, name="SPIN CONS", matrix_type=dbcsr_type_no_symmetry, &
1623  dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
1624  CALL dbcsr_create(matrix=dbcsr_sf, name="SPIN FLIP", matrix_type=dbcsr_type_no_symmetry, &
1625  dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
1626  CALL dbcsr_create(matrix=dbcsr_work, name="WORK", matrix_type=dbcsr_type_no_symmetry, &
1627  dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
1628  CALL dbcsr_create(matrix=dbcsr_prod, name="PROD", matrix_type=dbcsr_type_no_symmetry, &
1629  dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
1630  CALL dbcsr_create(matrix=dbcsr_ovlp, name="OVLP", matrix_type=dbcsr_type_no_symmetry, &
1631  dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
1632 
1633  col_blk_size = 1
1634  CALL dbcsr_create(matrix=dbcsr_tmp, name="TMP", matrix_type=dbcsr_type_no_symmetry, &
1635  dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
1636  CALL dbcsr_reserve_all_blocks(dbcsr_tmp)
1637 
1638  dbcsr_soc_package%dbcsr_sc => dbcsr_sc
1639  dbcsr_soc_package%dbcsr_sf => dbcsr_sf
1640  dbcsr_soc_package%dbcsr_work => dbcsr_work
1641  dbcsr_soc_package%dbcsr_ovlp => dbcsr_ovlp
1642  dbcsr_soc_package%dbcsr_prod => dbcsr_prod
1643  dbcsr_soc_package%dbcsr_tmp => dbcsr_tmp
1644 
1645  !Filling the coeffs matrices by copying from the stored fms
1646  CALL copy_fm_to_dbcsr(sc_coeffs, dbcsr_sc)
1647  CALL copy_fm_to_dbcsr(sf_coeffs, dbcsr_sf)
1648 
1649 ! Precompute what we can before looping over excited states.
1650  ! Need to compute the scalar: sum_i sum_sigma <phi^0_i,sigma|SOC|phi^0_i,sigma>, where all
1651  ! occupied MOs are taken into account
1652 
1653  !start with the alpha MOs
1654  CALL get_mo_set(mos(1), mo_coeff=mo_coeff, homo=homo)
1655  ALLOCATE (diag(homo))
1656  CALL cp_fm_get_info(mo_coeff, matrix_struct=vec_struct)
1657  CALL cp_fm_struct_create(prod_struct, context=blacs_env, para_env=para_env, &
1658  nrow_global=homo, ncol_global=homo)
1659  CALL cp_fm_create(vec_work, vec_struct)
1660  CALL cp_fm_create(prod_work, prod_struct)
1661 
1662  ! <alpha|SOC_z|alpha> => spin integration yields +1
1663  CALL cp_dbcsr_sm_fm_multiply(orb_soc_z, mo_coeff, vec_work, ncol=homo)
1664  CALL parallel_gemm('T', 'N', homo, homo, nao, 1.0_dp, mo_coeff, vec_work, 0.0_dp, prod_work)
1665  CALL cp_fm_get_diag(prod_work, diag)
1666  gs_sum = sum(diag)
1667 
1668  CALL cp_fm_release(vec_work)
1669  CALL cp_fm_release(prod_work)
1670  CALL cp_fm_struct_release(prod_struct)
1671  DEALLOCATE (diag)
1672  NULLIFY (vec_struct)
1673 
1674  ! Now do the same with the beta gs coeffs
1675  CALL get_mo_set(mos(2), mo_coeff=mo_coeff, homo=homo)
1676  ALLOCATE (diag(homo))
1677  CALL cp_fm_get_info(mo_coeff, matrix_struct=vec_struct)
1678  CALL cp_fm_struct_create(prod_struct, context=blacs_env, para_env=para_env, &
1679  nrow_global=homo, ncol_global=homo)
1680  CALL cp_fm_create(vec_work, vec_struct)
1681  CALL cp_fm_create(prod_work, prod_struct)
1682 
1683  ! <beta|SOC_z|beta> => spin integration yields -1
1684  CALL cp_dbcsr_sm_fm_multiply(orb_soc_z, mo_coeff, vec_work, ncol=homo)
1685  CALL parallel_gemm('T', 'N', homo, homo, nao, 1.0_dp, mo_coeff, vec_work, 0.0_dp, prod_work)
1686  CALL cp_fm_get_diag(prod_work, diag)
1687  gs_sum = gs_sum - sum(diag) ! -1 because of spin integration
1688 
1689  CALL cp_fm_release(vec_work)
1690  CALL cp_fm_release(prod_work)
1691  CALL cp_fm_struct_release(prod_struct)
1692  DEALLOCATE (diag)
1693 
1694  ! Need to compute: <phi^0_Isigma|SOC|phi^0_Jtau> for the donor MOs
1695 
1696  CALL cp_fm_struct_create(vec_struct, context=blacs_env, para_env=para_env, &
1697  nrow_global=nao, ncol_global=ndo_so)
1698  CALL cp_fm_struct_create(prod_struct, context=blacs_env, para_env=para_env, &
1699  nrow_global=ndo_so, ncol_global=ndo_so)
1700  CALL cp_fm_create(vec_soc_x, vec_struct) ! for SOC_x*gs_coeffs
1701  CALL cp_fm_create(vec_soc_y, vec_struct) ! for SOC_y*gs_coeffs
1702  CALL cp_fm_create(vec_soc_z, vec_struct) ! for SOC_z*gs_coeffs
1703  CALL cp_fm_create(prod_work, prod_struct)
1704  ALLOCATE (diag(ndo_so))
1705 
1706  ALLOCATE (domo_soc_x(ndo_so, ndo_so), domo_soc_y(ndo_so, ndo_so), domo_soc_z(ndo_so, ndo_so))
1707 
1708  CALL cp_dbcsr_sm_fm_multiply(orb_soc_x, gs_coeffs, vec_soc_x, ncol=ndo_so)
1709  CALL parallel_gemm('T', 'N', ndo_so, ndo_so, nao, 1.0_dp, gs_coeffs, vec_soc_x, 0.0_dp, prod_work)
1710  CALL cp_fm_get_submatrix(prod_work, domo_soc_x)
1711 
1712  CALL cp_dbcsr_sm_fm_multiply(orb_soc_y, gs_coeffs, vec_soc_y, ncol=ndo_so)
1713  CALL parallel_gemm('T', 'N', ndo_so, ndo_so, nao, 1.0_dp, gs_coeffs, vec_soc_y, 0.0_dp, prod_work)
1714  CALL cp_fm_get_submatrix(prod_work, domo_soc_y)
1715 
1716  CALL cp_dbcsr_sm_fm_multiply(orb_soc_z, gs_coeffs, vec_soc_z, ncol=ndo_so)
1717  CALL parallel_gemm('T', 'N', ndo_so, ndo_so, nao, 1.0_dp, gs_coeffs, vec_soc_z, 0.0_dp, prod_work)
1718  CALL cp_fm_get_submatrix(prod_work, domo_soc_z)
1719 
1720  ! some more useful work matrices
1721  CALL cp_fm_struct_create(work_struct, context=blacs_env, para_env=para_env, &
1722  nrow_global=nex, ncol_global=nex)
1723  CALL cp_fm_create(work_fm, work_struct)
1724 
1725 ! Looping over the excited states, computing the SOC and filling the perturbation matrix
1726 ! There are 3 loops to do: sc-sc, sc-sf and sf-sf
1727 ! The final perturbation matrix is Hermitian, only the upper diagonal is needed
1728 
1729  !need some work matrices for the GS stuff
1730  CALL cp_fm_struct_create(gsex_struct, context=blacs_env, para_env=para_env, &
1731  nrow_global=nex*ndo_so, ncol_global=ndo_so)
1732  CALL cp_fm_create(gsex_fm, gsex_struct)
1733  ALLOCATE (gsex_block(ndo_so, ndo_so))
1734 
1735 ! Start with ground-state/spin-conserving SOC:
1736  ! <Psi_0|SOC|Psi_Jsc> = sum_k,sigma <phi^0_k,sigma|SOC|phi^Jsc_k,sigma>
1737 
1738  !compute -sc_coeffs*SOC_Z*gs_coeffs, minus sign because SOC_z antisymmetric
1739  CALL parallel_gemm('T', 'N', nex*ndo_so, ndo_so, nao, -1.0_dp, sc_coeffs, vec_soc_z, 0.0_dp, gsex_fm)
1740 
1741  DO isc = 1, nsc
1742  CALL cp_fm_get_submatrix(fm=gsex_fm, target_m=gsex_block, start_row=(isc - 1)*ndo_so + 1, &
1743  start_col=1, n_rows=ndo_so, n_cols=ndo_so)
1744  diag(:) = get_diag(gsex_block)
1745  soc = sum(diag(1:ndo_mo)) - sum(diag(ndo_mo + 1:ndo_so)) !minus sign because of spin integration
1746 
1747  !purely imaginary contribution
1748  CALL cp_fm_set_element(img_fm, 1, 1 + isc, soc)
1749  END DO !isc
1750 
1751 ! Then ground-state/spin-flip SOC:
1752  !<Psi_0|SOC|Psi_Jsf> = sum_k,sigma <phi^0_k,sigma|SOC|phi^Jsc_k,tau> sigma != tau
1753 
1754  !compute -sc_coeffs*SOC_x*gs_coeffs, imaginary contribution
1755  CALL parallel_gemm('T', 'N', nex*ndo_so, ndo_so, nao, -1.0_dp, sc_coeffs, vec_soc_x, 0.0_dp, gsex_fm)
1756 
1757  DO isf = 1, nsf
1758  CALL cp_fm_get_submatrix(fm=gsex_fm, target_m=gsex_block, start_row=(isf - 1)*ndo_so + 1, &
1759  start_col=1, n_rows=ndo_so, n_cols=ndo_so)
1760  diag(:) = get_diag(gsex_block)
1761  soc = sum(diag) !alpha and beta parts are simply added due to spin integration
1762  CALL cp_fm_set_element(img_fm, 1, 1 + nsc + isf, soc)
1763  END DO !isf
1764 
1765  !compute -sc_coeffs*SOC_y*gs_coeffs, real contribution
1766  CALL parallel_gemm('T', 'N', nex*ndo_so, ndo_so, nao, -1.0_dp, sc_coeffs, vec_soc_y, 0.0_dp, gsex_fm)
1767 
1768  DO isf = 1, nsf
1769  CALL cp_fm_get_submatrix(fm=gsex_fm, target_m=gsex_block, start_row=(isf - 1)*ndo_so + 1, &
1770  start_col=1, n_rows=ndo_so, n_cols=ndo_so)
1771  diag(:) = get_diag(gsex_block)
1772  soc = sum(diag(1:ndo_mo)) ! alpha-beta
1773  soc = soc - sum(diag(ndo_mo + 1:ndo_so)) !beta-alpha
1774  CALL cp_fm_set_element(real_fm, 1, 1 + nsc + isf, soc)
1775  END DO !isf
1776 
1777  !ground-state cleanup
1778  CALL cp_fm_release(gsex_fm)
1779  CALL cp_fm_release(vec_soc_x)
1780  CALL cp_fm_release(vec_soc_y)
1781  CALL cp_fm_release(vec_soc_z)
1782  CALL cp_fm_release(prod_work)
1783  CALL cp_fm_struct_release(gsex_struct)
1784  CALL cp_fm_struct_release(prod_struct)
1785  CALL cp_fm_struct_release(vec_struct)
1786  DEALLOCATE (gsex_block)
1787 
1788 ! Then spin-conserving/spin-conserving SOC
1789 ! <Psi_Isc|SOC|Psi_Jsc> =
1790 ! sum_k,sigma [<psi^Isc_k,sigma|SOC|psi^Jsc_k,sigma> + <psi^Isc_k,sigma|psi^Jsc_k,sigma> * gs_sum]
1791 ! - sum_k,l,sigma <psi^0_k,sigma|SOC|psi^0_l,sigma> * <psi^Isc_l,sigma|psi^Jsc_k,sigma>
1792 
1793  !Same spin integration => only SOC_z matters, and the contribution is purely imaginary
1794  CALL dbcsr_multiply('N', 'N', 1.0_dp, orb_soc_z, dbcsr_sc, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
1795  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sc, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
1796 
1797  !the overlap as well
1798  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s(1)%matrix, dbcsr_sc, 0.0_dp, dbcsr_work, &
1799  filter_eps=eps_filter)
1800  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sc, dbcsr_work, 0.0_dp, dbcsr_ovlp, filter_eps=eps_filter)
1801 
1802  CALL os_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_soc_z, pref_diaga=1.0_dp, &
1803  pref_diagb=-1.0_dp, pref_tracea=-1.0_dp, pref_traceb=1.0_dp, &
1804  pref_diags=gs_sum, symmetric=.true.)
1805 
1806  CALL copy_dbcsr_to_fm(dbcsr_tmp, work_fm)
1807  CALL cp_fm_to_fm_submat(msource=work_fm, mtarget=img_fm, nrow=nex, ncol=nex, s_firstrow=1, &
1808  s_firstcol=1, t_firstrow=2, t_firstcol=2)
1809 
1810 ! Then spin-flip/spin-flip SOC
1811 ! <Psi_Isf|SOC|Psi_Jsf> =
1812 ! sum_k,sigma [<psi^Isf_k,tau|SOC|psi^Jsf_k,tau> + <psi^Isf_k,tau|psi^Jsf_k,tau> * gs_sum]
1813 ! - sum_k,l,sigma <psi^0_k,sigma|SOC|psi^0_l,sigma> * <psi^Isf_l,tau|psi^Jsf_k,tau> , tau != sigma
1814 
1815  !Same spin integration => only SOC_z matters, and the contribution is purely imaginary
1816  CALL dbcsr_multiply('N', 'N', 1.0_dp, orb_soc_z, dbcsr_sf, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
1817  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sf, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
1818 
1819  !the overlap as well
1820  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s(1)%matrix, dbcsr_sf, 0.0_dp, &
1821  dbcsr_work, filter_eps=eps_filter)
1822  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sf, dbcsr_work, 0.0_dp, dbcsr_ovlp, filter_eps=eps_filter)
1823 
1824  !note: the different prefactors are derived from the fact that because of spin-flip, we have
1825  !alpha-gs and beta-lr interaction
1826  CALL os_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_soc_z, pref_diaga=-1.0_dp, &
1827  pref_diagb=1.0_dp, pref_tracea=-1.0_dp, pref_traceb=1.0_dp, &
1828  pref_diags=gs_sum, symmetric=.true.)
1829 
1830  CALL copy_dbcsr_to_fm(dbcsr_tmp, work_fm)
1831  CALL cp_fm_to_fm_submat(msource=work_fm, mtarget=img_fm, nrow=nex, ncol=nex, s_firstrow=1, &
1832  s_firstcol=1, t_firstrow=1 + nsc + 1, t_firstcol=1 + nsc + 1)
1833 
1834 ! Finally the spin-conserving/spin-flip interaction
1835 ! <Psi_Isc|SOC|Psi_Jsf> = sum_k,sigma <psi^Isc_k,sigma|SOC|psi^Isf_k,tau>
1836 ! - sum_k,l,sigma <psi^0_k,tau|SOC|psi^0_l,sigma
1837 
1838  tas(1) = ndo_mo + 1; tbs(1) = 1
1839  tas(2) = 1; tbs(2) = ndo_mo + 1
1840 
1841  !the overlap
1842  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s(1)%matrix, dbcsr_sf, 0.0_dp, &
1843  dbcsr_work, filter_eps=eps_filter)
1844  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sc, dbcsr_work, 0.0_dp, dbcsr_ovlp, filter_eps=eps_filter)
1845 
1846  !start with the imaginary contribution
1847  CALL dbcsr_multiply('N', 'N', 1.0_dp, orb_soc_x, dbcsr_sc, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
1848  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sf, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
1849 
1850  CALL os_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_soc_x, pref_diaga=1.0_dp, &
1851  pref_diagb=1.0_dp, pref_tracea=-1.0_dp, pref_traceb=-1.0_dp, &
1852  tracea_start=tas, traceb_start=tbs)
1853 
1854  CALL copy_dbcsr_to_fm(dbcsr_tmp, work_fm)
1855  CALL cp_fm_to_fm_submat(msource=work_fm, mtarget=img_fm, nrow=nex, ncol=nex, s_firstrow=1, &
1856  s_firstcol=1, t_firstrow=2, t_firstcol=1 + nsc + 1)
1857 
1858  !follow up with the real contribution
1859  CALL dbcsr_multiply('N', 'N', 1.0_dp, orb_soc_y, dbcsr_sf, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
1860  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sc, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
1861 
1862  CALL os_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_soc_y, pref_diaga=1.0_dp, &
1863  pref_diagb=-1.0_dp, pref_tracea=1.0_dp, pref_traceb=-1.0_dp, &
1864  tracea_start=tas, traceb_start=tbs)
1865 
1866  CALL copy_dbcsr_to_fm(dbcsr_tmp, work_fm)
1867  CALL cp_fm_to_fm_submat(msource=work_fm, mtarget=real_fm, nrow=nex, ncol=nex, s_firstrow=1, &
1868  s_firstcol=1, t_firstrow=2, t_firstcol=1 + nsc + 1)
1869 
1870 ! Setting up the complex Hermitian perturbed matrix
1871  CALL cp_cfm_create(pert_cfm, full_struct)
1872  CALL cp_fm_to_cfm(real_fm, img_fm, pert_cfm)
1873 
1874  CALL cp_fm_release(real_fm)
1875  CALL cp_fm_release(img_fm)
1876 
1877 ! Diagonalize the perturbed matrix
1878  ALLOCATE (tmp_evals(ntot))
1879  CALL cp_cfm_create(evecs_cfm, full_struct)
1880  CALL cp_cfm_heevd(pert_cfm, evecs_cfm, tmp_evals)
1881 
1882  !shift the energies such that the GS has zero and store all that in soc_evals (\wo the GS)
1883  ALLOCATE (donor_state%soc_evals(ntot - 1))
1884  donor_state%soc_evals(:) = tmp_evals(2:ntot) - tmp_evals(1)
1885 
1886 ! The SOC dipole oscillator strengths
1887  CALL compute_soc_dipole_fosc(evecs_cfm, dbcsr_soc_package, donor_state, xas_tdp_env, &
1888  xas_tdp_control, qs_env, gs_coeffs=gs_coeffs)
1889 
1890 ! And quadrupole
1891  IF (xas_tdp_control%do_quad) THEN
1892  CALL compute_soc_quadrupole_fosc(evecs_cfm, dbcsr_soc_package, donor_state, xas_tdp_env, &
1893  xas_tdp_control, qs_env, gs_coeffs=gs_coeffs)
1894  END IF
1895 
1896 ! Clean-up
1897  CALL cp_cfm_release(pert_cfm)
1898  CALL cp_cfm_release(evecs_cfm)
1899  CALL cp_fm_struct_release(full_struct)
1900  IF (do_roks) THEN
1901  CALL cp_fm_release(gs_coeffs)
1902  DEALLOCATE (gs_coeffs)
1903  END IF
1904  CALL dbcsr_distribution_release(coeffs_dist)
1905  CALL dbcsr_distribution_release(prod_dist)
1906  CALL dbcsr_release(dbcsr_sc)
1907  CALL dbcsr_release(dbcsr_sf)
1908  CALL dbcsr_release(dbcsr_prod)
1909  CALL dbcsr_release(dbcsr_ovlp)
1910  CALL dbcsr_release(dbcsr_tmp)
1911  CALL dbcsr_release(dbcsr_work)
1912  CALL cp_fm_release(work_fm)
1913  CALL cp_fm_struct_release(work_struct)
1914  DEALLOCATE (coeffs_dist, prod_dist, col_dist, col_blk_size, row_dist_new)
1915  DEALLOCATE (dbcsr_sc, dbcsr_sf, dbcsr_work, dbcsr_prod, dbcsr_ovlp, dbcsr_tmp)
1916 
1917  CALL timestop(handle)
1918 
1919  END SUBROUTINE include_os_soc
1920 
1921 ! **************************************************************************************************
1922 !> \brief Includes the SOC effects on the precomputed restricted closed-shell singlet and triplet
1923 !> excitations. This is a perturbative treatmnent
1924 !> \param donor_state ...
1925 !> \param xas_tdp_env ...
1926 !> \param xas_tdp_control ...
1927 !> \param qs_env ...
1928 !> \note Using AMEWs, build an hermitian matrix with all excited states SOC coupling + the
1929 !> excitation energies on the diagonal. Then diagonalize it to get the new excitation
1930 !> energies and corresponding linear combinations of lr_coeffs.
1931 !> The AMEWs are normalized
1932 !> Only for spin-restricted calculations
1933 !> The ms=-1,+1 triplets are not explicitely computed in the first place. Assume they have
1934 !> the same energy as the ms=0 triplets and apply the spin raising and lowering operators
1935 !> on the latter to get their AMEWs => this is the qusi-degenerate perturbation theory
1936 !> approach by Neese (QDPT)
1937 ! **************************************************************************************************
1938  SUBROUTINE include_rcs_soc(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
1939 
1940  TYPE(donor_state_type), POINTER :: donor_state
1941  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1942  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1943  TYPE(qs_environment_type), POINTER :: qs_env
1944 
1945  CHARACTER(len=*), PARAMETER :: routinen = 'include_rcs_soc'
1946 
1947  INTEGER :: group, handle, iex, isg, itp, nao, &
1948  ndo_mo, nex, npcols, nprows, nsg, &
1949  ntot, ntp
1950  INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_dist, row_blk_size, &
1951  row_dist, row_dist_new
1952  INTEGER, DIMENSION(:, :), POINTER :: pgrid
1953  REAL(dp) :: eps_filter, soc_gst, sqrt2
1954  REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag, tmp_evals
1955  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: domo_soc_x, domo_soc_y, domo_soc_z, &
1956  gstp_block
1957  REAL(dp), DIMENSION(:), POINTER :: sg_evals, tp_evals
1958  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1959  TYPE(cp_cfm_type) :: evecs_cfm, hami_cfm
1960  TYPE(cp_fm_struct_type), POINTER :: full_struct, gstp_struct, prod_struct, &
1961  vec_struct, work_struct
1962  TYPE(cp_fm_type) :: gstp_fm, img_fm, prod_fm, real_fm, &
1963  tmp_fm, vec_soc_x, vec_soc_y, &
1964  vec_soc_z, work_fm
1965  TYPE(cp_fm_type), POINTER :: gs_coeffs, sg_coeffs, tp_coeffs
1966  TYPE(dbcsr_distribution_type), POINTER :: coeffs_dist, dbcsr_dist, prod_dist
1967  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1968  TYPE(dbcsr_soc_package_type) :: dbcsr_soc_package
1969  TYPE(dbcsr_type), POINTER :: dbcsr_ovlp, dbcsr_prod, dbcsr_sg, &
1970  dbcsr_tmp, dbcsr_tp, dbcsr_work, &
1971  orb_soc_x, orb_soc_y, orb_soc_z
1972  TYPE(mp_para_env_type), POINTER :: para_env
1973 
1974  NULLIFY (sg_coeffs, tp_coeffs, gs_coeffs, sg_evals, tp_evals, full_struct)
1975  NULLIFY (para_env, blacs_env, prod_struct, vec_struct, orb_soc_y, orb_soc_z)
1976  NULLIFY (matrix_s, orb_soc_x)
1977  NULLIFY (work_struct, dbcsr_dist, coeffs_dist, prod_dist, pgrid)
1978  NULLIFY (col_dist, row_dist, col_blk_size, row_blk_size, row_dist_new, gstp_struct)
1979  NULLIFY (dbcsr_tp, dbcsr_sg, dbcsr_prod, dbcsr_work, dbcsr_ovlp, dbcsr_tmp)
1980 
1981  CALL timeset(routinen, handle)
1982 
1983 ! Initialization
1984  cpassert(ASSOCIATED(xas_tdp_control))
1985  gs_coeffs => donor_state%gs_coeffs
1986  sg_coeffs => donor_state%sg_coeffs
1987  tp_coeffs => donor_state%tp_coeffs
1988  sg_evals => donor_state%sg_evals
1989  tp_evals => donor_state%tp_evals
1990  nsg = SIZE(sg_evals)
1991  ntp = SIZE(tp_evals)
1992  ntot = 1 + nsg + 3*ntp
1993  ndo_mo = donor_state%ndo_mo
1994  CALL get_qs_env(qs_env, matrix_s=matrix_s)
1995  CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
1996  orb_soc_x => xas_tdp_env%orb_soc(1)%matrix
1997  orb_soc_y => xas_tdp_env%orb_soc(2)%matrix
1998  orb_soc_z => xas_tdp_env%orb_soc(3)%matrix
1999  !by construction nsg == ntp, keep those separate for more code clarity though
2000  cpassert(nsg == ntp)
2001  nex = nsg
2002  eps_filter = xas_tdp_control%eps_filter
2003 
2004 ! Creating the real part and imaginary part of the final SOC fm
2005  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
2006  CALL cp_fm_struct_create(full_struct, context=blacs_env, para_env=para_env, &
2007  nrow_global=ntot, ncol_global=ntot)
2008  CALL cp_fm_create(real_fm, full_struct)
2009  CALL cp_fm_create(img_fm, full_struct)
2010 
2011 ! Put the excitation energies on the diagonal of the real matrix
2012  DO isg = 1, nsg
2013  CALL cp_fm_set_element(real_fm, 1 + isg, 1 + isg, sg_evals(isg))
2014  END DO
2015  DO itp = 1, ntp
2016  ! first T^-1, then T^0, then T^+1
2017  CALL cp_fm_set_element(real_fm, 1 + itp + nsg, 1 + itp + nsg, tp_evals(itp))
2018  CALL cp_fm_set_element(real_fm, 1 + itp + ntp + nsg, 1 + itp + ntp + nsg, tp_evals(itp))
2019  CALL cp_fm_set_element(real_fm, 1 + itp + 2*ntp + nsg, 1 + itp + 2*ntp + nsg, tp_evals(itp))
2020  END DO
2021 
2022 ! Create the dbcsr machinery (for fast MM, the core of this routine)
2023  CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist)
2024  CALL dbcsr_distribution_get(dbcsr_dist, group=group, row_dist=row_dist, pgrid=pgrid, &
2025  npcols=npcols, nprows=nprows)
2026  ALLOCATE (col_dist(nex), row_dist_new(nex))
2027  DO iex = 1, nex
2028  col_dist(iex) = modulo(npcols - iex, npcols)
2029  row_dist_new(iex) = modulo(nprows - iex, nprows)
2030  END DO
2031  ALLOCATE (coeffs_dist, prod_dist)
2032  CALL dbcsr_distribution_new(coeffs_dist, group=group, pgrid=pgrid, row_dist=row_dist, &
2033  col_dist=col_dist)
2034  CALL dbcsr_distribution_new(prod_dist, group=group, pgrid=pgrid, row_dist=row_dist_new, &
2035  col_dist=col_dist)
2036 
2037  !Create the matrices
2038  ALLOCATE (col_blk_size(nex))
2039  col_blk_size = ndo_mo
2040  CALL dbcsr_get_info(matrix_s(1)%matrix, row_blk_size=row_blk_size)
2041 
2042  ALLOCATE (dbcsr_sg, dbcsr_tp, dbcsr_work, dbcsr_ovlp, dbcsr_tmp, dbcsr_prod)
2043  CALL dbcsr_create(matrix=dbcsr_sg, name="SINGLETS", matrix_type=dbcsr_type_no_symmetry, &
2044  dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
2045  CALL dbcsr_create(matrix=dbcsr_tp, name="TRIPLETS", matrix_type=dbcsr_type_no_symmetry, &
2046  dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
2047  CALL dbcsr_create(matrix=dbcsr_work, name="WORK", matrix_type=dbcsr_type_no_symmetry, &
2048  dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
2049  CALL dbcsr_create(matrix=dbcsr_prod, name="PROD", matrix_type=dbcsr_type_no_symmetry, &
2050  dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
2051  CALL dbcsr_create(matrix=dbcsr_ovlp, name="OVLP", matrix_type=dbcsr_type_no_symmetry, &
2052  dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
2053 
2054  col_blk_size = 1
2055  CALL dbcsr_create(matrix=dbcsr_tmp, name="TMP", matrix_type=dbcsr_type_no_symmetry, &
2056  dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
2057  CALL dbcsr_reserve_all_blocks(dbcsr_tmp)
2058 
2059  dbcsr_soc_package%dbcsr_sg => dbcsr_sg
2060  dbcsr_soc_package%dbcsr_tp => dbcsr_tp
2061  dbcsr_soc_package%dbcsr_work => dbcsr_work
2062  dbcsr_soc_package%dbcsr_ovlp => dbcsr_ovlp
2063  dbcsr_soc_package%dbcsr_prod => dbcsr_prod
2064  dbcsr_soc_package%dbcsr_tmp => dbcsr_tmp
2065 
2066  !Filling the coeffs matrices by copying from the stored fms
2067  CALL copy_fm_to_dbcsr(sg_coeffs, dbcsr_sg)
2068  CALL copy_fm_to_dbcsr(tp_coeffs, dbcsr_tp)
2069 
2070 ! Create the work and helper fms
2071  CALL cp_fm_get_info(gs_coeffs, matrix_struct=vec_struct)
2072  CALL cp_fm_struct_create(prod_struct, context=blacs_env, para_env=para_env, &
2073  nrow_global=ndo_mo, ncol_global=ndo_mo)
2074  CALL cp_fm_create(prod_fm, prod_struct)
2075  CALL cp_fm_create(vec_soc_x, vec_struct)
2076  CALL cp_fm_create(vec_soc_y, vec_struct)
2077  CALL cp_fm_create(vec_soc_z, vec_struct)
2078  CALL cp_fm_struct_create(work_struct, context=blacs_env, para_env=para_env, &
2079  nrow_global=nex, ncol_global=nex)
2080  CALL cp_fm_create(work_fm, work_struct)
2081  CALL cp_fm_create(tmp_fm, work_struct)
2082  ALLOCATE (diag(ndo_mo))
2083 
2084 ! Precompute everything we can before looping over excited states
2085  sqrt2 = sqrt(2.0_dp)
2086 
2087  ! The subset of the donor MOs matrix elements: <phi_I^0|Hsoc|phi_J^0> (kept as global array, small)
2088  ALLOCATE (domo_soc_x(ndo_mo, ndo_mo), domo_soc_y(ndo_mo, ndo_mo), domo_soc_z(ndo_mo, ndo_mo))
2089 
2090  CALL cp_dbcsr_sm_fm_multiply(orb_soc_x, gs_coeffs, vec_soc_x, ncol=ndo_mo)
2091  CALL parallel_gemm('T', 'N', ndo_mo, ndo_mo, nao, 1.0_dp, gs_coeffs, vec_soc_x, 0.0_dp, prod_fm)
2092  CALL cp_fm_get_submatrix(prod_fm, domo_soc_x)
2093 
2094  CALL cp_dbcsr_sm_fm_multiply(orb_soc_y, gs_coeffs, vec_soc_y, ncol=ndo_mo)
2095  CALL parallel_gemm('T', 'N', ndo_mo, ndo_mo, nao, 1.0_dp, gs_coeffs, vec_soc_y, 0.0_dp, prod_fm)
2096  CALL cp_fm_get_submatrix(prod_fm, domo_soc_y)
2097 
2098  CALL cp_dbcsr_sm_fm_multiply(orb_soc_z, gs_coeffs, vec_soc_z, ncol=ndo_mo)
2099  CALL parallel_gemm('T', 'N', ndo_mo, ndo_mo, nao, 1.0_dp, gs_coeffs, vec_soc_z, 0.0_dp, prod_fm)
2100  CALL cp_fm_get_submatrix(prod_fm, domo_soc_z)
2101 
2102 ! Only have SOC between singlet-triplet triplet-triplet and ground_state-triplet, the resulting
2103 ! matrix is Hermitian i.e. the real part is symmetric and the imaginary part is anti-symmetric.
2104 ! Can only fill upper half
2105 
2106  !Start with the ground state/triplet SOC, SOC*gs_coeffs already computed above
2107  !note: we are computing <0|H|T>, but have SOC*gs_coeffs instead of gs_coeffs*SOC in store. Since
2108  ! the SOC Hamiltonian is anti-symmetric, a - signs pops up in the gemms below
2109 
2110  CALL cp_fm_struct_create(gstp_struct, context=blacs_env, para_env=para_env, &
2111  nrow_global=ntp*ndo_mo, ncol_global=ndo_mo)
2112  CALL cp_fm_create(gstp_fm, gstp_struct)
2113  ALLOCATE (gstp_block(ndo_mo, ndo_mo))
2114 
2115  !gs-triplet with Ms=+-1, imaginary part
2116  CALL parallel_gemm('T', 'N', ndo_mo*ntp, ndo_mo, nao, -1.0_dp, tp_coeffs, vec_soc_x, 0.0_dp, gstp_fm)
2117 
2118  DO itp = 1, ntp
2119  CALL cp_fm_get_submatrix(fm=gstp_fm, target_m=gstp_block, start_row=(itp - 1)*ndo_mo + 1, &
2120  start_col=1, n_rows=ndo_mo, n_cols=ndo_mo)
2121  diag(:) = get_diag(gstp_block)
2122  soc_gst = sum(diag)
2123  CALL cp_fm_set_element(img_fm, 1, 1 + nsg + itp, -1.0_dp*soc_gst) ! <0|H_x|T^-1>
2124  CALL cp_fm_set_element(img_fm, 1, 1 + nsg + 2*ntp + itp, soc_gst) ! <0|H_x|T^+1>
2125  END DO
2126 
2127  !gs-triplet with Ms=+-1, real part
2128  CALL parallel_gemm('T', 'N', ndo_mo*ntp, ndo_mo, nao, -1.0_dp, tp_coeffs, vec_soc_y, 0.0_dp, gstp_fm)
2129 
2130  DO itp = 1, ntp
2131  CALL cp_fm_get_submatrix(fm=gstp_fm, target_m=gstp_block, start_row=(itp - 1)*ndo_mo + 1, &
2132  start_col=1, n_rows=ndo_mo, n_cols=ndo_mo)
2133  diag(:) = get_diag(gstp_block)
2134  soc_gst = sum(diag)
2135  CALL cp_fm_set_element(real_fm, 1, 1 + nsg + itp, -1.0_dp*soc_gst) ! <0|H_y|T^-1>
2136  CALL cp_fm_set_element(real_fm, 1, 1 + nsg + 2*ntp + itp, -1.0_dp*soc_gst) ! <0|H_y|T^+1>
2137  END DO
2138 
2139  !gs-triplet with Ms=0, purely imaginary
2140  CALL parallel_gemm('T', 'N', ndo_mo*ntp, ndo_mo, nao, -1.0_dp, tp_coeffs, vec_soc_z, 0.0_dp, gstp_fm)
2141 
2142  DO itp = 1, ntp
2143  CALL cp_fm_get_submatrix(fm=gstp_fm, target_m=gstp_block, start_row=(itp - 1)*ndo_mo + 1, &
2144  start_col=1, n_rows=ndo_mo, n_cols=ndo_mo)
2145  diag(:) = get_diag(gstp_block)
2146  soc_gst = sqrt2*sum(diag)
2147  CALL cp_fm_set_element(img_fm, 1, 1 + nsg + ntp + itp, soc_gst)
2148  END DO
2149 
2150  !gs clean-up
2151  CALL cp_fm_release(prod_fm)
2152  CALL cp_fm_release(vec_soc_x)
2153  CALL cp_fm_release(vec_soc_y)
2154  CALL cp_fm_release(vec_soc_z)
2155  CALL cp_fm_release(gstp_fm)
2156  CALL cp_fm_struct_release(gstp_struct)
2157  CALL cp_fm_struct_release(prod_struct)
2158  DEALLOCATE (gstp_block)
2159 
2160  !Now do the singlet-triplet SOC
2161  !start by computing the singlet-triplet overlap
2162  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s(1)%matrix, dbcsr_tp, 0.0_dp, &
2163  dbcsr_work, filter_eps=eps_filter)
2164  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sg, dbcsr_work, 0.0_dp, dbcsr_ovlp, filter_eps=eps_filter)
2165 
2166  !singlet-triplet with Ms=+-1, imaginary part
2167  CALL dbcsr_multiply('N', 'N', 1.0_dp, orb_soc_x, dbcsr_tp, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
2168  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sg, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
2169 
2170  CALL rcs_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_soc_x, &
2171  pref_trace=-1.0_dp, pref_overall=-0.5_dp*sqrt2)
2172 
2173  !<S|H_x|T^-1>
2174  CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
2175  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=img_fm, nrow=nex, ncol=nex, &
2176  s_firstrow=1, s_firstcol=1, t_firstrow=2, &
2177  t_firstcol=1 + nsg + 1)
2178 
2179  !<S|H_x|T^+1> takes a minus sign
2180  CALL cp_fm_scale(-1.0_dp, tmp_fm)
2181  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=img_fm, nrow=nex, ncol=nex, &
2182  s_firstrow=1, s_firstcol=1, t_firstrow=2, &
2183  t_firstcol=1 + nsg + 2*ntp + 1)
2184 
2185  !singlet-triplet with Ms=+-1, real part
2186  CALL dbcsr_multiply('N', 'N', 1.0_dp, orb_soc_y, dbcsr_tp, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
2187  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sg, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
2188 
2189  CALL rcs_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_soc_y, &
2190  pref_trace=-1.0_dp, pref_overall=-0.5_dp*sqrt2)
2191 
2192  !<S|H_y|T^-1>
2193  CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
2194  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=real_fm, nrow=nex, ncol=nex, &
2195  s_firstrow=1, s_firstcol=1, t_firstrow=2, &
2196  t_firstcol=1 + nsg + 1)
2197 
2198  !<S|H_y|T^+1>
2199  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=real_fm, nrow=nex, ncol=nex, &
2200  s_firstrow=1, s_firstcol=1, t_firstrow=2, &
2201  t_firstcol=1 + nsg + 2*ntp + 1)
2202 
2203  !singlet-triplet with Ms=0, purely imaginary
2204  CALL dbcsr_multiply('N', 'N', 1.0_dp, orb_soc_z, dbcsr_tp, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
2205  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sg, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
2206 
2207  CALL rcs_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_soc_z, &
2208  pref_trace=-1.0_dp, pref_overall=1.0_dp)
2209 
2210  !<S|H_z|T^0>
2211  CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
2212  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=img_fm, nrow=nex, ncol=nex, &
2213  s_firstrow=1, s_firstcol=1, t_firstrow=2, &
2214  t_firstcol=1 + nsg + ntp + 1)
2215 
2216  !Now the triplet-triplet SOC
2217  !start by computing the overlap
2218  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s(1)%matrix, dbcsr_tp, 0.0_dp, &
2219  dbcsr_work, filter_eps=eps_filter)
2220  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_tp, dbcsr_work, 0.0_dp, dbcsr_ovlp, filter_eps=eps_filter)
2221 
2222  !Ms=0 to Ms=+-1 SOC, imaginary part
2223  CALL dbcsr_multiply('N', 'N', 1.0_dp, orb_soc_x, dbcsr_tp, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
2224  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_tp, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
2225 
2226  CALL rcs_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_soc_x, &
2227  pref_trace=1.0_dp, pref_overall=-0.5_dp*sqrt2)
2228 
2229  !<T^0|H_x|T^+1>
2230  CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
2231  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=img_fm, nrow=nex, ncol=nex, &
2232  s_firstrow=1, s_firstcol=1, t_firstrow=1 + nsg + ntp + 1, &
2233  t_firstcol=1 + nsg + 2*ntp + 1)
2234 
2235  !<T^-1|H_x|T^0>, takes a minus sign and a transpose (because computed <T^0|H_x|T^-1>)
2236  CALL cp_fm_transpose(tmp_fm, work_fm)
2237  CALL cp_fm_scale(-1.0_dp, work_fm)
2238  CALL cp_fm_to_fm_submat(msource=work_fm, mtarget=img_fm, nrow=nex, ncol=nex, &
2239  s_firstrow=1, s_firstcol=1, t_firstrow=1 + nsg + 1, &
2240  t_firstcol=1 + nsg + ntp + 1)
2241 
2242  !Ms=0 to Ms=+-1 SOC, real part
2243  CALL dbcsr_multiply('N', 'N', 1.0_dp, orb_soc_y, dbcsr_tp, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
2244  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_tp, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
2245 
2246  CALL rcs_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_soc_y, &
2247  pref_trace=1.0_dp, pref_overall=0.5_dp*sqrt2)
2248 
2249  !<T^0|H_y|T^+1>
2250  CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
2251  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=real_fm, nrow=nex, ncol=nex, &
2252  s_firstrow=1, s_firstcol=1, t_firstrow=1 + nsg + ntp + 1, &
2253  t_firstcol=1 + nsg + 2*ntp + 1)
2254 
2255  !<T^-1|H_y|T^0>, takes a minus sign and a transpose
2256  CALL cp_fm_transpose(tmp_fm, work_fm)
2257  CALL cp_fm_scale(-1.0_dp, work_fm)
2258  CALL cp_fm_to_fm_submat(msource=work_fm, mtarget=real_fm, nrow=nex, ncol=nex, &
2259  s_firstrow=1, s_firstcol=1, t_firstrow=1 + nsg + 1, &
2260  t_firstcol=1 + nsg + ntp + 1)
2261 
2262  !Ms=1 to Ms=1 and Ms=-1 to Ms=-1 SOC, purely imaginary
2263  CALL dbcsr_multiply('N', 'N', 1.0_dp, orb_soc_z, dbcsr_tp, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
2264  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_tp, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
2265 
2266  CALL rcs_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_soc_z, &
2267  pref_trace=1.0_dp, pref_overall=1.0_dp)
2268 
2269  !<T^+1|H_z|T^+1>
2270  CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
2271  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=img_fm, nrow=nex, ncol=nex, &
2272  s_firstrow=1, s_firstcol=1, t_firstrow=1 + nsg + 2*ntp + 1, &
2273  t_firstcol=1 + nsg + 2*ntp + 1)
2274 
2275  !<T^-1|H_z|T^-1>, takes a minus sign
2276  CALL cp_fm_scale(-1.0_dp, tmp_fm)
2277  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=img_fm, nrow=nex, ncol=nex, &
2278  s_firstrow=1, s_firstcol=1, t_firstrow=1 + nsg + 1, &
2279  t_firstcol=1 + nsg + 1)
2280 
2281 ! Intermediate clean-up
2282  CALL cp_fm_struct_release(work_struct)
2283  CALL cp_fm_release(work_fm)
2284  CALL cp_fm_release(tmp_fm)
2285  DEALLOCATE (diag, domo_soc_x, domo_soc_y, domo_soc_z)
2286 
2287 ! Set-up the complex hermitian perturbation matrix
2288  CALL cp_cfm_create(hami_cfm, full_struct)
2289  CALL cp_fm_to_cfm(real_fm, img_fm, hami_cfm)
2290 
2291  CALL cp_fm_release(real_fm)
2292  CALL cp_fm_release(img_fm)
2293 
2294 ! Diagonalize the Hamiltonian
2295  ALLOCATE (tmp_evals(ntot))
2296  CALL cp_cfm_create(evecs_cfm, full_struct)
2297  CALL cp_cfm_heevd(hami_cfm, evecs_cfm, tmp_evals)
2298 
2299  ! Adjust the energies so the GS has zero, and store in the donor_state (without the GS)
2300  ALLOCATE (donor_state%soc_evals(ntot - 1))
2301  donor_state%soc_evals(:) = tmp_evals(2:ntot) - tmp_evals(1)
2302 
2303 ! Compute the dipole oscillator strengths
2304  CALL compute_soc_dipole_fosc(evecs_cfm, dbcsr_soc_package, donor_state, xas_tdp_env, &
2305  xas_tdp_control, qs_env)
2306 
2307 ! And the quadrupole (if needed)
2308  IF (xas_tdp_control%do_quad) THEN
2309  CALL compute_soc_quadrupole_fosc(evecs_cfm, dbcsr_soc_package, donor_state, xas_tdp_env, &
2310  xas_tdp_control, qs_env)
2311  END IF
2312 
2313 ! Clean-up
2314  CALL cp_fm_struct_release(full_struct)
2315  CALL cp_cfm_release(hami_cfm)
2316  CALL cp_cfm_release(evecs_cfm)
2317  CALL dbcsr_distribution_release(coeffs_dist)
2318  CALL dbcsr_distribution_release(prod_dist)
2319  CALL dbcsr_release(dbcsr_sg)
2320  CALL dbcsr_release(dbcsr_tp)
2321  CALL dbcsr_release(dbcsr_prod)
2322  CALL dbcsr_release(dbcsr_ovlp)
2323  CALL dbcsr_release(dbcsr_tmp)
2324  CALL dbcsr_release(dbcsr_work)
2325  DEALLOCATE (coeffs_dist, prod_dist, col_dist, col_blk_size, row_dist_new)
2326  DEALLOCATE (dbcsr_sg, dbcsr_tp, dbcsr_work, dbcsr_prod, dbcsr_ovlp, dbcsr_tmp)
2327 
2328  CALL timestop(handle)
2329 
2330  END SUBROUTINE include_rcs_soc
2331 
2332 ! **************************************************************************************************
2333 !> \brief Computes the matrix elements of a one-body operator (given wrt AOs) in the basis of the
2334 !> excited state AMEWs with ground state, for the open-shell case
2335 !> \param amew_op the operator in the basis of the AMEWs (array because could have x,y,z components)
2336 !> \param ao_op the operator in the basis of the atomic orbitals
2337 !> \param gs_coeffs the coefficient of the GS donor MOs. Ecplicitely passed because of special
2338 !> format in the ROKS case (see include_os_soc routine)
2339 !> \param dbcsr_soc_package inhertited from the main SOC routine
2340 !> \param donor_state ...
2341 !> \param eps_filter ...
2342 !> \param qs_env ...
2343 !> \note The ordering of the AMEWs is consistent with SOC and is gs, sc, sf
2344 !> We assume that the operator is spin-independent => only <0|0>, <0|sc>, <sc|sc> and <sf|sf>
2345 !> yield non-zero matrix elements
2346 !> Only for open-shell calculations
2347 ! **************************************************************************************************
2348  SUBROUTINE get_os_amew_op(amew_op, ao_op, gs_coeffs, dbcsr_soc_package, donor_state, &
2349  eps_filter, qs_env)
2350 
2351  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
2352  INTENT(OUT) :: amew_op
2353  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ao_op
2354  TYPE(cp_fm_type), INTENT(IN) :: gs_coeffs
2355  TYPE(dbcsr_soc_package_type) :: dbcsr_soc_package
2356  TYPE(donor_state_type), POINTER :: donor_state
2357  REAL(dp), INTENT(IN) :: eps_filter
2358  TYPE(qs_environment_type), POINTER :: qs_env
2359 
2360  INTEGER :: dim_op, homo, i, isc, nao, ndo_mo, &
2361  ndo_so, nex, nsc, nsf, ntot
2362  REAL(dp) :: op
2363  REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag, gsgs_op
2364  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: domo_op, gsex_block, tmp
2365  TYPE(cp_blacs_env_type), POINTER :: blacs_env
2366  TYPE(cp_fm_struct_type), POINTER :: full_struct, gsex_struct, prod_struct, &
2367  tmp_struct, vec_struct
2368  TYPE(cp_fm_type) :: gsex_fm, prod_work, tmp_fm, vec_work, &
2369  work_fm
2370  TYPE(cp_fm_type), POINTER :: mo_coeff, sc_coeffs, sf_coeffs
2371  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2372  TYPE(dbcsr_type), POINTER :: ao_op_i, dbcsr_ovlp, dbcsr_prod, &
2373  dbcsr_sc, dbcsr_sf, dbcsr_tmp, &
2374  dbcsr_work
2375  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2376  TYPE(mp_para_env_type), POINTER :: para_env
2377 
2378  NULLIFY (matrix_s, para_env, blacs_env, full_struct, vec_struct, prod_struct, mos)
2379  NULLIFY (mo_coeff, ao_op_i, tmp_struct)
2380  NULLIFY (dbcsr_sc, dbcsr_sf, dbcsr_ovlp, dbcsr_work, dbcsr_tmp, dbcsr_prod)
2381 
2382 ! Iinitialization
2383  dim_op = SIZE(ao_op)
2384  sc_coeffs => donor_state%sc_coeffs
2385  sf_coeffs => donor_state%sf_coeffs
2386  nsc = SIZE(donor_state%sc_evals)
2387  nsf = SIZE(donor_state%sf_evals)
2388  nex = nsc
2389  ntot = 1 + nsc + nsf
2390  ndo_mo = donor_state%ndo_mo
2391  ndo_so = 2*donor_state%ndo_mo !open-shell => nspins = 2
2392  CALL get_qs_env(qs_env, matrix_s=matrix_s, para_env=para_env, blacs_env=blacs_env, mos=mos)
2393  CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
2394 
2395  dbcsr_sc => dbcsr_soc_package%dbcsr_sc
2396  dbcsr_sf => dbcsr_soc_package%dbcsr_sf
2397  dbcsr_work => dbcsr_soc_package%dbcsr_work
2398  dbcsr_tmp => dbcsr_soc_package%dbcsr_tmp
2399  dbcsr_prod => dbcsr_soc_package%dbcsr_prod
2400  dbcsr_ovlp => dbcsr_soc_package%dbcsr_ovlp
2401 
2402 ! Create the amew_op matrix set
2403  CALL cp_fm_struct_create(full_struct, context=blacs_env, para_env=para_env, &
2404  nrow_global=ntot, ncol_global=ntot)
2405  ALLOCATE (amew_op(dim_op))
2406  DO i = 1, dim_op
2407  CALL cp_fm_create(amew_op(i), full_struct)
2408  END DO
2409 
2410 ! Before looping, need to evaluate sum_j,sigma <phi^0_j,sgima|op|phi^0_j,sigma>, for each dimension
2411 ! of the operator
2412  ALLOCATE (gsgs_op(dim_op))
2413 
2414  !start with the alpha MOs
2415  CALL get_mo_set(mos(1), mo_coeff=mo_coeff, homo=homo)
2416  ALLOCATE (diag(homo))
2417  CALL cp_fm_get_info(mo_coeff, matrix_struct=vec_struct)
2418  CALL cp_fm_struct_create(prod_struct, context=blacs_env, para_env=para_env, &
2419  nrow_global=homo, ncol_global=homo)
2420  CALL cp_fm_create(vec_work, vec_struct)
2421  CALL cp_fm_create(prod_work, prod_struct)
2422 
2423  DO i = 1, dim_op
2424 
2425  ao_op_i => ao_op(i)%matrix
2426 
2427  CALL cp_dbcsr_sm_fm_multiply(ao_op_i, mo_coeff, vec_work, ncol=homo)
2428  CALL parallel_gemm('T', 'N', homo, homo, nao, 1.0_dp, mo_coeff, vec_work, 0.0_dp, prod_work)
2429  CALL cp_fm_get_diag(prod_work, diag)
2430  gsgs_op(i) = sum(diag)
2431 
2432  END DO !i
2433 
2434  CALL cp_fm_release(vec_work)
2435  CALL cp_fm_release(prod_work)
2436  CALL cp_fm_struct_release(prod_struct)
2437  DEALLOCATE (diag)
2438  NULLIFY (vec_struct)
2439 
2440  !then beta orbitals
2441  CALL get_mo_set(mos(2), mo_coeff=mo_coeff, homo=homo)
2442  ALLOCATE (diag(homo))
2443  CALL cp_fm_get_info(mo_coeff, matrix_struct=vec_struct)
2444  CALL cp_fm_struct_create(prod_struct, context=blacs_env, para_env=para_env, &
2445  nrow_global=homo, ncol_global=homo)
2446  CALL cp_fm_create(vec_work, vec_struct)
2447  CALL cp_fm_create(prod_work, prod_struct)
2448 
2449  DO i = 1, dim_op
2450 
2451  ao_op_i => ao_op(i)%matrix
2452 
2453  CALL cp_dbcsr_sm_fm_multiply(ao_op_i, mo_coeff, vec_work, ncol=homo)
2454  CALL parallel_gemm('T', 'N', homo, homo, nao, 1.0_dp, mo_coeff, vec_work, 0.0_dp, prod_work)
2455  CALL cp_fm_get_diag(prod_work, diag)
2456  gsgs_op(i) = gsgs_op(i) + sum(diag)
2457 
2458  END DO !i
2459 
2460  CALL cp_fm_release(vec_work)
2461  CALL cp_fm_release(prod_work)
2462  CALL cp_fm_struct_release(prod_struct)
2463  DEALLOCATE (diag)
2464  NULLIFY (vec_struct)
2465 
2466 ! Before looping over excited AMEWs, define some work matrices and structures
2467  CALL cp_fm_struct_create(vec_struct, context=blacs_env, para_env=para_env, &
2468  nrow_global=nao, ncol_global=ndo_so)
2469  CALL cp_fm_struct_create(prod_struct, context=blacs_env, para_env=para_env, &
2470  nrow_global=ndo_so, ncol_global=ndo_so)
2471  CALL cp_fm_struct_create(gsex_struct, context=blacs_env, para_env=para_env, &
2472  nrow_global=ndo_so*nex, ncol_global=ndo_so)
2473  CALL cp_fm_struct_create(tmp_struct, context=blacs_env, para_env=para_env, &
2474  nrow_global=nex, ncol_global=nex)
2475  CALL cp_fm_create(vec_work, vec_struct) !for op*|phi>
2476  CALL cp_fm_create(prod_work, prod_struct) !for any <phi|op|phi>
2477  CALL cp_fm_create(work_fm, full_struct)
2478  CALL cp_fm_create(gsex_fm, gsex_struct)
2479  CALL cp_fm_create(tmp_fm, tmp_struct)
2480  ALLOCATE (diag(ndo_so))
2481  ALLOCATE (domo_op(ndo_so, ndo_so))
2482  ALLOCATE (tmp(ndo_so, ndo_so))
2483  ALLOCATE (gsex_block(ndo_so, ndo_so))
2484 
2485 ! Loop over the dimensions of the operator
2486  DO i = 1, dim_op
2487 
2488  ao_op_i => ao_op(i)%matrix
2489 
2490  !put the gs-gs contribution
2491  CALL cp_fm_set_element(amew_op(i), 1, 1, gsgs_op(i))
2492 
2493  ! Precompute what we can before looping over excited states
2494  ! Need the operator in the donor MOs basis <phi^0_I,sigma|op_i|phi^0_J,tau>
2495  CALL cp_dbcsr_sm_fm_multiply(ao_op_i, gs_coeffs, vec_work, ncol=ndo_so)
2496  CALL parallel_gemm('T', 'N', ndo_so, ndo_so, nao, 1.0_dp, gs_coeffs, vec_work, 0.0_dp, prod_work)
2497  CALL cp_fm_get_submatrix(prod_work, domo_op)
2498 
2499  ! Do the ground-state/spin-conserving operator
2500  CALL parallel_gemm('T', 'N', ndo_so*nsc, ndo_so, nao, 1.0_dp, sc_coeffs, vec_work, 0.0_dp, gsex_fm)
2501  DO isc = 1, nsc
2502  CALL cp_fm_get_submatrix(fm=gsex_fm, target_m=gsex_block, start_row=(isc - 1)*ndo_so + 1, &
2503  start_col=1, n_rows=ndo_so, n_cols=ndo_so)
2504  diag(:) = get_diag(gsex_block)
2505  op = sum(diag)
2506  CALL cp_fm_set_element(amew_op(i), 1, 1 + isc, op)
2507  END DO !isc
2508 
2509  ! The spin-conserving/spin-conserving operator
2510  !overlap
2511  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s(1)%matrix, dbcsr_sc, 0.0_dp, &
2512  dbcsr_work, filter_eps=eps_filter)
2513  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sc, dbcsr_work, 0.0_dp, dbcsr_ovlp, filter_eps=eps_filter)
2514 
2515  !operator in SC LR-orbital basis
2516  CALL dbcsr_multiply('N', 'N', 1.0_dp, ao_op_i, dbcsr_sc, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
2517  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sc, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
2518 
2519  CALL os_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_op, pref_diaga=1.0_dp, &
2520  pref_diagb=1.0_dp, pref_tracea=-1.0_dp, pref_traceb=-1.0_dp, &
2521  pref_diags=gsgs_op(i), symmetric=.true.)
2522 
2523  CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
2524  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=amew_op(i), nrow=nex, ncol=nex, &
2525  s_firstrow=1, s_firstcol=1, t_firstrow=2, t_firstcol=2)
2526 
2527  ! The spin-flip/spin-flip operator
2528  !overlap
2529  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s(1)%matrix, dbcsr_sf, 0.0_dp, &
2530  dbcsr_work, filter_eps=eps_filter)
2531  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sf, dbcsr_work, 0.0_dp, dbcsr_ovlp, filter_eps=eps_filter)
2532 
2533  !operator in SF LR-orbital basis
2534  CALL dbcsr_multiply('N', 'N', 1.0_dp, ao_op_i, dbcsr_sf, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
2535  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sf, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
2536 
2537  !need to reorganize the domo_op array by swapping the alpha-alpha and the beta-beta quarter
2538  tmp(1:ndo_mo, 1:ndo_mo) = domo_op(ndo_mo + 1:ndo_so, ndo_mo + 1:ndo_so)
2539  tmp(ndo_mo + 1:ndo_so, ndo_mo + 1:ndo_so) = domo_op(1:ndo_mo, 1:ndo_mo)
2540 
2541  CALL os_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, tmp, pref_diaga=1.0_dp, &
2542  pref_diagb=1.0_dp, pref_tracea=-1.0_dp, pref_traceb=-1.0_dp, &
2543  pref_diags=gsgs_op(i), symmetric=.true.)
2544 
2545  CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
2546  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=amew_op(i), nrow=nex, ncol=nex, &
2547  s_firstrow=1, s_firstcol=1, t_firstrow=1 + nsc + 1, t_firstcol=1 + nsc + 1)
2548 
2549  !Symmetry => only upper diag explicitly built
2550  CALL cp_fm_upper_to_full(amew_op(i), work_fm)
2551 
2552  END DO !i
2553 
2554 ! Clean-up
2555  CALL cp_fm_struct_release(full_struct)
2556  CALL cp_fm_struct_release(prod_struct)
2557  CALL cp_fm_struct_release(vec_struct)
2558  CALL cp_fm_struct_release(tmp_struct)
2559  CALL cp_fm_struct_release(gsex_struct)
2560  CALL cp_fm_release(work_fm)
2561  CALL cp_fm_release(tmp_fm)
2562  CALL cp_fm_release(vec_work)
2563  CALL cp_fm_release(prod_work)
2564  CALL cp_fm_release(gsex_fm)
2565 
2566  END SUBROUTINE get_os_amew_op
2567 
2568 ! **************************************************************************************************
2569 !> \brief Computes the matrix elements of a one-body operator (given wrt AOs) in the basis of the
2570 !> excited state AMEWs with ground state, singlet and triplet with Ms = -1,0,+1
2571 !> \param amew_op the operator in the basis of the AMEWs (array because could have x,y,z components)
2572 !> \param ao_op the operator in the basis of the atomic orbitals
2573 !> \param dbcsr_soc_package inherited from the main SOC routine
2574 !> \param donor_state ...
2575 !> \param eps_filter for dbcsr multiplication
2576 !> \param qs_env ...
2577 !> \note The ordering of the AMEWs is consistent with SOC and is gs, sg, tp(-1), tp(0). tp(+1)
2578 !> We assume that the operator is spin-independent => only <0|0>, <0|S>, <S|S> and <T|T>
2579 !> yield non-zero matrix elements
2580 !> Only for spin-restricted calculations
2581 ! **************************************************************************************************
2582  SUBROUTINE get_rcs_amew_op(amew_op, ao_op, dbcsr_soc_package, donor_state, eps_filter, qs_env)
2583 
2584  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
2585  INTENT(OUT) :: amew_op
2586  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ao_op
2587  TYPE(dbcsr_soc_package_type) :: dbcsr_soc_package
2588  TYPE(donor_state_type), POINTER :: donor_state
2589  REAL(dp), INTENT(IN) :: eps_filter
2590  TYPE(qs_environment_type), POINTER :: qs_env
2591 
2592  INTEGER :: dim_op, homo, i, isg, nao, ndo_mo, nex, &
2593  nsg, ntot, ntp
2594  REAL(dp) :: op, sqrt2
2595  REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag, gs_diag, gsgs_op
2596  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: domo_op, sggs_block
2597  TYPE(cp_blacs_env_type), POINTER :: blacs_env
2598  TYPE(cp_fm_struct_type), POINTER :: full_struct, gsgs_struct, prod_struct, &
2599  sggs_struct, std_struct, tmp_struct, &
2600  vec_struct
2601  TYPE(cp_fm_type) :: gs_fm, prod_fm, sggs_fm, tmp_fm, vec_op, &
2602  work_fm
2603  TYPE(cp_fm_type), POINTER :: gs_coeffs, mo_coeff, sg_coeffs
2604  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
2605  TYPE(dbcsr_type), POINTER :: ao_op_i, dbcsr_ovlp, dbcsr_prod, &
2606  dbcsr_sg, dbcsr_tmp, dbcsr_tp, &
2607  dbcsr_work
2608  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2609  TYPE(mp_para_env_type), POINTER :: para_env
2610 
2611  NULLIFY (gs_coeffs, sg_coeffs, matrix_s, full_struct, prod_struct, vec_struct, blacs_env)
2612  NULLIFY (para_env, mo_coeff, mos, gsgs_struct, std_struct, tmp_struct, sggs_struct)
2613  NULLIFY (ao_op_i, dbcsr_tp, dbcsr_sg, dbcsr_ovlp, dbcsr_work, dbcsr_tmp, dbcsr_prod)
2614 
2615 ! Initialization
2616  gs_coeffs => donor_state%gs_coeffs
2617  sg_coeffs => donor_state%sg_coeffs
2618  nsg = SIZE(donor_state%sg_evals)
2619  ntp = nsg; nex = nsg !all the same by construction, keep them separate for clarity
2620  ntot = 1 + nsg + 3*ntp
2621  ndo_mo = donor_state%ndo_mo
2622  CALL get_qs_env(qs_env, matrix_s=matrix_s, para_env=para_env, blacs_env=blacs_env, mos=mos)
2623  sqrt2 = sqrt(2.0_dp)
2624  dim_op = SIZE(ao_op)
2625 
2626  dbcsr_sg => dbcsr_soc_package%dbcsr_sg
2627  dbcsr_tp => dbcsr_soc_package%dbcsr_tp
2628  dbcsr_work => dbcsr_soc_package%dbcsr_work
2629  dbcsr_prod => dbcsr_soc_package%dbcsr_prod
2630  dbcsr_ovlp => dbcsr_soc_package%dbcsr_ovlp
2631  dbcsr_tmp => dbcsr_soc_package%dbcsr_tmp
2632 
2633 ! Create the amew_op matrix
2634  CALL cp_fm_struct_create(full_struct, context=blacs_env, para_env=para_env, &
2635  nrow_global=ntot, ncol_global=ntot)
2636  ALLOCATE (amew_op(dim_op))
2637  DO i = 1, dim_op
2638  CALL cp_fm_create(amew_op(i), full_struct)
2639  END DO !i
2640 
2641 ! Deal with the GS-GS contribution <0|0> = 2*sum_j <phi_j|op|phi_j>
2642  CALL get_mo_set(mos(1), mo_coeff=mo_coeff, nao=nao, homo=homo)
2643  CALL cp_fm_struct_create(gsgs_struct, context=blacs_env, para_env=para_env, &
2644  nrow_global=homo, ncol_global=homo)
2645  CALL cp_fm_get_info(mo_coeff, matrix_struct=std_struct)
2646  CALL cp_fm_create(gs_fm, gsgs_struct)
2647  CALL cp_fm_create(work_fm, std_struct)
2648  ALLOCATE (gsgs_op(dim_op))
2649  ALLOCATE (gs_diag(homo))
2650 
2651  DO i = 1, dim_op
2652 
2653  ao_op_i => ao_op(i)%matrix
2654 
2655  CALL cp_dbcsr_sm_fm_multiply(ao_op_i, mo_coeff, work_fm, ncol=homo)
2656  CALL parallel_gemm('T', 'N', homo, homo, nao, 1.0_dp, mo_coeff, work_fm, 0.0_dp, gs_fm)
2657  CALL cp_fm_get_diag(gs_fm, gs_diag)
2658  gsgs_op(i) = 2.0_dp*sum(gs_diag)
2659 
2660  END DO !i
2661 
2662  CALL cp_fm_release(gs_fm)
2663  CALL cp_fm_release(work_fm)
2664  CALL cp_fm_struct_release(gsgs_struct)
2665  DEALLOCATE (gs_diag)
2666 
2667 ! Create the work and helper fms
2668  CALL cp_fm_get_info(gs_coeffs, matrix_struct=vec_struct)
2669  CALL cp_fm_struct_create(prod_struct, context=blacs_env, para_env=para_env, &
2670  nrow_global=ndo_mo, ncol_global=ndo_mo)
2671  CALL cp_fm_create(prod_fm, prod_struct)
2672  CALL cp_fm_create(vec_op, vec_struct)
2673  CALL cp_fm_struct_create(tmp_struct, context=blacs_env, para_env=para_env, &
2674  nrow_global=nex, ncol_global=nex)
2675  CALL cp_fm_struct_create(sggs_struct, context=blacs_env, para_env=para_env, &
2676  nrow_global=ndo_mo*nsg, ncol_global=ndo_mo)
2677  CALL cp_fm_create(tmp_fm, tmp_struct)
2678  CALL cp_fm_create(work_fm, full_struct)
2679  CALL cp_fm_create(sggs_fm, sggs_struct)
2680  ALLOCATE (diag(ndo_mo))
2681  ALLOCATE (domo_op(ndo_mo, ndo_mo))
2682  ALLOCATE (sggs_block(ndo_mo, ndo_mo))
2683 
2684 ! Iterate over the dimensions of the operator
2685 ! Note: operator matrices are asusmed symmetric, can only do upper half
2686  DO i = 1, dim_op
2687 
2688  ao_op_i => ao_op(i)%matrix
2689 
2690  ! The GS-GS contribution
2691  CALL cp_fm_set_element(amew_op(i), 1, 1, gsgs_op(i))
2692 
2693  ! Compute the operator in the donor MOs basis
2694  CALL cp_dbcsr_sm_fm_multiply(ao_op_i, gs_coeffs, vec_op, ncol=ndo_mo)
2695  CALL parallel_gemm('T', 'N', ndo_mo, ndo_mo, nao, 1.0_dp, gs_coeffs, vec_op, 0.0_dp, prod_fm)
2696  CALL cp_fm_get_submatrix(prod_fm, domo_op)
2697 
2698  ! Compute the ground-state/singlet components. ao_op*gs_coeffs already stored in vec_op
2699  CALL parallel_gemm('T', 'N', ndo_mo*nsg, ndo_mo, nao, 1.0_dp, sg_coeffs, vec_op, 0.0_dp, sggs_fm)
2700  DO isg = 1, nsg
2701  CALL cp_fm_get_submatrix(fm=sggs_fm, target_m=sggs_block, start_row=(isg - 1)*ndo_mo + 1, &
2702  start_col=1, n_rows=ndo_mo, n_cols=ndo_mo)
2703  diag(:) = get_diag(sggs_block)
2704  op = sqrt2*sum(diag)
2705  CALL cp_fm_set_element(amew_op(i), 1, 1 + isg, op)
2706  END DO
2707 
2708  ! do the singlet-singlet components
2709  !start with the overlap
2710  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s(1)%matrix, dbcsr_sg, 0.0_dp, &
2711  dbcsr_work, filter_eps=eps_filter)
2712  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sg, dbcsr_work, 0.0_dp, dbcsr_ovlp, filter_eps=eps_filter)
2713 
2714  !then the operator in the LR orbital basis
2715  CALL dbcsr_multiply('N', 'N', 1.0_dp, ao_op_i, dbcsr_sg, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
2716  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sg, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
2717 
2718  !use the soc routine, it is compatible
2719  CALL rcs_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_op, pref_trace=-1.0_dp, &
2720  pref_overall=1.0_dp, pref_diags=gsgs_op(i), symmetric=.true.)
2721 
2722  CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
2723  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=amew_op(i), nrow=nex, ncol=nex, &
2724  s_firstrow=1, s_firstcol=1, t_firstrow=2, t_firstcol=2)
2725 
2726  ! compute the triplet-triplet components
2727  !the overlap
2728  CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s(1)%matrix, dbcsr_tp, 0.0_dp, &
2729  dbcsr_work, filter_eps=eps_filter)
2730  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_tp, dbcsr_work, 0.0_dp, dbcsr_ovlp, filter_eps=eps_filter)
2731 
2732  !the operator in the LR orbital basis
2733  CALL dbcsr_multiply('N', 'N', 1.0_dp, ao_op_i, dbcsr_sg, 0.0_dp, dbcsr_work, filter_eps=eps_filter)
2734  CALL dbcsr_multiply('T', 'N', 1.0_dp, dbcsr_sg, dbcsr_work, 0.0_dp, dbcsr_prod, filter_eps=eps_filter)
2735 
2736  CALL rcs_amew_soc_elements(dbcsr_tmp, dbcsr_prod, dbcsr_ovlp, domo_op, pref_trace=-1.0_dp, &
2737  pref_overall=1.0_dp, pref_diags=gsgs_op(i), symmetric=.true.)
2738 
2739  CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
2740  !<T^-1|op|T^-1>
2741  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=amew_op(i), nrow=nex, ncol=nex, &
2742  s_firstrow=1, s_firstcol=1, t_firstrow=1 + nsg + 1, t_firstcol=1 + nsg + 1)
2743  !<T^0|op|T^0>
2744  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=amew_op(i), nrow=nex, ncol=nex, &
2745  s_firstrow=1, s_firstcol=1, t_firstrow=1 + nsg + ntp + 1, &
2746  t_firstcol=1 + nsg + ntp + 1)
2747  !<T^-1|op|T^-1>
2748  CALL cp_fm_to_fm_submat(msource=tmp_fm, mtarget=amew_op(i), nrow=nex, ncol=nex, &
2749  s_firstrow=1, s_firstcol=1, t_firstrow=1 + nsg + 2*ntp + 1, &
2750  t_firstcol=1 + nsg + 2*ntp + 1)
2751 
2752  ! Symmetrize the matrix (only upper triangle built)
2753  CALL cp_fm_upper_to_full(amew_op(i), work_fm)
2754 
2755  END DO !i
2756 
2757 ! Clean-up
2758  CALL cp_fm_release(prod_fm)
2759  CALL cp_fm_release(work_fm)
2760  CALL cp_fm_release(tmp_fm)
2761  CALL cp_fm_release(vec_op)
2762  CALL cp_fm_release(sggs_fm)
2763  CALL cp_fm_struct_release(prod_struct)
2764  CALL cp_fm_struct_release(full_struct)
2765  CALL cp_fm_struct_release(tmp_struct)
2766  CALL cp_fm_struct_release(sggs_struct)
2767 
2768  END SUBROUTINE get_rcs_amew_op
2769 
2770 ! **************************************************************************************************
2771 !> \brief Computes the os SOC matrix elements between excited states AMEWs based on the LR orbitals
2772 !> \param amew_soc output dbcsr matrix with the SOC in the AMEW basis (needs to be fully resereved)
2773 !> \param lr_soc dbcsr matrix with the SOC wrt the LR orbitals
2774 !> \param lr_overlap dbcsr matrix with the excited states LR orbital overlap
2775 !> \param domo_soc the SOC in the basis of the donor MOs
2776 !> \param pref_diaga ...
2777 !> \param pref_diagb ...
2778 !> \param pref_tracea ...
2779 !> \param pref_traceb ...
2780 !> \param pref_diags see notes
2781 !> \param symmetric if the outcome is known to be symmetric, only elements with iex <= jex are done
2782 !> \param tracea_start the indices where to start in the trace part for alpha
2783 !> \param traceb_start the indices where to start in the trace part for beta
2784 !> \note For an excited states pair i,j, the AMEW SOC matrix element is:
2785 !> soc_ij = pref_diaga*SUM(alpha part of diag of lr_soc_ij)
2786 !> + pref_diagb*SUM(beta part of diag of lr_soc_ij)
2787 !> + pref_tracea*SUM(alpha part of lr_ovlp_ij*TRANSPOSE(domo_soc))
2788 !> + pref_traceb*SUM(beta part of lr_ovlp_ij*TRANSPOSE(domo_soc))
2789 !> optinally, one can add pref_diags*SUM(diag lr_ovlp_ij)
2790 ! **************************************************************************************************
2791  SUBROUTINE os_amew_soc_elements(amew_soc, lr_soc, lr_overlap, domo_soc, pref_diaga, &
2792  pref_diagb, pref_tracea, pref_traceb, pref_diags, &
2793  symmetric, tracea_start, traceb_start)
2794 
2795  TYPE(dbcsr_type) :: amew_soc, lr_soc, lr_overlap
2796  REAL(dp), DIMENSION(:, :) :: domo_soc
2797  REAL(dp) :: pref_diaga, pref_diagb, pref_tracea, &
2798  pref_traceb
2799  REAL(dp), OPTIONAL :: pref_diags
2800  LOGICAL, OPTIONAL :: symmetric
2801  INTEGER, DIMENSION(2), OPTIONAL :: tracea_start, traceb_start
2802 
2803  INTEGER :: blk, iex, jex, ndo_mo, ndo_so
2804  INTEGER, DIMENSION(2) :: tas, tbs
2805  LOGICAL :: do_diags, found, my_symm
2806  REAL(dp) :: soc_elem
2807  REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag
2808  REAL(dp), DIMENSION(:, :), POINTER :: pblock
2809  TYPE(dbcsr_iterator_type) :: iter
2810 
2811  ndo_so = SIZE(domo_soc, 1)
2812  ndo_mo = ndo_so/2
2813  ALLOCATE (diag(ndo_so))
2814  my_symm = .false.
2815  IF (PRESENT(symmetric)) my_symm = symmetric
2816  do_diags = .false.
2817  IF (PRESENT(pref_diags)) do_diags = .true.
2818 
2819  !by default, alpha part is (1:ndo_mo,1:ndo_mo) and beta is (ndo_mo+1:ndo_so,ndo_mo+1:ndo_so)
2820  !note: in some SF cases, that might change, mainly because the spin-flip LR-coeffs have
2821  !inverse order, that is: the beta-coeffs in the alpha spot and the alpha coeffs in the
2822  !beta spot
2823  tas = 1
2824  tbs = ndo_mo + 1
2825  IF (PRESENT(tracea_start)) tas = tracea_start
2826  IF (PRESENT(traceb_start)) tbs = traceb_start
2827 
2828  CALL dbcsr_set(amew_soc, 0.0_dp)
2829  !loop over the excited states pairs as the block of amew_soc (which are all reserved)
2830  CALL dbcsr_iterator_start(iter, amew_soc)
2831  DO WHILE (dbcsr_iterator_blocks_left(iter))
2832 
2833  CALL dbcsr_iterator_next_block(iter, row=iex, column=jex, blk=blk)
2834 
2835  IF (my_symm .AND. iex > jex) cycle
2836 
2837  !compute the soc matrix element
2838  soc_elem = 0.0_dp
2839  CALL dbcsr_get_block_p(lr_soc, iex, jex, pblock, found)
2840  IF (found) THEN
2841  diag(:) = get_diag(pblock)
2842  soc_elem = soc_elem + pref_diaga*sum(diag(1:ndo_mo)) + pref_diagb*(sum(diag(ndo_mo + 1:ndo_so)))
2843  END IF
2844 
2845  CALL dbcsr_get_block_p(lr_overlap, iex, jex, pblock, found)
2846  IF (found) THEN
2847  soc_elem = soc_elem &
2848  + pref_tracea*sum(pblock(tas(1):tas(1) + ndo_mo - 1, tas(2):tas(2) + ndo_mo - 1)* &
2849  domo_soc(tas(1):tas(1) + ndo_mo - 1, tas(2):tas(2) + ndo_mo - 1)) &
2850  + pref_traceb*sum(pblock(tbs(1):tbs(1) + ndo_mo - 1, tbs(2):tbs(2) + ndo_mo - 1)* &
2851  domo_soc(tbs(1):tbs(1) + ndo_mo - 1, tbs(2):tbs(2) + ndo_mo - 1))
2852 
2853  IF (do_diags) THEN
2854  diag(:) = get_diag(pblock)
2855  soc_elem = soc_elem + pref_diags*sum(diag)
2856  END IF
2857  END IF
2858 
2859  CALL dbcsr_get_block_p(amew_soc, iex, jex, pblock, found)
2860  pblock = soc_elem
2861 
2862  END DO
2863  CALL dbcsr_iterator_stop(iter)
2864 
2865  END SUBROUTINE os_amew_soc_elements
2866 
2867 ! **************************************************************************************************
2868 !> \brief Computes the rcs SOC matrix elements between excited states AMEWs based on the LR orbitals
2869 !> \param amew_soc output dbcsr matrix with the SOC in the AMEW basis (needs to be fully resereved)
2870 !> \param lr_soc dbcsr matrix with the SOC wrt the LR orbitals
2871 !> \param lr_overlap dbcsr matrix with the excited states LR orbital overlap
2872 !> \param domo_soc the SOC in the basis of the donor MOs
2873 !> \param pref_trace see notes
2874 !> \param pref_overall see notes
2875 !> \param pref_diags see notes
2876 !> \param symmetric if the outcome is known to be symmetric, only elements with iex <= jex are done
2877 !> \note For an excited states pair i,j, the AMEW SOC matrix element is:
2878 !> soc_ij = pref_overall*(SUM(diag(lr_soc_ij)) + pref_trace*SUM(lr_overlap_ij*TRANSPOSE(domo_soc)))
2879 !> optionally, the value pref_diags*SUM(diag(lr_overlap_ij)) can be added (before pref_overall)
2880 ! **************************************************************************************************
2881  SUBROUTINE rcs_amew_soc_elements(amew_soc, lr_soc, lr_overlap, domo_soc, pref_trace, &
2882  pref_overall, pref_diags, symmetric)
2883 
2884  TYPE(dbcsr_type) :: amew_soc, lr_soc, lr_overlap
2885  REAL(dp), DIMENSION(:, :) :: domo_soc
2886  REAL(dp) :: pref_trace, pref_overall
2887  REAL(dp), OPTIONAL :: pref_diags
2888  LOGICAL, OPTIONAL :: symmetric
2889 
2890  INTEGER :: blk, iex, jex
2891  LOGICAL :: do_diags, found, my_symm
2892  REAL(dp) :: soc_elem
2893  REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag
2894  REAL(dp), DIMENSION(:, :), POINTER :: pblock
2895  TYPE(dbcsr_iterator_type) :: iter
2896 
2897  ALLOCATE (diag(SIZE(domo_soc, 1)))
2898  my_symm = .false.
2899  IF (PRESENT(symmetric)) my_symm = symmetric
2900  do_diags = .false.
2901  IF (PRESENT(pref_diags)) do_diags = .true.
2902 
2903  CALL dbcsr_set(amew_soc, 0.0_dp)
2904  !loop over the excited states pairs as the block of amew_soc (which are all reserved)
2905  CALL dbcsr_iterator_start(iter, amew_soc)
2906  DO WHILE (dbcsr_iterator_blocks_left(iter))
2907 
2908  CALL dbcsr_iterator_next_block(iter, row=iex, column=jex, blk=blk)
2909 
2910  IF (my_symm .AND. iex > jex) cycle
2911 
2912  !compute the soc matrix element
2913  soc_elem = 0.0_dp
2914  CALL dbcsr_get_block_p(lr_soc, iex, jex, pblock, found)
2915  IF (found) THEN
2916  diag(:) = get_diag(pblock)
2917  soc_elem = soc_elem + sum(diag)
2918  END IF
2919 
2920  CALL dbcsr_get_block_p(lr_overlap, iex, jex, pblock, found)
2921  IF (found) THEN
2922  soc_elem = soc_elem + pref_trace*sum(pblock*transpose(domo_soc))
2923 
2924  IF (do_diags) THEN
2925  diag(:) = get_diag(pblock)
2926  soc_elem = soc_elem + pref_diags*sum(diag)
2927  END IF
2928  END IF
2929 
2930  CALL dbcsr_get_block_p(amew_soc, iex, jex, pblock, found)
2931  pblock = pref_overall*soc_elem
2932 
2933  END DO
2934  CALL dbcsr_iterator_stop(iter)
2935 
2936  END SUBROUTINE rcs_amew_soc_elements
2937 
2938 ! **************************************************************************************************
2939 !> \brief Computes the dipole oscillator strengths in the AMEWs basis for SOC
2940 !> \param soc_evecs_cfm the complex AMEWs coefficients
2941 !> \param dbcsr_soc_package ...
2942 !> \param donor_state ...
2943 !> \param xas_tdp_env ...
2944 !> \param xas_tdp_control ...
2945 !> \param qs_env ...
2946 !> \param gs_coeffs the ground state coefficients, given for open-shell because in ROKS, the gs_coeffs
2947 !> are stored slightly differently within SOC for efficiency and code uniquness
2948 ! **************************************************************************************************
2949  SUBROUTINE compute_soc_dipole_fosc(soc_evecs_cfm, dbcsr_soc_package, donor_state, xas_tdp_env, &
2950  xas_tdp_control, qs_env, gs_coeffs)
2951 
2952  TYPE(cp_cfm_type), INTENT(IN) :: soc_evecs_cfm
2953  TYPE(dbcsr_soc_package_type) :: dbcsr_soc_package
2954  TYPE(donor_state_type), POINTER :: donor_state
2955  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
2956  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
2957  TYPE(qs_environment_type), POINTER :: qs_env
2958  TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: gs_coeffs
2959 
2960  CHARACTER(len=*), PARAMETER :: routinen = 'compute_soc_dipole_fosc'
2961 
2962  COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: transdip
2963  INTEGER :: handle, i, nosc, ntot
2964  LOGICAL :: do_os, do_rcs
2965  REAL(dp), ALLOCATABLE, DIMENSION(:) :: osc_xyz
2966  REAL(dp), DIMENSION(:), POINTER :: soc_evals
2967  REAL(dp), DIMENSION(:, :), POINTER :: osc_str
2968  TYPE(cp_blacs_env_type), POINTER :: blacs_env
2969  TYPE(cp_cfm_type) :: dip_cfm, work1_cfm, work2_cfm
2970  TYPE(cp_fm_struct_type), POINTER :: dip_struct, full_struct
2971  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: amew_dip
2972  TYPE(mp_para_env_type), POINTER :: para_env
2973 
2974  NULLIFY (para_env, blacs_env, dip_struct, full_struct, osc_str)
2975  NULLIFY (soc_evals)
2976 
2977  CALL timeset(routinen, handle)
2978 
2979  !init
2980  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
2981  do_os = xas_tdp_control%do_spin_cons
2982  do_rcs = xas_tdp_control%do_singlet
2983  soc_evals => donor_state%soc_evals
2984  nosc = SIZE(soc_evals)
2985  ntot = nosc + 1 !because GS AMEW is in there
2986  ALLOCATE (donor_state%soc_osc_str(nosc, 4))
2987  osc_str => donor_state%soc_osc_str
2988  osc_str(:, :) = 0.0_dp
2989  IF (do_os .AND. .NOT. PRESENT(gs_coeffs)) cpabort("Need to pass gs_coeffs for open-shell")
2990 
2991  !get some work arrays/matrix
2992  CALL cp_fm_struct_create(dip_struct, context=blacs_env, para_env=para_env, &
2993  nrow_global=ntot, ncol_global=1)
2994  CALL cp_cfm_get_info(soc_evecs_cfm, matrix_struct=full_struct)
2995  CALL cp_cfm_create(dip_cfm, dip_struct)
2996  CALL cp_cfm_create(work1_cfm, full_struct)
2997  CALL cp_cfm_create(work2_cfm, full_struct)
2998  ALLOCATE (transdip(ntot, 1))
2999 
3000  !get the dipole in the AMEW basis
3001  IF (do_os) THEN
3002  CALL get_os_amew_op(amew_dip, xas_tdp_env%dipmat, gs_coeffs, dbcsr_soc_package, &
3003  donor_state, xas_tdp_control%eps_filter, qs_env)
3004  ELSE
3005  CALL get_rcs_amew_op(amew_dip, xas_tdp_env%dipmat, dbcsr_soc_package, donor_state, &
3006  xas_tdp_control%eps_filter, qs_env)
3007  END IF
3008 
3009  ALLOCATE (osc_xyz(nosc))
3010  DO i = 1, 3 !cartesian coord x, y, z
3011 
3012  !Convert the real dipole into the cfm format for calculations
3013  CALL cp_fm_to_cfm(msourcer=amew_dip(i), mtarget=work1_cfm)
3014 
3015  !compute amew_coeffs^dagger * amew_dip * amew_gs to get the transition moments
3016  CALL parallel_gemm('C', 'N', ntot, ntot, ntot, (1.0_dp, 0.0_dp), soc_evecs_cfm, work1_cfm, &
3017  (0.0_dp, 0.0_dp), work2_cfm)
3018  CALL parallel_gemm('N', 'N', ntot, 1, ntot, (1.0_dp, 0.0_dp), work2_cfm, soc_evecs_cfm, &
3019  (0.0_dp, 0.0_dp), dip_cfm)
3020 
3021  CALL cp_cfm_get_submatrix(dip_cfm, transdip)
3022 
3023  !transition dipoles are real numbers
3024  osc_xyz(:) = real(transdip(2:ntot, 1))**2 + aimag(transdip(2:ntot, 1))**2
3025  osc_str(:, 4) = osc_str(:, 4) + osc_xyz(:)
3026  osc_str(:, i) = osc_xyz(:)
3027 
3028  END DO !i
3029 
3030  !multiply with appropriate prefac depending in the rep
3031  DO i = 1, 4
3032  IF (xas_tdp_control%dipole_form == xas_dip_len) THEN
3033  osc_str(:, i) = 2.0_dp/3.0_dp*soc_evals(:)*osc_str(:, i)
3034  ELSE
3035  osc_str(:, i) = 2.0_dp/3.0_dp/soc_evals(:)*osc_str(:, i)
3036  END IF
3037  END DO
3038 
3039  !clean-up
3040  CALL cp_fm_struct_release(dip_struct)
3041  CALL cp_cfm_release(work1_cfm)
3042  CALL cp_cfm_release(work2_cfm)
3043  CALL cp_cfm_release(dip_cfm)
3044  DO i = 1, 3
3045  CALL cp_fm_release(amew_dip(i))
3046  END DO
3047  DEALLOCATE (amew_dip, transdip)
3048 
3049  CALL timestop(handle)
3050 
3051  END SUBROUTINE compute_soc_dipole_fosc
3052 
3053 ! **************************************************************************************************
3054 !> \brief Computes the quadrupole oscillator strengths in the AMEWs basis for SOC
3055 !> \param soc_evecs_cfm the complex AMEWs coefficients
3056 !> \param dbcsr_soc_package inherited from the main SOC routine
3057 !> \param donor_state ...
3058 !> \param xas_tdp_env ...
3059 !> \param xas_tdp_control ...
3060 !> \param qs_env ...
3061 !> \param gs_coeffs the ground state coefficients, given for open-shell because in ROKS, the gs_coeffs
3062 !> are stored slightly differently within SOC for efficiency and code uniquness
3063 ! **************************************************************************************************
3064  SUBROUTINE compute_soc_quadrupole_fosc(soc_evecs_cfm, dbcsr_soc_package, donor_state, &
3065  xas_tdp_env, xas_tdp_control, qs_env, gs_coeffs)
3066 
3067  TYPE(cp_cfm_type), INTENT(IN) :: soc_evecs_cfm
3068  TYPE(dbcsr_soc_package_type) :: dbcsr_soc_package
3069  TYPE(donor_state_type), POINTER :: donor_state
3070  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
3071  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
3072  TYPE(qs_environment_type), POINTER :: qs_env
3073  TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: gs_coeffs
3074 
3075  CHARACTER(len=*), PARAMETER :: routinen = 'compute_soc_quadrupole_fosc'
3076 
3077  COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: trace
3078  COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: transquad
3079  INTEGER :: handle, i, nosc, ntot
3080  LOGICAL :: do_os, do_rcs
3081  REAL(dp), DIMENSION(:), POINTER :: osc_str, soc_evals
3082  TYPE(cp_blacs_env_type), POINTER :: blacs_env
3083  TYPE(cp_cfm_type) :: quad_cfm, work1_cfm, work2_cfm
3084  TYPE(cp_fm_struct_type), POINTER :: full_struct, quad_struct
3085  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: amew_quad
3086  TYPE(mp_para_env_type), POINTER :: para_env
3087 
3088  NULLIFY (para_env, blacs_env, quad_struct, full_struct, osc_str)
3089  NULLIFY (soc_evals)
3090 
3091  CALL timeset(routinen, handle)
3092 
3093  !init
3094  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
3095  do_os = xas_tdp_control%do_spin_cons
3096  do_rcs = xas_tdp_control%do_singlet
3097  soc_evals => donor_state%soc_evals
3098  nosc = SIZE(soc_evals)
3099  ntot = nosc + 1 !because GS AMEW is in there
3100  ALLOCATE (donor_state%soc_quad_osc_str(nosc))
3101  osc_str => donor_state%soc_quad_osc_str
3102  osc_str(:) = 0.0_dp
3103  IF (do_os .AND. .NOT. PRESENT(gs_coeffs)) cpabort("Need to pass gs_coeffs for open-shell")
3104 
3105  !get some work arrays/matrix
3106  CALL cp_fm_struct_create(quad_struct, context=blacs_env, para_env=para_env, &
3107  nrow_global=ntot, ncol_global=1)
3108  CALL cp_cfm_get_info(soc_evecs_cfm, matrix_struct=full_struct)
3109  CALL cp_cfm_create(quad_cfm, quad_struct)
3110  CALL cp_cfm_create(work1_cfm, full_struct)
3111  CALL cp_cfm_create(work2_cfm, full_struct)
3112  ALLOCATE (transquad(ntot, 1))
3113  ALLOCATE (trace(nosc))
3114  trace = (0.0_dp, 0.0_dp)
3115 
3116  !get the quadrupole in the AMEWs basis
3117  IF (do_os) THEN
3118  CALL get_os_amew_op(amew_quad, xas_tdp_env%quadmat, gs_coeffs, dbcsr_soc_package, &
3119  donor_state, xas_tdp_control%eps_filter, qs_env)
3120  ELSE
3121  CALL get_rcs_amew_op(amew_quad, xas_tdp_env%quadmat, dbcsr_soc_package, donor_state, &
3122  xas_tdp_control%eps_filter, qs_env)
3123  END IF
3124 
3125  DO i = 1, 6 ! x2, xy, xz, y2, yz, z2
3126 
3127  !Convert the real quadrupole into a cfm for further calculation
3128  CALL cp_fm_to_cfm(msourcer=amew_quad(i), mtarget=work1_cfm)
3129 
3130  !compute amew_coeffs^dagger * amew_quad * amew_gs to get the transition moments
3131  CALL parallel_gemm('C', 'N', ntot, ntot, ntot, (1.0_dp, 0.0_dp), soc_evecs_cfm, work1_cfm, &
3132  (0.0_dp, 0.0_dp), work2_cfm)
3133  CALL parallel_gemm('N', 'N', ntot, 1, ntot, (1.0_dp, 0.0_dp), work2_cfm, soc_evecs_cfm, &
3134  (0.0_dp, 0.0_dp), quad_cfm)
3135 
3136  CALL cp_cfm_get_submatrix(quad_cfm, transquad)
3137 
3138  !if x2, y2 or z2, need to keep track of trace
3139  IF (i == 1 .OR. i == 4 .OR. i == 6) THEN
3140  osc_str(:) = osc_str(:) + real(transquad(2:ntot, 1))**2 + aimag(transquad(2:ntot, 1))**2
3141  trace(:) = trace(:) + transquad(2:ntot, 1)
3142 
3143  !if xy, xz, or yz, need to count twice (for yx, zx and zy)
3144  ELSE
3145  osc_str(:) = osc_str(:) + 2.0_dp*(real(transquad(2:ntot, 1))**2 + aimag(transquad(2:ntot, 1))**2)
3146  END IF
3147 
3148  END DO !i
3149 
3150  !remove a third of the trace
3151  osc_str(:) = osc_str(:) - 1._dp/3._dp*(real(trace(:))**2 + aimag(trace(:))**2)
3152 
3153  !multiply by the prefactor
3154  osc_str(:) = osc_str(:)*1._dp/20._dp*a_fine**2*soc_evals(:)**3
3155 
3156  !clean-up
3157  CALL cp_fm_struct_release(quad_struct)
3158  CALL cp_cfm_release(work1_cfm)
3159  CALL cp_cfm_release(work2_cfm)
3160  CALL cp_cfm_release(quad_cfm)
3161  CALL cp_fm_release(amew_quad)
3162  DEALLOCATE (transquad, trace)
3163 
3164  CALL timestop(handle)
3165 
3166  END SUBROUTINE compute_soc_quadrupole_fosc
3167 
3168 END MODULE xas_tdp_utils
3169 
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:117
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition: cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition: cp_cfm_diag.F:52
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition: cp_cfm_types.F:12
cp_cfm_types::cp_cfm_get_submatrix
subroutine, public cp_cfm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
Extract a sub-matrix from the full matrix: op(target_m)(1:n_rows,1:n_cols) = fm(start_row:start_row+n...
Definition: cp_cfm_types.F:333
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition: cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition: cp_cfm_types.F:159
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition: cp_cfm_types.F:817
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition: cp_cfm_types.F:607
cp_dbcsr_cholesky
Interface to (sca)lapack for the Cholesky based procedures.
Definition: cp_dbcsr_cholesky.F:17
cp_dbcsr_cholesky::cp_dbcsr_cholesky_decompose
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition: cp_dbcsr_cholesky.F:61
cp_dbcsr_cholesky::cp_dbcsr_cholesky_invert
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, upper_to_full)
used to replace the cholesky decomposition by the inverse
Definition: cp_dbcsr_cholesky.F:114
cp_dbcsr_diag
Interface to (sca)lapack for the Cholesky based procedures.
Definition: cp_dbcsr_diag.F:17
cp_dbcsr_diag::cp_dbcsr_power
subroutine, public cp_dbcsr_power(matrix, exponent, threshold, n_dependent, para_env, blacs_env, verbose, eigenvectors, eigenvalues)
...
Definition: cp_dbcsr_diag.F:233
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition: cp_dbcsr_operations.F:118
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition: cp_fm_basic_linalg.F:1816
cp_fm_basic_linalg::cp_fm_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition: cp_fm_basic_linalg.F:1644
cp_fm_basic_linalg::cp_fm_upper_to_full
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
Definition: cp_fm_basic_linalg.F:1693
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition: cp_fm_basic_linalg.F:1614
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
cp_fm_diag::cp_fm_geeig
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition: cp_fm_diag.F:1274
cp_fm_diag::choose_eigv_solver
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:208
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
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
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition: cp_fm_types.F:570
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:1016
cp_fm_types::cp_fm_to_fm_submat
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
Definition: cp_fm_types.F:1473
cp_fm_types::cp_fm_get_submatrix
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,...
Definition: cp_fm_types.F:901
cp_fm_types::cp_fm_set_element
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
Definition: cp_fm_types.F:700
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition: cp_log_handling.F:515
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::tddfpt_singlet
integer, parameter, public tddfpt_singlet
Definition: input_constants.F:575
input_constants::xas_dip_len
integer, parameter, public xas_dip_len
Definition: input_constants.F:682
input_constants::tddfpt_triplet
integer, parameter, public tddfpt_triplet
Definition: input_constants.F:575
input_constants::tddfpt_spin_flip
integer, parameter, public tddfpt_spin_flip
Definition: input_constants.F:575
input_constants::ot_precond_full_single
integer, parameter, public ot_precond_full_single
Definition: input_constants.F:513
input_constants::tddfpt_spin_cons
integer, parameter, public tddfpt_spin_cons
Definition: input_constants.F:575
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::get_diag
pure real(kind=dp) function, dimension(min(size(a, 1), size(a, 2))), public get_diag(a)
Return the diagonal elements of matrix a as a vector.
Definition: mathlib.F:493
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::a_fine
real(kind=dp), parameter, public a_fine
Definition: physcon.F:119
preconditioner_types
types of preconditioners
Definition: preconditioner_types.F:14
preconditioner_types::init_preconditioner
subroutine, public init_preconditioner(preconditioner_env, para_env, blacs_env)
...
Definition: preconditioner_types.F:86
preconditioner_types::destroy_preconditioner
subroutine, public destroy_preconditioner(preconditioner_env)
...
Definition: preconditioner_types.F:117
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition: preconditioner.F:15
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition: qs_mo_methods.F:14
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_ot_eigensolver
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
Definition: qs_ot_eigensolver.F:13
qs_ot_eigensolver::ot_eigensolver
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
Definition: qs_ot_eigensolver.F:73
xas_tdp_kernel
All the kernel specific subroutines for XAS TDP calculations.
Definition: xas_tdp_kernel.F:13
xas_tdp_kernel::kernel_coulomb_xc
subroutine, public kernel_coulomb_xc(coul_ker, xc_ker, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
Computes, if asked for it, the Coulomb and XC kernel matrices, in the usuall matrix format.
Definition: xas_tdp_kernel.F:81
xas_tdp_kernel::kernel_exchange
subroutine, public kernel_exchange(ex_ker, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
Computes the exact exchange kernel matrix using RI. Returns an array of 2 matrices,...
Definition: xas_tdp_kernel.F:595
xas_tdp_types
Define XAS TDP control type and associated create, release, etc subroutines, as well as XAS TDP envir...
Definition: xas_tdp_types.F:13
xas_tdp_utils
Utilities for X-ray absorption spectroscopy using TDDFPT.
Definition: xas_tdp_utils.F:13
xas_tdp_utils::include_rcs_soc
subroutine, public include_rcs_soc(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
Includes the SOC effects on the precomputed restricted closed-shell singlet and triplet excitations....
Definition: xas_tdp_utils.F:1939
xas_tdp_utils::setup_xas_tdp_prob
subroutine, public setup_xas_tdp_prob(donor_state, qs_env, xas_tdp_env, xas_tdp_control)
Builds the matrix that defines the XAS TDDFPT generalized eigenvalue problem to be solved for excitat...
Definition: xas_tdp_utils.F:168
xas_tdp_utils::solve_xas_tdp_prob
subroutine, public solve_xas_tdp_prob(donor_state, xas_tdp_control, xas_tdp_env, qs_env, ex_type)
Solves the XAS TDP generalized eigenvalue problem omega*C = matrix_tdp*C using standard full diagonal...
Definition: xas_tdp_utils.F:458
xas_tdp_utils::include_os_soc
subroutine, public include_os_soc(donor_state, xas_tdp_env, xas_tdp_control, qs_env)
Includes the SOC effects on the precomputed spin-conserving and spin-flip excitations from an open-sh...
Definition: xas_tdp_utils.F:1502
xas_tdp_utils::rcs_amew_soc_elements
subroutine, public rcs_amew_soc_elements(amew_soc, lr_soc, lr_overlap, domo_soc, pref_trace, pref_overall, pref_diags, symmetric)
Computes the rcs SOC matrix elements between excited states AMEWs based on the LR orbitals.
Definition: xas_tdp_utils.F:2883