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