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