(git:ed6f26b)
Loading...
Searching...
No Matches
qs_tddfpt2_soc.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7MODULE qs_tddfpt2_soc
8 USE basis_set_types, ONLY: get_gto_basis_set,&
9 gto_basis_set_type
10 USE cp_blacs_env, ONLY: cp_blacs_env_type
11 USE cp_cfm_diag, ONLY: cp_cfm_heevd
12 USE cp_cfm_types, ONLY: cp_cfm_create,&
13 cp_cfm_get_info,&
14 cp_cfm_get_submatrix,&
15 cp_cfm_release,&
16 cp_cfm_type,&
17 cp_fm_to_cfm
18 USE cp_control_types, ONLY: dft_control_type,&
19 qs_control_type,&
20 tddfpt2_control_type
21 USE cp_dbcsr_api, ONLY: &
22 dbcsr_copy, dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, &
23 dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_get_info, dbcsr_multiply, &
24 dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_type, dbcsr_type_no_symmetry
25 USE cp_dbcsr_contrib, ONLY: dbcsr_print,&
26 dbcsr_reserve_all_blocks
27 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 copy_fm_to_dbcsr,&
29 cp_dbcsr_sm_fm_multiply
30 USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
31 cp_fm_transpose,&
32 cp_fm_uplo_to_full
33 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
34 cp_fm_struct_release,&
35 cp_fm_struct_type
36 USE cp_fm_types, ONLY: &
37 cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, &
38 cp_fm_set_all, cp_fm_set_element, cp_fm_to_fm_submat, cp_fm_type
39 USE cp_log_handling, ONLY: cp_get_default_logger,&
40 cp_logger_get_default_unit_nr,&
41 cp_logger_type
42 USE cp_output_handling, ONLY: cp_print_key_unit_nr
43 USE input_constants, ONLY: tddfpt_dipole_length
44 USE input_section_types, ONLY: section_vals_get_subs_vals,&
45 section_vals_type,&
46 section_vals_val_get
47 USE kinds, ONLY: default_string_length,&
48 dp
49 USE lebedev, ONLY: deallocate_lebedev_grids,&
50 get_number_of_lebedev_grid,&
51 init_lebedev_grids,&
52 lebedev_grid
53 USE mathlib, ONLY: get_diag
54 USE memory_utilities, ONLY: reallocate
55 USE message_passing, ONLY: mp_para_env_type
56 USE orbital_pointers, ONLY: indso,&
57 nsoset
58 USE parallel_gemm_api, ONLY: parallel_gemm
59 USE particle_types, ONLY: particle_type
60 USE physcon, ONLY: evolt,&
61 wavenumbers
62 USE qs_environment_types, ONLY: get_qs_env,&
63 qs_environment_type
64 USE qs_grid_atom, ONLY: allocate_grid_atom,&
65 create_grid_atom,&
66 grid_atom_type
67 USE qs_harmonics_atom, ONLY: allocate_harmonics_atom,&
68 create_harmonics_atom,&
69 get_maxl_cg,&
70 harmonics_atom_type
71 USE qs_kind_types, ONLY: get_qs_kind,&
72 get_qs_kind_set,&
73 qs_kind_type
74 USE qs_mo_types, ONLY: get_mo_set,&
75 mo_set_type
76 USE qs_tddfpt2_soc_types, ONLY: soc_atom_create,&
77 soc_atom_env_type,&
78 soc_atom_release,&
79 soc_env_create,&
80 soc_env_release,&
81 soc_env_type
82 USE qs_tddfpt2_soc_utils, ONLY: dip_vel_op,&
83 resort_evects,&
84 soc_contract_evect,&
85 soc_dipole_operator
86 USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos
87 USE soc_pseudopotential_methods, ONLY: v_soc_xyz_from_pseudopotential
88 USE spherical_harmonics, ONLY: clebsch_gordon,&
89 clebsch_gordon_deallocate,&
90 clebsch_gordon_init
91 USE xas_tdp_atom, ONLY: compute_sphi_so,&
92 integrate_soc_atoms,&
93 truncate_radial_grid
94 USE xas_tdp_utils, ONLY: rcs_amew_soc_elements
95#include "./base/base_uses.f90"
96
97 IMPLICIT NONE
98
99 PRIVATE
100
101 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_soc'
102
103 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
104
105 INTEGER, PARAMETER, PRIVATE :: nderivs = 3
106 INTEGER, PARAMETER, PRIVATE :: maxspins = 2
107
108 !A helper type for SOC
109 TYPE dbcsr_soc_package_type
110 TYPE(dbcsr_type), POINTER :: dbcsr_sg => null()
111 TYPE(dbcsr_type), POINTER :: dbcsr_tp => null()
112 TYPE(dbcsr_type), POINTER :: dbcsr_sc => null()
113 TYPE(dbcsr_type), POINTER :: dbcsr_sf => null()
114 TYPE(dbcsr_type), POINTER :: dbcsr_prod => null()
115 TYPE(dbcsr_type), POINTER :: dbcsr_ovlp => null()
116 TYPE(dbcsr_type), POINTER :: dbcsr_tmp => null()
117 TYPE(dbcsr_type), POINTER :: dbcsr_work => null()
118 END TYPE dbcsr_soc_package_type
119
120 PUBLIC :: tddfpt_soc
121
122! **************************************************************************************************
123
124CONTAINS
125
126! **************************************************************************************************
127!> \brief Perform TDDFPT-SOC calculation.
128!> \param qs_env Quickstep environment
129!> \param evals_a eigenvalues for the singlet states
130!> \param evals_b eigenvalues for the triplet states
131!> \param evects_a eigenvectors for the singlet states
132!> \param evects_b eigenvectors for the triplet states
133!> \param gs_mos ground state orbitlas from the TDDFPT calculation
134!> \par History
135!> * 02.2023 created [Jan-Robert Vogt]
136!> \note Based on tddfpt2_methods and xas_tdp_utils.
137!> \note Only rcs for now, but written with the addition of os in mind
138! **************************************************************************************************
139
140 SUBROUTINE tddfpt_soc(qs_env, evals_a, evals_b, evects_a, evects_b, gs_mos)
141
142 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
143 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: evals_a, evals_b
144 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects_a, evects_b
145 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
146 POINTER :: gs_mos
147
148 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_soc'
149
150 INTEGER :: handle, istate, log_unit
151 LOGICAL :: do_os
152 TYPE(cp_logger_type), POINTER :: logger
153 TYPE(dft_control_type), POINTER :: dft_control
154
155 CALL timeset(routinen, handle)
156
157 logger => cp_get_default_logger()
158 IF (logger%para_env%is_source()) THEN
159 log_unit = cp_logger_get_default_unit_nr(logger, local=.true.)
160 ELSE
161 log_unit = -1
162 END IF
163
164 IF (log_unit > 0) THEN
165 WRITE (log_unit, "(1X,A)") "", &
166 "-------------------------------------------------------------------------------", &
167 "- START SOC CALCULATIONS -", &
168 "-------------------------------------------------------------------------------"
169 END IF
170
171 NULLIFY (dft_control)
172 CALL get_qs_env(qs_env, dft_control=dft_control)
173 do_os = dft_control%uks .OR. dft_control%roks
174 IF (do_os) cpabort("SOC only implemented for closed shell.")
175
176 IF (log_unit > 0) THEN
177 WRITE (log_unit, '(A)') "Starting from TDDFPT Excited States:"
178 WRITE (log_unit, '(A)') " STATE SINGLET/eV TRIPLET/eV"
179 DO istate = 1, SIZE(evals_a)
180 WRITE (log_unit, '(6X,I3,11X,F10.5,6X,F10.5)') istate, evals_a(istate)*evolt, evals_b(istate)*evolt
181 END DO
182 END IF
183
184 IF (log_unit > 0) WRITE (log_unit, '(A)') "Starting restricted closed shell:"
185 CALL tddfpt_soc_rcs(qs_env, evals_a, evals_b, evects_a, evects_b, log_unit, gs_mos)
186
187 IF (log_unit > 0) THEN
188 WRITE (log_unit, '(A,/,A)') "SOC Calculation terminated", &
189 "Returning to TDDFPT for Force calculation and deallocations"
190 END IF
191
192 CALL timestop(handle)
193
194 END SUBROUTINE
195
196! **************************************************************************************************
197!> \brief Will perform the soc-calculation for restricted-closed-shell systems
198!> \param qs_env Quickstep Enviroment
199!> \param evals_sing eigenvalues singlet states
200!> \param evals_trip eigenvalues triplet states
201!> \param evects_sing eigenvector singlet states
202!> \param evects_trip eigenvectors triplet states
203!> \param log_unit default log_unit for convinients
204!> \param gs_mos ground state MOs from TDDFPT
205!> \par History
206!> * created 02.2023 [Jan-Robert Vogt]
207!> \note Mostly copied and modified from xas_tdp_utils include_rcs_soc
208! **************************************************************************************************
209 SUBROUTINE tddfpt_soc_rcs(qs_env, evals_sing, evals_trip, evects_sing, evects_trip, log_unit, gs_mos)
210
211 TYPE(qs_environment_type), POINTER :: qs_env
212 REAL(kind=dp), DIMENSION(:), INTENT(IN), TARGET :: evals_sing, evals_trip
213 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects_sing, evects_trip
214 INTEGER :: log_unit
215 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
216 POINTER :: gs_mos
217
218 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_soc_rcs'
219
220 CHARACTER(len=3) :: mult
221 INTEGER :: dipole_form, group, handle, iex, ii, &
222 isg, istate, itp, jj, nactive, nao, &
223 nex, npcols, nprows, nsg, ntot, ntp
224 INTEGER, ALLOCATABLE, DIMENSION(:) :: evects_sort
225 INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_dist, row_blk_size, &
226 row_dist, row_dist_new
227 INTEGER, DIMENSION(:, :), POINTER :: pgrid
228 LOGICAL :: print_ev, print_some, print_splitting, &
229 print_wn
230 REAL(dp) :: eps_filter, soc_gst, sqrt2, unit_scale
231 REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag, tmp_evals
232 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: img_tmp, real_tmp
233 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gstp_block, mo_soc_x, mo_soc_y, mo_soc_z
234 TYPE(cp_blacs_env_type), POINTER :: blacs_env
235 TYPE(cp_cfm_type) :: evecs_cfm, hami_cfm
236 TYPE(cp_fm_struct_type), POINTER :: full_struct, gstp_struct, prod_struct, &
237 vec_struct, work_struct
238 TYPE(cp_fm_type) :: gstp_fm, img_fm, prod_fm, real_fm, &
239 tmp_fm, vec_soc_x, vec_soc_y, &
240 vec_soc_z, work_fm
241 TYPE(cp_fm_type), POINTER :: gs_coeffs, tp_coeffs
242 TYPE(cp_logger_type), POINTER :: logger
243 TYPE(dbcsr_distribution_type), POINTER :: coeffs_dist, dbcsr_dist, prod_dist
244 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
245 TYPE(dbcsr_soc_package_type) :: dbcsr_soc_package
246 TYPE(dbcsr_type), POINTER :: dbcsr_dummy, dbcsr_ovlp, dbcsr_prod, &
247 dbcsr_sg, dbcsr_tmp, dbcsr_tp, &
248 dbcsr_work, orb_soc_x, orb_soc_y, &
249 orb_soc_z
250 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
251 TYPE(mp_para_env_type), POINTER :: para_env
252 TYPE(section_vals_type), POINTER :: soc_print_section, soc_section, &
253 tddfpt_print_section
254 TYPE(soc_atom_env_type), POINTER :: soc_atom_env
255 TYPE(soc_env_type), TARGET :: soc_env
256
257 CALL timeset(routinen, handle)
258
259 NULLIFY (logger, pgrid, row_dist, row_dist_new)
260 NULLIFY (col_dist, col_blk_size, row_blk_size, matrix_s)
261 NULLIFY (gstp_struct, tddfpt_print_section, soc_print_section, soc_section)
262
263 logger => cp_get_default_logger()
264 tddfpt_print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
265 soc_print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT%SOC_PRINT")
266 soc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%SOC")
267
268 CALL section_vals_val_get(soc_section, "EPS_FILTER", r_val=eps_filter)
269
270 nsg = SIZE(evals_sing) ! Number of excited singlet states
271 ntp = SIZE(evals_trip) ! Number of excited triplet states
272 nex = nsg ! Number of excited states of each multiplicity
273 ntot = 1 + nsg + 3*ntp ! Number of (GS + S + T^-1 + T^0 + T^1)
274
275 ! Initzialize Working environment
276 CALL inititialize_soc(qs_env, soc_atom_env, soc_env, &
277 evects_sing, evects_trip, dipole_form, gs_mos)
278
279 NULLIFY (mos)
280 CALL get_qs_env(qs_env, mos=mos, para_env=para_env, blacs_env=blacs_env)
281 CALL get_mo_set(mos(1), nao=nao, homo=nactive)
282
283 ! this will create the H^SOC in an atomic basis
284 CALL integrate_soc_atoms(soc_env%orb_soc, qs_env=qs_env, soc_atom_env=soc_atom_env)
285 CALL soc_atom_release(soc_atom_env)
286
287 !! Point at H^SOC and MOs for better readablity
288 NULLIFY (tp_coeffs, orb_soc_x, orb_soc_y, orb_soc_z)
289 tp_coeffs => soc_env%b_coeff
290 soc_env%evals_a => evals_sing
291 soc_env%evals_b => evals_trip
292 orb_soc_x => soc_env%orb_soc(1)%matrix
293 orb_soc_y => soc_env%orb_soc(2)%matrix
294 orb_soc_z => soc_env%orb_soc(3)%matrix
295
296 !! Create a matrix-structure, which links all states in this calculation
297 NULLIFY (full_struct)
298 CALL cp_fm_struct_create(full_struct, context=blacs_env, para_env=para_env, nrow_global=ntot, ncol_global=ntot)
299 CALL cp_fm_create(real_fm, full_struct)
300 CALL cp_fm_create(img_fm, full_struct)
301 CALL cp_fm_set_all(real_fm, 0.0_dp)
302 CALL cp_fm_set_all(img_fm, 0.0_dp)
303
304 ! Put the excitation energies on the diagonal of the real matrix
305 DO isg = 1, nsg
306 CALL cp_fm_set_element(real_fm, 1 + isg, 1 + isg, evals_sing(isg))
307 END DO
308 DO itp = 1, ntp
309 ! first T^-1, then T^0, then T^+1
310 CALL cp_fm_set_element(real_fm, 1 + itp + nsg, 1 + itp + nsg, evals_trip(itp))
311 CALL cp_fm_set_element(real_fm, 1 + itp + ntp + nsg, 1 + itp + ntp + nsg, evals_trip(itp))
312 CALL cp_fm_set_element(real_fm, 1 + itp + 2*ntp + nsg, 1 + itp + 2*ntp + nsg, evals_trip(itp))
313 END DO
314
315 !!Create the dbcsr structures for this calculations
316 CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist)
317 CALL dbcsr_distribution_get(dbcsr_dist, group=group, row_dist=row_dist, pgrid=pgrid, &
318 npcols=npcols, nprows=nprows)
319 ALLOCATE (col_dist(nex), row_dist_new(nex)) ! Split for each excitation
320 DO iex = 1, nex
321 col_dist(iex) = modulo(npcols - iex, npcols)
322 row_dist_new(iex) = modulo(nprows - iex, nprows)
323 END DO
324 ALLOCATE (coeffs_dist, prod_dist)
325 CALL dbcsr_distribution_new(coeffs_dist, group=group, pgrid=pgrid, row_dist=row_dist, &
326 col_dist=col_dist)
327 CALL dbcsr_distribution_new(prod_dist, group=group, pgrid=pgrid, row_dist=row_dist_new, &
328 col_dist=col_dist)
329
330 !! Create the matrices
331 ALLOCATE (col_blk_size(nex))
332 col_blk_size = nactive
333
334 CALL get_qs_env(qs_env, matrix_s=matrix_s)
335 CALL dbcsr_get_info(matrix_s(1)%matrix, row_blk_size=row_blk_size)
336
337 !! The Eigenvectors for Sg und Tp will be dived into their diffrent components again
338 ALLOCATE (dbcsr_sg, dbcsr_tp, dbcsr_work, dbcsr_ovlp, dbcsr_tmp, dbcsr_prod)
339 CALL dbcsr_create(matrix=dbcsr_sg, name="SINGLETS", matrix_type=dbcsr_type_no_symmetry, &
340 dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
341 CALL dbcsr_create(matrix=dbcsr_tp, name="TRIPLETS", matrix_type=dbcsr_type_no_symmetry, &
342 dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
343 CALL dbcsr_create(matrix=dbcsr_work, name="WORK", matrix_type=dbcsr_type_no_symmetry, &
344 dist=coeffs_dist, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
345 CALL dbcsr_create(matrix=dbcsr_prod, name="PROD", matrix_type=dbcsr_type_no_symmetry, &
346 dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
347 CALL dbcsr_create(matrix=dbcsr_ovlp, name="OVLP", matrix_type=dbcsr_type_no_symmetry, &
348 dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
349
350 col_blk_size = 1
351 CALL dbcsr_create(matrix=dbcsr_tmp, name="TMP", matrix_type=dbcsr_type_no_symmetry, &
352 dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
353 CALL dbcsr_reserve_all_blocks(dbcsr_tmp)
354
355 IF (debug_this_module) THEN
356 ALLOCATE (dbcsr_dummy)
357 CALL dbcsr_create(matrix=dbcsr_dummy, name="DUMMY", matrix_type=dbcsr_type_no_symmetry, &
358 dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
359 CALL dbcsr_reserve_all_blocks(dbcsr_dummy)
360 END IF
361
362 !! This work dbcsr matrix will be packed together for easy transfer to other subroutines
363 dbcsr_soc_package%dbcsr_sg => dbcsr_sg
364 dbcsr_soc_package%dbcsr_tp => dbcsr_tp
365 dbcsr_soc_package%dbcsr_work => dbcsr_work
366 dbcsr_soc_package%dbcsr_ovlp => dbcsr_ovlp
367 dbcsr_soc_package%dbcsr_prod => dbcsr_prod
368 dbcsr_soc_package%dbcsr_tmp => dbcsr_tmp
369
370 !Filling the coeffs matrices by copying from the stored fms
371 CALL copy_fm_to_dbcsr(soc_env%a_coeff, dbcsr_sg)
372 CALL copy_fm_to_dbcsr(soc_env%b_coeff, dbcsr_tp)
373
374 !Create the work and helper fms
375 !CALL get_mo_set(mos(1),mo_coeff=gs_coeffs)
376 gs_coeffs => gs_mos(1)%mos_occ
377 CALL cp_fm_get_info(gs_coeffs, matrix_struct=vec_struct)
378 CALL cp_fm_create(vec_soc_x, vec_struct)
379 CALL cp_fm_create(vec_soc_y, vec_struct)
380 CALL cp_fm_create(vec_soc_z, vec_struct)
381
382 CALL cp_fm_struct_create(prod_struct, context=blacs_env, para_env=para_env, &
383 nrow_global=nactive, ncol_global=nactive)
384 CALL cp_fm_create(prod_fm, prod_struct)
385
386 CALL cp_fm_struct_create(work_struct, context=blacs_env, para_env=para_env, &
387 nrow_global=nex, ncol_global=nex)
388 CALL cp_fm_create(work_fm, work_struct)
389 CALL cp_fm_create(tmp_fm, work_struct)
390
391 IF (log_unit > 0 .AND. debug_this_module) THEN
392 WRITE (log_unit, '(A)') "Starting Precomputations"
393 END IF
394
395 !! Begin with the precomputation
396 !! Prefactor due to rcs.
397 !! The excitation is presented as a linear combination of
398 !! alpha and beta, where dalpha=-dbeta for triplet exc.
399 !! and dalpha=dbeta for the singlet case
400 sqrt2 = sqrt(2.0_dp)
401
402 !! Precompute the <phi_i^0|H^SOC|phi_j^0> matrix elements
403 ALLOCATE (diag(nactive))
404 ALLOCATE (mo_soc_x(nactive, nactive), mo_soc_y(nactive, nactive), mo_soc_z(nactive, nactive))
405
406 !! This will be the GS|H|GS contribution needed for all couplings
407 CALL cp_dbcsr_sm_fm_multiply(orb_soc_x, gs_coeffs, vec_soc_x, ncol=nactive)
408 CALL parallel_gemm('T', 'N', nactive, &
409 nactive, nao, 1.0_dp, &
410 gs_coeffs, vec_soc_x, &
411 0.0_dp, prod_fm)
412 CALL cp_fm_get_submatrix(prod_fm, mo_soc_x)
413
414 CALL cp_dbcsr_sm_fm_multiply(orb_soc_y, gs_coeffs, vec_soc_y, ncol=nactive)
415 CALL parallel_gemm('T', 'N', nactive, &
416 nactive, nao, 1.0_dp, &
417 gs_coeffs, vec_soc_y, &
418 0.0_dp, prod_fm)
419 CALL cp_fm_get_submatrix(prod_fm, mo_soc_y)
420
421 CALL cp_dbcsr_sm_fm_multiply(orb_soc_z, gs_coeffs, vec_soc_z, ncol=nactive)
422 CALL parallel_gemm('T', 'N', nactive, &
423 nactive, nao, 1.0_dp, &
424 gs_coeffs, vec_soc_z, &
425 0.0_dp, prod_fm)
426 CALL cp_fm_get_submatrix(prod_fm, mo_soc_z)
427
428 ! Only have SOC between singlet-triplet triplet-triplet and ground_state-triplet, the resulting
429 ! matrix is Hermitian i.e. the real part is symmetric and the imaginary part is anti-symmetric.
430 ! Can only fill upper half
431 IF (log_unit > 0 .AND. debug_this_module) THEN
432 WRITE (log_unit, '(A)') "Starting Ground-State contributions..."
433 END IF
434 !Start with the ground state/triplet SOC, SOC*gs_coeffs already computed above
435 !note: we are computing <0|H|T>, but have SOC*gs_coeffs instead of gs_coeffs*SOC in store.
436 !Since the SOC Hamiltonian is anti-symmetric, a - signs pops up in the gemms below
437 CALL cp_fm_struct_create(gstp_struct, context=blacs_env, para_env=para_env, &
438 nrow_global=ntp*nactive, ncol_global=nactive)
439 CALL cp_fm_create(gstp_fm, gstp_struct)
440 ALLOCATE (gstp_block(nactive, nactive))
441
442 !gs-triplet with Ms=+-1, imaginary part
443 ! <T+-1|H_x|GS>
444 ! -1 to change to <GS|H|T+-1>
445 CALL parallel_gemm('T', 'N', nactive*ntp, &
446 nactive, nao, -1.0_dp, &
447 tp_coeffs, vec_soc_x, &
448 0.0_dp, gstp_fm)
449
450 !! Seperate them into the different states again (nactive x nactive)
451 DO itp = 1, ntp
452 CALL cp_fm_get_submatrix(fm=gstp_fm, target_m=gstp_block, &
453 start_row=(itp - 1)*nactive + 1, start_col=1, &
454 n_rows=nactive, n_cols=nactive)
455 diag(:) = get_diag(gstp_block)
456 soc_gst = sum(diag)
457 !! <0|H_x|T^-1>
458 CALL cp_fm_set_element(img_fm, 1, 1 + nsg + itp, -1_dp*soc_gst)
459 !! <0|H_x|T^+1>
460 CALL cp_fm_set_element(img_fm, 1, 1 + nsg + 2*ntp + itp, soc_gst)
461 END DO
462
463 IF (debug_this_module .AND. (log_unit > 0)) THEN
464 WRITE (log_unit, "(A,2F18.6)") "<0|H_x|T^-1> in Hartree / cm^-1", soc_gst, soc_gst*wavenumbers
465 WRITE (log_unit, "(A,2F18.6)") "<0|H_x|T^+1> in Hartree / cm^-1", soc_gst, soc_gst*wavenumbers
466 END IF
467
468 !gs-triplet with Ms=+-1, real part
469 ! <T+-1|H_y|GS>
470 CALL parallel_gemm('T', 'N', nactive*ntp, &
471 nactive, nao, -1.0_dp, &
472 tp_coeffs, vec_soc_y, &
473 0.0_dp, gstp_fm)
474
475 DO itp = 1, ntp
476 CALL cp_fm_get_submatrix(fm=gstp_fm, target_m=gstp_block, start_row=(itp - 1)*nactive + 1, &
477 start_col=1, n_rows=nactive, n_cols=nactive)
478 diag(:) = get_diag(gstp_block)
479 soc_gst = sum(diag)
480 ! <0|H_y|T^-1>
481 CALL cp_fm_set_element(real_fm, 1, 1 + nsg + itp, -1.0_dp*soc_gst)
482 ! <0|H_y|T^+1>
483 CALL cp_fm_set_element(real_fm, 1, 1 + nsg + 2*ntp + itp, -1.0_dp*soc_gst)
484 END DO
485
486 IF (debug_this_module .AND. (log_unit > 0)) THEN
487 WRITE (log_unit, "(A,2F18.6)") "<0|H_y|T^-1> in Hartree / cm^-1", soc_gst, soc_gst*wavenumbers
488 WRITE (log_unit, "(A,2F18.6)") "<0|H_y|T^+1> in Hartree / cm^-1", soc_gst, soc_gst*wavenumbers
489 END IF
490
491 !gs-triplet with Ms=0, purely imaginary
492 !< T0|H_z|GS>
493 CALL parallel_gemm('T', 'N', nactive*ntp, &
494 nactive, nao, -1.0_dp, &
495 tp_coeffs, vec_soc_z, &
496 0.0_dp, gstp_fm)
497
498 DO itp = 1, ntp
499 CALL cp_fm_get_submatrix(fm=gstp_fm, &
500 target_m=gstp_block, &
501 start_row=(itp - 1)*nactive + 1, &
502 start_col=1, &
503 n_rows=nactive, &
504 n_cols=nactive)
505 diag(:) = get_diag(gstp_block)
506 soc_gst = sqrt2*sum(diag)
507 CALL cp_fm_set_element(img_fm, 1, 1 + nsg + ntp + itp, soc_gst)
508 END DO
509
510 IF (debug_this_module .AND. (log_unit > 0)) THEN
511 WRITE (log_unit, "(A,2F18.6)") "<0|H_z|T> in Hartree / cm^-1", soc_gst, soc_gst*wavenumbers
512 END IF
513
514 !! After all::
515 !! T-1 :: -<0|H_x|T^-1> - <0|H_y|T^-1> at 1+nsg+itp
516 !! T+1 :: <0|H_x|T^-1> - <0|H_y|T^-1> at 1+nsp+2tnp+itp
517 !! T0 :: < T0|H_z|GS> at 1+nsg+ntp+itp
518
519 !gs clean-up
520 CALL cp_fm_release(prod_fm)
521 CALL cp_fm_release(vec_soc_x)
522 CALL cp_fm_release(vec_soc_y)
523 CALL cp_fm_release(vec_soc_z)
524 CALL cp_fm_release(gstp_fm)
525 CALL cp_fm_struct_release(gstp_struct)
526 CALL cp_fm_struct_release(prod_struct)
527 DEALLOCATE (gstp_block)
528
529 IF (log_unit > 0 .AND. debug_this_module) THEN
530 WRITE (log_unit, '(A)') "Starting Singlet-Triplet contributions..."
531 END IF
532
533 !Now do the singlet-triplet SOC
534 !start by computing the singlet-triplet overlap
535 ! <S|T>
536 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
537 matrix_s(1)%matrix, &
538 dbcsr_tp, 0.0_dp, &
539 dbcsr_work, &
540 filter_eps=eps_filter)
541 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
542 dbcsr_sg, &
543 dbcsr_work, &
544 0.0_dp, dbcsr_ovlp, &
545 filter_eps=eps_filter)
546
547 !singlet-triplet with Ms=+-1, imaginary part
548 ! First precalculate <S|H_x|T+-1>
549 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
550 orb_soc_x, &
551 dbcsr_tp, 0.0_dp, &
552 dbcsr_work, &
553 filter_eps=eps_filter)
554 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
555 dbcsr_sg, &
556 dbcsr_work, 0.0_dp, &
557 dbcsr_prod, &
558 filter_eps=eps_filter)
559
560 !! This will lead to:
561 !! -1/sqrt(2)(<S|H_x|T> - <S|T> <GS|H_x|GS>)
562 CALL rcs_amew_soc_elements(dbcsr_tmp, &
563 dbcsr_prod, &
564 dbcsr_ovlp, &
565 mo_soc_x, &
566 pref_trace=-1.0_dp, &
567 pref_overall=-0.5_dp*sqrt2)
568
569 !<S|H_x|T^-1>
570 !! Convert to fm for transfer to img_fm
571 CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
572 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
573 mtarget=img_fm, &
574 nrow=nex, &
575 ncol=nex, &
576 s_firstrow=1, &
577 s_firstcol=1, &
578 t_firstrow=2, &
579 t_firstcol=1 + nsg + 1)
580
581 IF (debug_this_module) THEN
582 WRITE (log_unit, "(/,A)") "<S|H_x|T^-1> in Hartree / cm^-1"
583 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
584 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
585 CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
586 END IF
587
588 !<S|H_x|T^+1> takes a minus sign
589 CALL cp_fm_scale(-1.0_dp, tmp_fm)
590 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
591 mtarget=img_fm, &
592 nrow=nex, &
593 ncol=nex, &
594 s_firstrow=1, &
595 s_firstcol=1, &
596 t_firstrow=2, &
597 t_firstcol=1 + nsg + 2*ntp + 1)
598
599 IF (debug_this_module) THEN
600 WRITE (log_unit, "(/,A)") "<S|H_x|T^+1> in Hartree / cm^-1"
601 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
602 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
603 CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
604 END IF
605
606 !!singlet-triplet with Ms=+-1, real part
607 !! Precompute <S|H_y|T>
608 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
609 orb_soc_y, dbcsr_tp, &
610 0.0_dp, dbcsr_work, &
611 filter_eps=eps_filter)
612 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
613 dbcsr_sg, dbcsr_work, &
614 0.0_dp, dbcsr_prod, &
615 filter_eps=eps_filter)
616
617 !! This will lead to -1/sqrt(2)(<S|H_y|T> - <S|T> <GS|H_y|GS>)
618 CALL rcs_amew_soc_elements(dbcsr_tmp, &
619 dbcsr_prod, &
620 dbcsr_ovlp, &
621 mo_soc_y, &
622 pref_trace=-1.0_dp, &
623 pref_overall=-0.5_dp*sqrt2)
624
625 !<S|H_y|T^-1>
626 CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
627 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
628 mtarget=real_fm, &
629 nrow=nex, &
630 ncol=nex, &
631 s_firstrow=1, &
632 s_firstcol=1, &
633 t_firstrow=2, &
634 t_firstcol=1 + nsg + 1)
635
636 IF (debug_this_module) THEN
637 WRITE (log_unit, "(/,A)") "<S|H_y|T^-1> in Hartree / cm^-1"
638 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
639 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
640 CALL dbcsr_print(dbcsr_tmp, unit_nr=log_unit)
641 END IF
642
643 !<S|H_y|T^+1>
644 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
645 mtarget=real_fm, &
646 nrow=nex, &
647 ncol=nex, &
648 s_firstrow=1, &
649 s_firstcol=1, &
650 t_firstrow=2, &
651 t_firstcol=1 + nsg + 2*ntp + 1)
652
653 IF (debug_this_module) THEN
654 WRITE (log_unit, "(/,A)") "<S|H_y|T^+1> in Hartree / cm^-1"
655 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
656 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
657 CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
658 END IF
659
660 !singlet-triplet with Ms=0, purely imaginary
661 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
662 orb_soc_z, dbcsr_tp, &
663 0.0_dp, dbcsr_work, &
664 filter_eps=eps_filter)
665 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
666 dbcsr_sg, dbcsr_work, &
667 0.0_dp, dbcsr_prod, &
668 filter_eps=eps_filter)
669
670 CALL rcs_amew_soc_elements(dbcsr_tmp, &
671 dbcsr_prod, &
672 dbcsr_ovlp, &
673 mo_soc_z, &
674 pref_trace=-1.0_dp, &
675 pref_overall=1.0_dp)
676
677 !<S|H_z|T^0>
678 CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
679 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
680 mtarget=img_fm, &
681 nrow=nex, &
682 ncol=nex, &
683 s_firstrow=1, &
684 s_firstcol=1, &
685 t_firstrow=2, &
686 t_firstcol=1 + nsg + ntp + 1)
687
688 IF (debug_this_module) THEN
689 WRITE (log_unit, "(/,A)") "<S|H_z|T^0> in Hartree / cm^-1"
690 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
691 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
692 CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
693 END IF
694
695 IF (log_unit > 0 .AND. debug_this_module) THEN
696 WRITE (log_unit, '(A)') "Starting Triplet-Triplet contributions..."
697 END IF
698
699 !Now the triplet-triplet SOC
700 !start by computing the overlap
701 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
702 matrix_s(1)%matrix, &
703 dbcsr_tp, 0.0_dp, &
704 dbcsr_work, &
705 filter_eps=eps_filter)
706 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
707 dbcsr_tp, dbcsr_work, &
708 0.0_dp, dbcsr_ovlp, &
709 filter_eps=eps_filter)
710
711 !Ms=0 to Ms=+-1 SOC, imaginary part
712 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
713 orb_soc_x, dbcsr_tp, &
714 0.0_dp, dbcsr_work, &
715 filter_eps=eps_filter)
716 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
717 dbcsr_tp, dbcsr_work, &
718 0.0_dp, dbcsr_prod, &
719 filter_eps=eps_filter)
720
721 CALL rcs_amew_soc_elements(dbcsr_tmp, &
722 dbcsr_prod, &
723 dbcsr_ovlp, &
724 mo_soc_x, &
725 pref_trace=1.0_dp, &
726 pref_overall=-0.5_dp*sqrt2)
727
728 !<T^0|H_x|T^+1>
729 CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
730 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
731 mtarget=img_fm, &
732 nrow=nex, &
733 ncol=nex, &
734 s_firstrow=1, &
735 s_firstcol=1, &
736 t_firstrow=1 + nsg + ntp + 1, &
737 t_firstcol=1 + nsg + 2*ntp + 1)
738
739 IF (debug_this_module) THEN
740 WRITE (log_unit, "(/,A)") "<T^0|H_x|T^+1> in Hartree / cm^-1"
741 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
742 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
743 CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
744 END IF
745
746 !<T^-1|H_x|T^0>, takes a minus sign and a transpose (because computed <T^0|H_x|T^-1>)
747 CALL cp_fm_transpose(tmp_fm, work_fm)
748 CALL cp_fm_scale(-1.0_dp, work_fm)
749 CALL cp_fm_to_fm_submat(msource=work_fm, &
750 mtarget=img_fm, &
751 nrow=nex, &
752 ncol=nex, &
753 s_firstrow=1, &
754 s_firstcol=1, &
755 t_firstrow=1 + nsg + 1, &
756 t_firstcol=1 + nsg + ntp + 1)
757
758 IF (debug_this_module) THEN
759 WRITE (log_unit, "(/,A)") "<T^-1|H_x|T^0> in Hartree / cm^-1"
760 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
761 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
762 CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
763 END IF
764
765 !Ms=0 to Ms=+-1 SOC, real part
766 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
767 orb_soc_y, dbcsr_tp, &
768 0.0_dp, dbcsr_work, &
769 filter_eps=eps_filter)
770 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
771 dbcsr_tp, dbcsr_work, &
772 0.0_dp, dbcsr_prod, &
773 filter_eps=eps_filter)
774
775 CALL rcs_amew_soc_elements(dbcsr_tmp, &
776 dbcsr_prod, &
777 dbcsr_ovlp, &
778 mo_soc_y, &
779 pref_trace=1.0_dp, &
780 pref_overall=0.5_dp*sqrt2)
781
782 !<T^0|H_y|T^+1>
783 CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
784 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
785 mtarget=real_fm, &
786 nrow=nex, &
787 ncol=nex, &
788 s_firstrow=1, &
789 s_firstcol=1, t_firstrow=1 + nsg + ntp + 1, &
790 t_firstcol=1 + nsg + 2*ntp + 1)
791
792 IF (debug_this_module) THEN
793 WRITE (log_unit, "(/,A)") "<T^0|H_y|T^+1> in Hartree / cm^-1"
794 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
795 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
796 CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
797 END IF
798
799 !<T^-1|H_y|T^0>, takes a minus sign and a transpose
800 CALL cp_fm_transpose(tmp_fm, work_fm)
801 CALL cp_fm_scale(-1.0_dp, work_fm)
802 CALL cp_fm_to_fm_submat(msource=work_fm, &
803 mtarget=real_fm, &
804 nrow=nex, &
805 ncol=nex, &
806 s_firstrow=1, &
807 s_firstcol=1, t_firstrow=1 + nsg + 1, &
808 t_firstcol=1 + nsg + ntp + 1)
809
810 IF (debug_this_module) THEN
811 WRITE (log_unit, "(/,A)") "<T^0|H_y|T^+1> in Hartree / cm^-1"
812 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
813 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
814 CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
815 END IF
816
817 !Ms=1 to Ms=1 and Ms=-1 to Ms=-1 SOC, purely imaginary
818 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
819 orb_soc_z, dbcsr_tp, &
820 0.0_dp, dbcsr_work, &
821 filter_eps=eps_filter)
822 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
823 dbcsr_tp, dbcsr_work, &
824 0.0_dp, dbcsr_prod, &
825 filter_eps=eps_filter)
826
827 CALL rcs_amew_soc_elements(dbcsr_tmp, &
828 dbcsr_prod, &
829 dbcsr_ovlp, &
830 mo_soc_z, &
831 pref_trace=1.0_dp, &
832 pref_overall=1.0_dp)
833
834 !<T^+1|H_z|T^+1>
835 CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
836 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
837 mtarget=img_fm, &
838 nrow=nex, &
839 ncol=nex, &
840 s_firstrow=1, &
841 s_firstcol=1, &
842 t_firstrow=1 + nsg + 2*ntp + 1, &
843 t_firstcol=1 + nsg + 2*ntp + 1)
844
845 IF (debug_this_module) THEN
846 WRITE (log_unit, "(/,A)") "<T^+1|H_z|T^+1> in Hartree / cm^-1"
847 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
848 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
849 CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
850 END IF
851
852 !<T^-1|H_z|T^-1>, takes a minus sign
853 CALL cp_fm_scale(-1.0_dp, tmp_fm)
854 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
855 mtarget=img_fm, &
856 nrow=nex, &
857 ncol=nex, &
858 s_firstrow=1, &
859 s_firstcol=1, &
860 t_firstrow=1 + nsg + 1, &
861 t_firstcol=1 + nsg + 1)
862
863 IF (debug_this_module) THEN
864 WRITE (log_unit, "(/,A)") "<T^-1|H_z|T^-1> in Hartree / cm^-1"
865 CALL dbcsr_copy(dbcsr_dummy, dbcsr_tmp)
866 CALL dbcsr_scale(dbcsr_dummy, wavenumbers)
867 CALL dbcsr_print(dbcsr_dummy, unit_nr=log_unit)
868 END IF
869
870 !Intermediate clean-up
871 CALL cp_fm_struct_release(work_struct)
872 CALL cp_fm_release(work_fm)
873 CALL cp_fm_release(tmp_fm)
874 DEALLOCATE (diag, mo_soc_x, mo_soc_y, mo_soc_z)
875
876 IF (log_unit > 0) THEN
877 WRITE (log_unit, '(A)') "Diagonlzing SOC-Matrix"
878 END IF
879
880 !Set-up the complex hermitian perturbation matrix
881 CALL cp_cfm_create(hami_cfm, full_struct)
882 CALL cp_fm_to_cfm(real_fm, img_fm, hami_cfm)
883
884 !!Optinal Output: SOC-Matrix
885 IF (logger%para_env%is_source()) THEN
886 log_unit = cp_print_key_unit_nr(logger, &
887 tddfpt_print_section, &
888 "SOC_PRINT", &
889 extension=".socme", &
890 file_form="FORMATTED", &
891 file_action="WRITE", &
892 file_status="UNKNOWN")
893 ELSE
894 log_unit = -1
895 END IF
896
897 !! Get the requested energy unit and optional outputs
898 CALL section_vals_val_get(soc_print_section, "UNIT_eV", l_val=print_ev)
899 CALL section_vals_val_get(soc_print_section, "UNIT_wn", l_val=print_wn)
900 CALL section_vals_val_get(soc_print_section, "SPLITTING", l_val=print_splitting)
901 CALL section_vals_val_get(soc_print_section, "SOME", l_val=print_some)
902
903 !! save convertion unit into different varible for cleaner code
904 IF (print_wn) THEN
905 print_ev = .false.
906 unit_scale = wavenumbers
907 ELSE IF (print_ev) THEN
908 unit_scale = evolt
909 ELSE
910 cpwarn("No output unit was selected, printout will be in Hartree!")
911 unit_scale = 1
912 END IF
913
914 !! This needs something to deal with a parallel run
915 !! This output will be needed for NEWTONX for ISC!!
916 IF (print_some) THEN
917 IF (log_unit > 0) THEN
918 WRITE (log_unit, '(A)') "____________________________________________________"
919 WRITE (log_unit, '(A)') " FULL SOC-Matrix "
920 IF (print_ev) THEN
921 WRITE (log_unit, '(A)') "STATE STATE REAL-PART[eV] IMG-PART[eV] "
922 ELSE IF (print_wn) THEN
923 WRITE (log_unit, '(A)') "STATE STATE REAL-PART[cm-1] IMG-PART[cm-1]"
924 ELSE
925 WRITE (log_unit, '(A)') "STATE STATE REAL-PART[Hartree] IMG-PART[Hartree]"
926 END IF
927 WRITE (log_unit, '(A)') "____________________________________________________"
928 END IF
929 ALLOCATE (real_tmp(ntot, ntot), img_tmp(ntot, ntot))
930 real_tmp = 0_dp
931 img_tmp = 0_dp
932
933 CALL cp_fm_get_submatrix(real_fm, real_tmp, 1, 1, ntot, ntot)
934 CALL cp_fm_get_submatrix(img_fm, img_tmp, 1, 1, ntot, ntot)
935
936 IF (log_unit > 0) THEN
937 DO jj = 1, ntot
938 DO ii = 1, ntot
939 WRITE (unit=log_unit, fmt="(I3,5X,I3,5X,f16.8,5X,f16.8)") ii, jj, &
940 real_tmp(ii, jj)*unit_scale, &
941 img_tmp(ii, jj)*unit_scale
942 END DO !! jj
943 END DO !! ii
944 WRITE (log_unit, '(A)') " END SOC-MATRIX "
945 WRITE (log_unit, '(A)') "____________________________________________________"
946 WRITE (log_unit, '(A)') "____________________________________________________"
947 END IF !(log_unit)
948 DEALLOCATE (real_tmp, img_tmp)
949 END IF !print_some
950
951 CALL cp_fm_release(real_fm)
952 CALL cp_fm_release(img_fm)
953
954 ! Diagonalize the Hamiltonian
955 ALLOCATE (tmp_evals(ntot))
956 CALL cp_cfm_create(evecs_cfm, full_struct)
957 CALL cp_cfm_heevd(hami_cfm, evecs_cfm, tmp_evals)
958
959 ! Adjust the energies so the GS has zero, and store in the donor_state (without the GS)
960 ALLOCATE (soc_env%soc_evals(ntot - 1))
961 soc_env%soc_evals(:) = tmp_evals(2:ntot) - tmp_evals(1)
962
963 !! We may be interested in the ground state stabilisation energy
964 IF (logger%para_env%is_source()) THEN
965 log_unit = cp_logger_get_default_unit_nr(logger, local=.true.)
966 ELSE
967 log_unit = -1
968 END IF
969
970 IF (log_unit > 0) THEN
971 WRITE (log_unit, '(A27,6X,F18.10)') "Ground state stabilisation:", (soc_env%soc_evals(1)*unit_scale)
972 IF (debug_this_module) THEN
973 WRITE (log_unit, '(A)') "Calculation Dipole Moments"
974 END IF
975 END IF
976
977 !Compute the dipole oscillator strengths
978 soc_env%evals_a => evals_sing
979 soc_env%evals_b => evals_trip
980 CALL soc_dipol(qs_env, dbcsr_soc_package, soc_env, &
981 evecs_cfm, eps_filter, dipole_form, gs_mos, &
982 evects_sing, evects_trip)
983
984 !Create final output
985 !! Output unit is choosen by the keyword "UNIT_eV" and "UNIT_wn" in "SOC_PRINT" section
986 IF (logger%para_env%is_source()) THEN
987 log_unit = cp_logger_get_default_unit_nr(logger, local=.true.)
988 ELSE
989 log_unit = -1
990 END IF
991
992 IF (log_unit > 0) THEN
993 WRITE (log_unit, '(A)') "------------------------------------------------------------------------------"
994 WRITE (log_unit, '(A)') "- SOC CORRECTED SPECTRUM -"
995 WRITE (log_unit, '(A)') "------------------------------------------------------------------------------"
996 IF (print_ev) THEN
997 WRITE (log_unit, '(A)') "- STATE SOC-corrected exc. energies [eV] Oscillator strengths [a.u.] -"
998 ELSE IF (print_wn) THEN
999 WRITE (log_unit, '(A)') "- STATE SOC-corrected exc. energies [cm-1] Oscillator strengths [a.u.]-"
1000 ELSE
1001 WRITE (log_unit, '(A)') "- STATE SOC-corrected exc.energies [Hartree] Oscillator strengths [a.u.]-"
1002 END IF
1003 WRITE (log_unit, '(A)') "------------------------------------------------------------------------------"
1004 DO istate = 1, SIZE(soc_env%soc_evals)
1005 WRITE (log_unit, '(6X,I5,11X,2F16.5)') istate, soc_env%soc_evals(istate)*unit_scale, &
1006 soc_env%soc_osc(istate)
1007 END DO
1008 !! This is used for regtesting
1009 WRITE (log_unit, '(A,F18.8)') "TDDFT+SOC : CheckSum (eV) ", sum(soc_env%soc_evals)*evolt
1010 END IF
1011
1012 !! We may be interested in the differenz between the SOC and
1013 !! TDDFPT Eigenvalues
1014 !! Activated with the "SPLITTING" keyword in "SOC_PRINT" section
1015 IF (print_splitting .OR. debug_this_module) THEN
1016 CALL resort_evects(evecs_cfm, evects_sort)
1017 IF (log_unit > 0) THEN
1018 WRITE (log_unit, '(A)') "-----------------------------------------------------------------------------"
1019 WRITE (log_unit, '(A)') "- Splittings -"
1020 WRITE (log_unit, '(A)') "-----------------------------------------------------------------------------"
1021 IF (print_ev) THEN
1022 WRITE (log_unit, '(A)') "- STATE exc.energies[eV] splittings[eV] -"
1023 ELSE IF (print_wn) THEN
1024 WRITE (log_unit, '(A)') "- STATE exc.energies[cm-1] splittings[cm-1] -"
1025 ELSE
1026 WRITE (log_unit, '(A)') "- STATE exc.energies[Hartree] splittings[Hartree] -"
1027 END IF
1028 WRITE (log_unit, '(A)') "-----------------------------------------------------------------------------"
1029 mult = " "
1030 DO iex = 1, SIZE(soc_env%soc_evals)
1031 IF (evects_sort(iex + 1) .GT. nsg + ntp*2 + 1) THEN
1032 mult = "T+1"
1033 ELSE IF (evects_sort(iex + 1) .GT. nsg + ntp + 1) THEN
1034 mult = "T0"
1035 ELSE IF (evects_sort(iex + 1) .GT. nsg + 1) THEN
1036 mult = "T-1"
1037 ELSE
1038 mult = "S "
1039 END IF
1040 IF (mult == "T-1") THEN
1041 WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
1042 evects_sort(iex + 1) - (1 + nsg), soc_env%soc_evals(iex)*unit_scale, &
1043 (soc_env%soc_evals(iex) - evals_trip(evects_sort(iex + 1) - (1 + nsg)))*unit_scale
1044 ELSE IF (mult == "T0") THEN
1045 WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
1046 evects_sort(iex + 1) - (1 + nsg + ntp), soc_env%soc_evals(iex)*unit_scale, &
1047 (soc_env%soc_evals(iex) - evals_trip(evects_sort(iex + 1) - (1 + nsg + ntp)))*unit_scale
1048 ELSE IF (mult == "T+1") THEN
1049 WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
1050 evects_sort(iex + 1) - (1 + nsg + ntp*2), soc_env%soc_evals(iex)*unit_scale, &
1051 (soc_env%soc_evals(iex) - evals_trip(evects_sort(iex + 1) - (1 + nsg + ntp*2)))*unit_scale
1052 ELSE IF (mult == "S ") THEN
1053 WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
1054 evects_sort(iex + 1), soc_env%soc_evals(iex)*unit_scale, &
1055 (soc_env%soc_evals(iex) - evals_sing(evects_sort(iex + 1) - 1))*unit_scale
1056 END IF
1057 END DO
1058 WRITE (log_unit, '(A)') "----------------------------------------------------------------------------"
1059 DEALLOCATE (evects_sort)
1060 END IF
1061 END IF
1062
1063 !! clean up
1064 CALL soc_env_release(soc_env)
1065 CALL cp_fm_struct_release(full_struct)
1066 CALL cp_cfm_release(hami_cfm)
1067 CALL cp_cfm_release(evecs_cfm)
1068 CALL dbcsr_distribution_release(coeffs_dist)
1069 CALL dbcsr_distribution_release(prod_dist)
1070 CALL dbcsr_release(dbcsr_sg)
1071 CALL dbcsr_release(dbcsr_tp)
1072 CALL dbcsr_release(dbcsr_prod)
1073 CALL dbcsr_release(dbcsr_ovlp)
1074 CALL dbcsr_release(dbcsr_tmp)
1075 CALL dbcsr_release(dbcsr_work)
1076 IF (debug_this_module) THEN
1077 CALL dbcsr_release(dbcsr_dummy)
1078 DEALLOCATE (dbcsr_dummy)
1079 END IF
1080
1081 DEALLOCATE (dbcsr_work, dbcsr_prod, dbcsr_ovlp, dbcsr_tmp)
1082 DEALLOCATE (coeffs_dist, prod_dist, col_dist, col_blk_size, row_dist_new)
1083 DEALLOCATE (dbcsr_sg, dbcsr_tp, tmp_evals)
1084
1085 CALL timestop(handle)
1086
1087 END SUBROUTINE tddfpt_soc_rcs
1088
1089! **************************************************************************************************
1090!> \brief Some additional initalizations
1091!> \param qs_env Quickstep environment
1092!> \param soc_atom_env contains information to build the ZORA-Hamiltionian
1093!> \param soc_env contails details and outputs for SOC Correction
1094!> \param evects_a singlet Eigenvectors
1095!> \param evects_b triplet eigenvectors
1096!> \param dipole_form choosen dipole-form from TDDFPT-Section
1097!> \param gs_mos ...
1098! **************************************************************************************************
1099 SUBROUTINE inititialize_soc(qs_env, soc_atom_env, soc_env, &
1100 evects_a, evects_b, &
1101 dipole_form, gs_mos)
1102 !! Different Types
1103 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1104 TYPE(soc_atom_env_type), INTENT(OUT), POINTER :: soc_atom_env
1105 TYPE(soc_env_type), INTENT(OUT), TARGET :: soc_env
1106 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects_a, evects_b
1107 INTEGER :: dipole_form
1108 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1109 INTENT(in) :: gs_mos
1110
1111 CHARACTER(LEN=*), PARAMETER :: routinen = 'inititialize_soc'
1112
1113 CHARACTER(len=default_string_length), &
1114 ALLOCATABLE, DIMENSION(:, :) :: grid_info
1115 CHARACTER(len=default_string_length), &
1116 DIMENSION(:), POINTER :: k_list
1117 INTEGER :: handle, i, ikind, irep, ispin, nactive, &
1118 nao, natom, nkind, nrep, nspins, &
1119 nstates
1120 LOGICAL :: all_potential_present, do_os
1121 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1122 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1123 TYPE(cp_fm_type), POINTER :: a_coeff, b_coeff
1124 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1125 TYPE(dft_control_type), POINTER :: dft_control
1126 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1127 TYPE(mp_para_env_type), POINTER :: para_env
1128 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1129 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1130 TYPE(section_vals_type), POINTER :: soc_section
1131 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
1132
1133 CALL timeset(routinen, handle)
1134
1135 NULLIFY (qs_kind_set, particle_set, blacs_env, para_env, &
1136 a_coeff, b_coeff, tddfpt_control, soc_section)
1137
1138 !! Open_shell is not implemented yet
1139 do_os = .false.
1140
1141 !! soc_env will pass most of the larger matrices to the different modules
1142 CALL soc_env_create(soc_env)
1143 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, &
1144 mos=mos, natom=natom, &
1145 particle_set=particle_set, blacs_env=blacs_env, &
1146 para_env=para_env, matrix_s=matrix_s, &
1147 dft_control=dft_control)
1148
1149 !! The Dipole form shall be consistent with the tddfpt-module!
1150 tddfpt_control => dft_control%tddfpt2_control
1151 dipole_form = tddfpt_control%dipole_form
1152
1153 soc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%SOC")
1154
1155 !! We may want to use a bigger grid then default
1156 CALL section_vals_val_get(soc_section, "GRID", n_rep_val=nrep)
1157 ALLOCATE (grid_info(nrep, 3))
1158 DO irep = 1, nrep
1159 CALL section_vals_val_get(soc_section, &
1160 "GRID", &
1161 i_rep_val=irep, &
1162 c_vals=k_list)
1163 grid_info(irep, :) = k_list
1164 END DO
1165
1166 CALL soc_dipole_operator(soc_env, tddfpt_control, qs_env, gs_mos)
1167
1168 nspins = SIZE(evects_a, 1)
1169 DO ispin = 1, nspins
1170 CALL cp_fm_get_info(evects_a(ispin, 1), nrow_global=nao, ncol_global=nactive)
1171 END DO
1172
1173 !! We need to change and create structures for the exitation vectors.
1174 !! All excited states shall be written in one matrix, oneafter the other
1175 nstates = SIZE(evects_a, 2)
1176
1177 a_coeff => soc_env%a_coeff
1178 b_coeff => soc_env%b_coeff
1179
1180 NULLIFY (fm_struct)
1181 CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
1182 nrow_global=nao, ncol_global=nspins*nactive*nstates)
1183 CALL cp_fm_create(a_coeff, fm_struct)
1184 CALL cp_fm_create(b_coeff, fm_struct)
1185 CALL cp_fm_struct_release(fm_struct)
1186
1187 CALL soc_contract_evect(evects_a, a_coeff)
1188 CALL soc_contract_evect(evects_b, b_coeff)
1189
1190 ALLOCATE (soc_env%orb_soc(3))
1191 DO i = 1, 3
1192 ALLOCATE (soc_env%orb_soc(i)%matrix)
1193 END DO
1194
1195 ! Make soc_atom_env for H^SOC calculations
1196 ! Compute the contraction coefficients for spherical orbitals
1197 CALL soc_atom_create(soc_atom_env)
1198 IF (do_os) THEN
1199 soc_atom_env%nspins = 2
1200 ELSE
1201 soc_atom_env%nspins = 1
1202 END IF
1203
1204 nkind = SIZE(qs_kind_set)
1205 ALLOCATE (soc_atom_env%grid_atom_set(nkind))
1206 ALLOCATE (soc_atom_env%harmonics_atom_set(nkind))
1207 ALLOCATE (soc_atom_env%orb_sphi_so(nkind))
1208
1209 CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present)
1210 IF (.NOT. all_potential_present) THEN
1211 CALL v_soc_xyz_from_pseudopotential(qs_env, soc_atom_env%soc_pp)
1212 END IF
1213
1214 !! Prepare the atomic grid we need to calculate the SOC Operator in an Atomic basis
1215 !! copied from init_xas_atom_grid_harmo for convenience
1216 CALL init_atom_grid(soc_atom_env, grid_info, qs_env)
1217
1218 DO ikind = 1, nkind
1219 NULLIFY (soc_atom_env%orb_sphi_so(ikind)%array)
1220 CALL compute_sphi_so(ikind, "ORB", soc_atom_env%orb_sphi_so(ikind)%array, qs_env)
1221 END DO !ikind
1222
1223 DEALLOCATE (grid_info)
1224
1225 CALL timestop(handle)
1226
1227 END SUBROUTINE inititialize_soc
1228
1229! **************************************************************************************************
1230!> \brief Initializes the atomic grids and harmonics for the atomic calculations
1231!> \param soc_atom_env ...
1232!> \param grid_info ...
1233!> \param qs_env ...
1234!> \note Copied and modified from init_xas_atom_grid_harmo
1235! **************************************************************************************************
1236 SUBROUTINE init_atom_grid(soc_atom_env, grid_info, qs_env)
1237
1238 TYPE(soc_atom_env_type) :: soc_atom_env
1239 CHARACTER(len=default_string_length), &
1240 ALLOCATABLE, DIMENSION(:, :) :: grid_info
1241 TYPE(qs_environment_type) :: qs_env
1242
1243 CHARACTER(LEN=*), PARAMETER :: routinen = 'init_atom_grid'
1244
1245 CHARACTER(LEN=default_string_length) :: kind_name
1246 INTEGER :: handle, igrid, ikind, il, iso, iso1, iso2, l1, l1l2, l2, la, lc1, lc2, lcleb, ll, &
1247 llmax, lp, m1, m2, max_s_harm, max_s_set, maxl, maxlgto, maxs, mm, mp, na, nr, &
1248 quadrature, stat
1249 REAL(dp) :: kind_radius
1250 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga
1251 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
1252 TYPE(dft_control_type), POINTER :: dft_control
1253 TYPE(grid_atom_type), POINTER :: grid_atom
1254 TYPE(gto_basis_set_type), POINTER :: tmp_basis
1255 TYPE(harmonics_atom_type), POINTER :: harmonics
1256 TYPE(qs_control_type), POINTER :: qs_control
1257 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1258
1259 CALL timeset(routinen, handle)
1260
1261 NULLIFY (my_cg, qs_kind_set, dft_control, qs_control)
1262 NULLIFY (grid_atom, harmonics, tmp_basis)
1263
1264 !! Initialization of some integer for the CG coeff generation
1265 CALL get_qs_env(qs_env, &
1266 qs_kind_set=qs_kind_set, &
1267 dft_control=dft_control)
1268
1269 CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="ORB")
1270 qs_control => dft_control%qs_control
1271
1272 !!To be sure:: Is the basis grid present?
1273 IF (maxlgto < 0) THEN
1274 cpabort("tmp_basis is not Present!")
1275 END IF
1276
1277 !maximum expansion
1278 llmax = 2*maxlgto
1279 max_s_harm = nsoset(llmax)
1280 max_s_set = nsoset(maxlgto)
1281 lcleb = llmax
1282
1283 ! Allocate and compute the CG coeffs (copied from init_rho_atom)
1284 CALL clebsch_gordon_init(lcleb)
1285 CALL reallocate(my_cg, 1, max_s_set, 1, max_s_set, 1, max_s_harm)
1286
1287 ALLOCATE (rga(lcleb, 2))
1288 DO lc1 = 0, maxlgto
1289 DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
1290 l1 = indso(1, iso1)
1291 m1 = indso(2, iso1)
1292 DO lc2 = 0, maxlgto
1293 DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
1294 l2 = indso(1, iso2)
1295 m2 = indso(2, iso2)
1296 CALL clebsch_gordon(l1, m1, l2, m2, rga)
1297 IF (l1 + l2 > llmax) THEN
1298 l1l2 = llmax
1299 ELSE
1300 l1l2 = l1 + l2
1301 END IF
1302 mp = m1 + m2
1303 mm = m1 - m2
1304 IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
1305 mp = -abs(mp)
1306 mm = -abs(mm)
1307 ELSE
1308 mp = abs(mp)
1309 mm = abs(mm)
1310 END IF
1311 DO lp = mod(l1 + l2, 2), l1l2, 2
1312 il = lp/2 + 1
1313 IF (abs(mp) <= lp) THEN
1314 IF (mp >= 0) THEN
1315 iso = nsoset(lp - 1) + lp + 1 + mp
1316 ELSE
1317 iso = nsoset(lp - 1) + lp + 1 - abs(mp)
1318 END IF
1319 my_cg(iso1, iso2, iso) = rga(il, 1)
1320 END IF
1321 IF (mp /= mm .AND. abs(mm) <= lp) THEN
1322 IF (mm >= 0) THEN
1323 iso = nsoset(lp - 1) + lp + 1 + mm
1324 ELSE
1325 iso = nsoset(lp - 1) + lp + 1 - abs(mm)
1326 END IF
1327 my_cg(iso1, iso2, iso) = rga(il, 2)
1328 END IF
1329 END DO
1330 END DO ! iso2
1331 END DO ! lc2
1332 END DO ! iso1
1333 END DO ! lc1
1334 DEALLOCATE (rga)
1335 CALL clebsch_gordon_deallocate()
1336
1337 ! Create the Lebedev grids and compute the spherical harmonics
1338 CALL init_lebedev_grids()
1339 quadrature = qs_control%gapw_control%quadrature
1340
1341 DO ikind = 1, SIZE(soc_atom_env%grid_atom_set)
1342
1343 ! Allocate the grid and the harmonics for this kind
1344 NULLIFY (soc_atom_env%grid_atom_set(ikind)%grid_atom)
1345 NULLIFY (soc_atom_env%harmonics_atom_set(ikind)%harmonics_atom)
1346 CALL allocate_grid_atom(soc_atom_env%grid_atom_set(ikind)%grid_atom)
1347 CALL allocate_harmonics_atom( &
1348 soc_atom_env%harmonics_atom_set(ikind)%harmonics_atom)
1349
1350 NULLIFY (grid_atom, harmonics)
1351 grid_atom => soc_atom_env%grid_atom_set(ikind)%grid_atom
1352 harmonics => soc_atom_env%harmonics_atom_set(ikind)%harmonics_atom
1353
1354 ! Initialize some integers
1355 CALL get_qs_kind(qs_kind_set(ikind), &
1356 ngrid_rad=nr, &
1357 ngrid_ang=na, &
1358 name=kind_name)
1359
1360 ! take the grid dimension given as input, if none, take the GAPW ones above
1361 DO igrid = 1, SIZE(grid_info, 1)
1362 IF (grid_info(igrid, 1) == kind_name) THEN
1363 READ (grid_info(igrid, 2), *, iostat=stat) na
1364 IF (stat .NE. 0) cpabort("The 'na' value for the GRID keyword must be an integer")
1365 READ (grid_info(igrid, 3), *, iostat=stat) nr
1366 IF (stat .NE. 0) cpabort("The 'nr' value for the GRID keyword must be an integer")
1367 EXIT
1368 END IF
1369 END DO !igrid
1370
1371 ll = get_number_of_lebedev_grid(n=na)
1372 na = lebedev_grid(ll)%n
1373 la = lebedev_grid(ll)%l
1374 grid_atom%ng_sphere = na
1375 grid_atom%nr = nr
1376
1377 !Create the harmonics with the ORB-Basis
1378 CALL get_qs_kind(qs_kind_set(ikind), &
1379 basis_set=tmp_basis, &
1380 basis_type="ORB")
1381 CALL get_gto_basis_set(gto_basis_set=tmp_basis, &
1382 maxl=maxl, &
1383 kind_radius=kind_radius)
1384
1385 CALL create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
1386 CALL truncate_radial_grid(grid_atom, kind_radius)
1387
1388 maxs = nsoset(maxl)
1389 CALL create_harmonics_atom(harmonics, &
1390 my_cg, na, &
1391 llmax, maxs, &
1392 max_s_harm, ll, &
1393 grid_atom%wa, &
1394 grid_atom%azi, &
1395 grid_atom%pol)
1396 CALL get_maxl_cg(harmonics, tmp_basis, llmax, max_s_harm)
1397 END DO !ikind
1398
1399 CALL deallocate_lebedev_grids()
1400 DEALLOCATE (my_cg)
1401
1402 CALL timestop(handle)
1403
1404 END SUBROUTINE init_atom_grid
1405
1406! **************************************************************************************************
1407!> \brief ...
1408!> \param qs_env ...
1409!> \param dbcsr_soc_package ...
1410!> \param soc_env ...
1411!> \param soc_evecs_cfm ...
1412!> \param eps_filter ...
1413!> \param dipole_form ...
1414!> \param gs_mos ...
1415!> \param evects_sing ...
1416!> \param evects_trip ...
1417! **************************************************************************************************
1418 SUBROUTINE soc_dipol(qs_env, dbcsr_soc_package, soc_env, &
1419 soc_evecs_cfm, eps_filter, dipole_form, gs_mos, &
1420 evects_sing, evects_trip)
1421 TYPE(qs_environment_type), POINTER :: qs_env
1422 TYPE(dbcsr_soc_package_type) :: dbcsr_soc_package
1423 TYPE(soc_env_type), TARGET :: soc_env
1424 TYPE(cp_cfm_type), INTENT(in) :: soc_evecs_cfm
1425 REAL(dp), INTENT(IN) :: eps_filter
1426 INTEGER :: dipole_form
1427 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1428 POINTER :: gs_mos
1429 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects_sing, evects_trip
1430
1431 CHARACTER(len=*), PARAMETER :: routinen = 'soc_dipol'
1432
1433 COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: transdip
1434 INTEGER :: handle, i, nosc, ntot
1435 REAL(dp), DIMENSION(:), POINTER :: osc_str, soc_evals
1436 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1437 TYPE(cp_cfm_type) :: dip_cfm, work1_cfm, work2_cfm
1438 TYPE(cp_fm_struct_type), POINTER :: dip_struct, full_struct
1439 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: amew_dip
1440 TYPE(mp_para_env_type), POINTER :: para_env
1441
1442 CALL timeset(routinen, handle)
1443
1444 NULLIFY (para_env, blacs_env, dip_struct, full_struct, osc_str)
1445 NULLIFY (soc_evals)
1446
1447 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1448
1449 soc_evals => soc_env%soc_evals
1450 nosc = SIZE(soc_evals)
1451 ntot = nosc + 1
1452 ALLOCATE (soc_env%soc_osc(nosc))
1453 osc_str => soc_env%soc_osc
1454 osc_str(:) = 0.0_dp
1455
1456 !get some work arrays/matrix
1457 CALL cp_fm_struct_create(dip_struct, &
1458 context=blacs_env, para_env=para_env, &
1459 nrow_global=ntot, ncol_global=1)
1460
1461 CALL cp_cfm_get_info(soc_evecs_cfm, matrix_struct=full_struct)
1462 CALL cp_cfm_create(dip_cfm, dip_struct)
1463 CALL cp_cfm_create(work1_cfm, full_struct)
1464 CALL cp_cfm_create(work2_cfm, full_struct)
1465
1466 ALLOCATE (transdip(ntot, 1))
1467
1468 CALL get_rcs_amew_op(amew_dip, soc_env, dbcsr_soc_package, &
1469 eps_filter, qs_env, gs_mos, &
1470 evects_sing, evects_trip)
1471
1472 DO i = 1, 3 !cartesian coord x, y, z
1473
1474 !Convert the real dipole into the cfm format for calculations
1475 CALL cp_fm_to_cfm(msourcer=amew_dip(i), mtarget=work1_cfm)
1476
1477 !compute amew_coeffs^dagger * amew_dip * amew_gs to get the transition moments
1478 CALL parallel_gemm('C', 'N', ntot, &
1479 ntot, ntot, &
1480 (1.0_dp, 0.0_dp), &
1481 soc_evecs_cfm, &
1482 work1_cfm, &
1483 (0.0_dp, 0.0_dp), &
1484 work2_cfm)
1485 CALL parallel_gemm('N', 'N', ntot, &
1486 1, ntot, &
1487 (1.0_dp, 0.0_dp), &
1488 work2_cfm, &
1489 soc_evecs_cfm, &
1490 (0.0_dp, 0.0_dp), &
1491 dip_cfm)
1492
1493 CALL cp_cfm_get_submatrix(dip_cfm, transdip)
1494
1495 !transition dipoles are real numbers
1496 osc_str(:) = osc_str(:) + real(transdip(2:ntot, 1))**2 + aimag(transdip(2:ntot, 1))**2
1497
1498 END DO !i
1499
1500 !multiply with appropriate prefac depending in the rep
1501 IF (dipole_form == tddfpt_dipole_length) THEN
1502 osc_str(:) = 2.0_dp/3.0_dp*soc_evals(:)*osc_str(:)
1503 ELSE
1504 osc_str(:) = 2.0_dp/3.0_dp*soc_evals(:)*osc_str(:)
1505 END IF
1506
1507 !clean-up
1508 CALL cp_fm_struct_release(dip_struct)
1509 CALL cp_cfm_release(work1_cfm)
1510 CALL cp_cfm_release(work2_cfm)
1511 CALL cp_cfm_release(dip_cfm)
1512 DO i = 1, 3
1513 CALL cp_fm_release(amew_dip(i))
1514 END DO
1515 DEALLOCATE (amew_dip, transdip)
1516
1517 CALL timestop(handle)
1518
1519 END SUBROUTINE soc_dipol
1520
1521!**********************************************************************************
1522!> \brief: This will create the Dipole operator within the amew basis
1523!> \param amew_op Output dipole operator
1524!> \param soc_env ...
1525!> \param dbcsr_soc_package Some work matrices
1526!> \param eps_filter ...
1527!> \param qs_env ...
1528!> \param gs_mos ...
1529!> \param evects_sing ...
1530!> \param evects_trip ...
1531! **************************************************************************************************
1532 SUBROUTINE get_rcs_amew_op(amew_op, soc_env, dbcsr_soc_package, eps_filter, qs_env, gs_mos, &
1533 evects_sing, evects_trip)
1534
1535 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
1536 INTENT(OUT) :: amew_op
1537 TYPE(soc_env_type), TARGET :: soc_env
1538 TYPE(dbcsr_soc_package_type) :: dbcsr_soc_package
1539 REAL(dp), INTENT(IN) :: eps_filter
1540 TYPE(qs_environment_type), POINTER :: qs_env
1541 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1542 POINTER :: gs_mos
1543 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects_sing, evects_trip
1544
1545 CHARACTER(len=*), PARAMETER :: routinen = 'get_rcs_amew_op'
1546
1547 INTEGER :: dim_op, handle, i, isg, nactive, nao, &
1548 nex, nsg, ntot, ntp
1549 LOGICAL :: vel_rep
1550 REAL(dp) :: op, sqrt2
1551 REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag, gs_diag, gsgs_op
1552 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: mo_op, sggs_block
1553 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1554 TYPE(cp_fm_struct_type), POINTER :: full_struct, gsgs_struct, sggs_struct, &
1555 std_struct, tmp_struct
1556 TYPE(cp_fm_type) :: gs_fm, sggs_fm, tmp_fm, vec_op, work_fm
1557 TYPE(cp_fm_type), POINTER :: gs_coeffs
1558 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ao_op, matrix_s
1559 TYPE(dbcsr_type), POINTER :: ao_op_i, dbcsr_ovlp, dbcsr_prod, &
1560 dbcsr_sg, dbcsr_tmp, dbcsr_tp, &
1561 dbcsr_work
1562 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1563 TYPE(mp_para_env_type), POINTER :: para_env
1564
1565 CALL timeset(routinen, handle)
1566
1567 NULLIFY (mos, gs_coeffs)
1568 NULLIFY (matrix_s)
1569 NULLIFY (para_env, blacs_env)
1570 NULLIFY (full_struct, gsgs_struct, std_struct, tmp_struct, sggs_struct)
1571 NULLIFY (ao_op_i, ao_op)
1572 NULLIFY (dbcsr_tp, dbcsr_sg, dbcsr_ovlp, dbcsr_work, dbcsr_tmp, dbcsr_prod)
1573
1574 ! Initialization
1575 !! Fist case for length, secound for velocity-representation
1576 IF (ASSOCIATED(soc_env%dipmat)) THEN
1577 ao_op => soc_env%dipmat
1578 vel_rep = .false.
1579 ELSE
1580 ao_op => soc_env%dipmat_ao
1581 vel_rep = .true.
1582 END IF
1583 nsg = SIZE(soc_env%evals_a)
1584 ntp = nsg; nex = nsg !all the same by construction, keep them separate for clarity
1585 ntot = 1 + nsg + 3*ntp
1586
1587 CALL get_qs_env(qs_env, matrix_s=matrix_s, mos=mos, para_env=para_env, blacs_env=blacs_env)
1588 sqrt2 = sqrt(2.0_dp)
1589 dim_op = SIZE(ao_op)
1590
1591 dbcsr_sg => dbcsr_soc_package%dbcsr_sg
1592 dbcsr_tp => dbcsr_soc_package%dbcsr_tp
1593 dbcsr_work => dbcsr_soc_package%dbcsr_work
1594 dbcsr_prod => dbcsr_soc_package%dbcsr_prod
1595 dbcsr_ovlp => dbcsr_soc_package%dbcsr_ovlp
1596 dbcsr_tmp => dbcsr_soc_package%dbcsr_tmp
1597
1598 ! Create the amew_op matrix
1599 CALL cp_fm_struct_create(full_struct, &
1600 context=blacs_env, para_env=para_env, &
1601 nrow_global=ntot, ncol_global=ntot)
1602 ALLOCATE (amew_op(dim_op))
1603 DO i = 1, dim_op
1604 CALL cp_fm_create(amew_op(i), full_struct)
1605 END DO !i
1606
1607 ! Deal with the GS-GS contribution <0|0> = 2*sum_j <phi_j|op|phi_j>
1608 CALL get_mo_set(mos(1), nao=nao, homo=nactive)
1609 gs_coeffs => gs_mos(1)%mos_occ
1610 CALL cp_fm_get_info(gs_coeffs, matrix_struct=std_struct)
1611
1612 CALL cp_fm_struct_create(gsgs_struct, &
1613 context=blacs_env, para_env=para_env, &
1614 nrow_global=nactive, ncol_global=nactive)
1615
1616 CALL cp_fm_create(gs_fm, gsgs_struct) !nactive nactive
1617 CALL cp_fm_create(work_fm, std_struct) !nao nactive
1618
1619 ALLOCATE (gsgs_op(dim_op))
1620 ALLOCATE (gs_diag(nactive))
1621
1622 DO i = 1, dim_op
1623
1624 ao_op_i => ao_op(i)%matrix
1625 CALL cp_dbcsr_sm_fm_multiply(ao_op_i, gs_coeffs, work_fm, ncol=nactive)
1626 CALL parallel_gemm('T', 'N', nactive, nactive, nao, &
1627 1.0_dp, gs_coeffs, work_fm, 0.0_dp, gs_fm)
1628 CALL cp_fm_get_diag(gs_fm, gs_diag)
1629 gsgs_op(i) = 2.0_dp*sum(gs_diag)
1630
1631 END DO !i
1632 DEALLOCATE (gs_diag)
1633
1634 CALL cp_fm_release(work_fm)
1635
1636 ! Create the work and helper fms
1637 CALL cp_fm_create(vec_op, std_struct) !nao nactive
1638 CALL cp_fm_struct_release(gsgs_struct)
1639
1640 CALL cp_fm_struct_create(tmp_struct, &
1641 context=blacs_env, para_env=para_env, &
1642 nrow_global=nex, ncol_global=nex)
1643 CALL cp_fm_struct_create(sggs_struct, &
1644 context=blacs_env, para_env=para_env, &
1645 nrow_global=nactive*nsg, ncol_global=nactive)
1646
1647 CALL cp_fm_create(tmp_fm, tmp_struct)
1648 CALL cp_fm_create(work_fm, full_struct)
1649 CALL cp_fm_create(sggs_fm, sggs_struct)
1650
1651 ALLOCATE (diag(nactive))
1652 ALLOCATE (mo_op(nactive, nactive))
1653 ALLOCATE (sggs_block(nactive, nactive))
1654
1655 ! Iterate over the dimensions of the operator
1656 ! Note: operator matrices are asusmed symmetric, can only do upper half
1657 DO i = 1, dim_op
1658
1659 ao_op_i => ao_op(i)%matrix
1660
1661 ! The GS-GS contribution
1662 CALL cp_fm_set_element(amew_op(i), 1, 1, gsgs_op(i))
1663
1664 ! Compute the operator in the MO basis
1665 CALL cp_dbcsr_sm_fm_multiply(ao_op_i, gs_coeffs, vec_op, ncol=nactive)
1666 CALL cp_fm_get_submatrix(gs_fm, mo_op) ! instead of prod_fm
1667
1668 ! Compute the ground-state/singlet components. ao_op*gs_coeffs already stored in vec_op
1669 IF (vel_rep) THEN
1670 CALL dip_vel_op(soc_env, qs_env, evects_sing, dbcsr_prod, i, .false., &
1671 gs_coeffs, sggs_fm)
1672 ELSE
1673 CALL parallel_gemm('T', 'N', nactive*nsg, nactive, nao, &
1674 1.0_dp, soc_env%a_coeff, vec_op, 0.0_dp, sggs_fm)
1675 END IF
1676
1677 DO isg = 1, nsg
1678 CALL cp_fm_get_submatrix(fm=sggs_fm, &
1679 target_m=sggs_block, &
1680 start_row=(isg - 1)*nactive + 1, &
1681 start_col=1, &
1682 n_rows=nactive, &
1683 n_cols=nactive)
1684 diag(:) = get_diag(sggs_block)
1685 op = sqrt2*sum(diag)
1686 CALL cp_fm_set_element(amew_op(i), 1, 1 + isg, op)
1687 END DO !isg
1688
1689 ! do the singlet-singlet components
1690 !start with the overlap
1691 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
1692 matrix_s(1)%matrix, dbcsr_sg, &
1693 0.0_dp, dbcsr_work, &
1694 filter_eps=eps_filter)
1695 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
1696 dbcsr_sg, dbcsr_work, &
1697 0.0_dp, dbcsr_ovlp, &
1698 filter_eps=eps_filter)
1699
1700 IF (vel_rep) THEN
1701 CALL dip_vel_op(soc_env, qs_env, evects_sing, dbcsr_prod, i, .false.)
1702 ELSE
1703 !then the operator in the LR orbital basis
1704 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
1705 ao_op_i, dbcsr_sg, &
1706 0.0_dp, dbcsr_work, &
1707 filter_eps=eps_filter)
1708 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
1709 dbcsr_sg, dbcsr_work, &
1710 0.0_dp, dbcsr_prod, &
1711 filter_eps=eps_filter)
1712 END IF
1713
1714 !use the soc routine, it is compatible
1715 CALL rcs_amew_soc_elements(dbcsr_tmp, &
1716 dbcsr_prod, &
1717 dbcsr_ovlp, &
1718 mo_op, &
1719 pref_trace=-1.0_dp, &
1720 pref_overall=1.0_dp, &
1721 pref_diags=gsgs_op(i), &
1722 symmetric=.true.)
1723
1724 CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
1725 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
1726 mtarget=amew_op(i), &
1727 nrow=nex, &
1728 ncol=nex, &
1729 s_firstrow=1, &
1730 s_firstcol=1, &
1731 t_firstrow=2, &
1732 t_firstcol=2)
1733
1734 ! compute the triplet-triplet components
1735 !the overlap
1736 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
1737 matrix_s(1)%matrix, &
1738 dbcsr_tp, 0.0_dp, &
1739 dbcsr_work, &
1740 filter_eps=eps_filter)
1741 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
1742 dbcsr_tp, dbcsr_work, &
1743 0.0_dp, dbcsr_ovlp, &
1744 filter_eps=eps_filter)
1745
1746 !the operator in the LR orbital basis
1747 IF (vel_rep) THEN
1748 CALL dip_vel_op(soc_env, qs_env, evects_trip, dbcsr_prod, i, .true.)
1749 ELSE
1750 CALL dbcsr_multiply('N', 'N', 1.0_dp, &
1751 ao_op_i, dbcsr_tp, &
1752 0.0_dp, dbcsr_work, &
1753 filter_eps=eps_filter)
1754 CALL dbcsr_multiply('T', 'N', 1.0_dp, &
1755 dbcsr_tp, dbcsr_work, &
1756 0.0_dp, dbcsr_prod, &
1757 filter_eps=eps_filter)
1758 END IF
1759
1760 CALL rcs_amew_soc_elements(dbcsr_tmp, &
1761 dbcsr_prod, &
1762 dbcsr_ovlp, &
1763 mo_op, &
1764 pref_trace=-1.0_dp, &
1765 pref_overall=1.0_dp, &
1766 pref_diags=gsgs_op(i), &
1767 symmetric=.true.)
1768
1769 CALL copy_dbcsr_to_fm(dbcsr_tmp, tmp_fm)
1770 !<T^-1|op|T^-1>
1771 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
1772 mtarget=amew_op(i), &
1773 nrow=nex, &
1774 ncol=nex, &
1775 s_firstrow=1, &
1776 s_firstcol=1, &
1777 t_firstrow=1 + nsg + 1, &
1778 t_firstcol=1 + nsg + 1)
1779
1780 !<T^0|op|T^0>
1781 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
1782 mtarget=amew_op(i), &
1783 nrow=nex, &
1784 ncol=nex, &
1785 s_firstrow=1, &
1786 s_firstcol=1, &
1787 t_firstrow=1 + nsg + ntp + 1, &
1788 t_firstcol=1 + nsg + ntp + 1)
1789 !<T^+1|op|T^+1>
1790 CALL cp_fm_to_fm_submat(msource=tmp_fm, &
1791 mtarget=amew_op(i), &
1792 nrow=nex, &
1793 ncol=nex, &
1794 s_firstrow=1, &
1795 s_firstcol=1, &
1796 t_firstrow=1 + nsg + 2*ntp + 1, &
1797 t_firstcol=1 + nsg + 2*ntp + 1)
1798
1799 ! Symmetrize the matrix (only upper triangle built)
1800 CALL cp_fm_uplo_to_full(amew_op(i), work_fm)
1801
1802 END DO !i
1803
1804 ! Clean-up
1805 CALL cp_fm_release(gs_fm)
1806 CALL cp_fm_release(work_fm)
1807 CALL cp_fm_release(tmp_fm)
1808 CALL cp_fm_release(vec_op)
1809 CALL cp_fm_release(sggs_fm)
1810 CALL cp_fm_struct_release(full_struct)
1811 CALL cp_fm_struct_release(tmp_struct)
1812 CALL cp_fm_struct_release(sggs_struct)
1813
1814 DEALLOCATE (sggs_block)
1815 DEALLOCATE (diag, gsgs_op, mo_op)
1816 CALL timestop(handle)
1817
1818 END SUBROUTINE get_rcs_amew_op
1819
1820END MODULE qs_tddfpt2_soc
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
memory_utilities::reallocate
Definition memory_utilities.F:27
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
spherical_harmonics::clebsch_gordon
Definition spherical_harmonics.F:53
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Definition basis_set_types.F:624
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:58
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_get_submatrix
subroutine, public cp_cfm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
Extract a sub-matrix from the full matrix: op(target_m)(1:n_rows,1:n_cols) = fm(start_row:start_row+n...
Definition cp_cfm_types.F:329
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:155
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:813
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:603
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:629
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1165
cp_dbcsr_api::dbcsr_distribution_new
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
Definition cp_dbcsr_api.F:608
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1073
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:794
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:564
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_print
subroutine, public dbcsr_print(matrix, variable_name, unit_nr)
Prints given matrix in matlab format (only present blocks).
Definition cp_dbcsr_contrib.F:682
cp_dbcsr_contrib::dbcsr_reserve_all_blocks
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
Definition cp_dbcsr_contrib.F:295
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:552
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:214
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition cp_fm_basic_linalg.F:1695
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1747
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1665
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:563
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1009
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition cp_fm_types.F:1472
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:528
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:894
cp_fm_types::cp_fm_set_element
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
Definition cp_fm_types.F:693
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
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_dipole_length
integer, parameter, public tddfpt_dipole_length
Definition input_constants.F:700
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
lebedev
Generation of the spherical Lebedev grids. All Lebedev grids were generated with a precision of at le...
Definition lebedev.F:57
lebedev::deallocate_lebedev_grids
subroutine, public deallocate_lebedev_grids()
...
Definition lebedev.F:324
lebedev::lebedev_grid
type(oh_grid), dimension(nlg), target, public lebedev_grid
Definition lebedev.F:85
lebedev::get_number_of_lebedev_grid
integer function, public get_number_of_lebedev_grid(l, n)
Get the number of the Lebedev grid, which has the requested angular momentum quantnum number l or siz...
Definition lebedev.F:114
lebedev::init_lebedev_grids
subroutine, public init_lebedev_grids()
Load the coordinates and weights of the nonredundant Lebedev grid points.
Definition lebedev.F:344
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::get_diag
pure real(kind=dp) function, dimension(min(size(a, 1), size(a, 2))), public get_diag(a)
Return the diagonal elements of matrix a as a vector.
Definition mathlib.F:493
mathlib::diag
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
Definition mathlib.F:1496
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition orbital_pointers.F:42
orbital_pointers::indso
integer, dimension(:, :), allocatable, public indso
Definition orbital_pointers.F:43
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
physcon::wavenumbers
real(kind=dp), parameter, public wavenumbers
Definition physcon.F:192
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_grid_atom
Definition qs_grid_atom.F:8
qs_grid_atom::allocate_grid_atom
subroutine, public allocate_grid_atom(grid_atom)
Initialize components of the grid_atom_type structure.
Definition qs_grid_atom.F:73
qs_grid_atom::create_grid_atom
subroutine, public create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
Definition qs_grid_atom.F:183
qs_harmonics_atom
Definition qs_harmonics_atom.F:8
qs_harmonics_atom::allocate_harmonics_atom
subroutine, public allocate_harmonics_atom(harmonics)
Allocate a spherical harmonics set for the atom grid.
Definition qs_harmonics_atom.F:66
qs_harmonics_atom::get_maxl_cg
subroutine, public get_maxl_cg(harmonics, orb_basis, llmax, max_s_harm)
...
Definition qs_harmonics_atom.F:403
qs_harmonics_atom::create_harmonics_atom
subroutine, public create_harmonics_atom(harmonics, my_cg, na, llmax, maxs, max_s_harm, ll, wa, azi, pol)
...
Definition qs_harmonics_atom.F:158
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:880
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_tddfpt2_soc_types
Definition qs_tddfpt2_soc_types.F:7
qs_tddfpt2_soc_types::soc_env_create
subroutine, public soc_env_create(soc_env)
...
Definition qs_tddfpt2_soc_types.F:87
qs_tddfpt2_soc_types::soc_atom_release
subroutine, public soc_atom_release(soc_atom_env)
...
Definition qs_tddfpt2_soc_types.F:176
qs_tddfpt2_soc_types::soc_atom_create
subroutine, public soc_atom_create(soc_atom_env)
...
Definition qs_tddfpt2_soc_types.F:159
qs_tddfpt2_soc_types::soc_env_release
subroutine, public soc_env_release(soc_env)
...
Definition qs_tddfpt2_soc_types.F:101
qs_tddfpt2_soc_utils
Utilities absorption spectroscopy using TDDFPT with SOC.
Definition qs_tddfpt2_soc_utils.F:13
qs_tddfpt2_soc_utils::soc_dipole_operator
subroutine, public soc_dipole_operator(soc_env, tddfpt_control, qs_env, gs_mos)
Build the atomic dipole operator.
Definition qs_tddfpt2_soc_utils.F:91
qs_tddfpt2_soc_utils::resort_evects
subroutine, public resort_evects(evects_cfm, sort)
Used to find out, which state has which spin-multiplicity.
Definition qs_tddfpt2_soc_utils.F:518
qs_tddfpt2_soc_utils::soc_contract_evect
subroutine, public soc_contract_evect(fm_start, fm_res)
...
Definition qs_tddfpt2_soc_utils.F:452
qs_tddfpt2_soc_utils::dip_vel_op
subroutine, public dip_vel_op(soc_env, qs_env, evec_fm, op, ideriv, tp, gs_coeffs, sggs_fm)
This routine will construct the dipol operator within velocity representation.
Definition qs_tddfpt2_soc_utils.F:350
qs_tddfpt2_soc
Definition qs_tddfpt2_soc.F:7
qs_tddfpt2_soc::tddfpt_soc
subroutine, public tddfpt_soc(qs_env, evals_a, evals_b, evects_a, evects_b, gs_mos)
Perform TDDFPT-SOC calculation.
Definition qs_tddfpt2_soc.F:141
qs_tddfpt2_types
Definition qs_tddfpt2_types.F:8
soc_pseudopotential_methods
Definition soc_pseudopotential_methods.F:8
soc_pseudopotential_methods::v_soc_xyz_from_pseudopotential
subroutine, public v_soc_xyz_from_pseudopotential(qs_env, mat_v_soc_xyz)
V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,...
Definition soc_pseudopotential_methods.F:70
spherical_harmonics
Calculate spherical harmonics.
Definition spherical_harmonics.F:27
spherical_harmonics::clebsch_gordon_init
subroutine, public clebsch_gordon_init(l)
...
Definition spherical_harmonics.F:238
spherical_harmonics::clebsch_gordon_deallocate
subroutine, public clebsch_gordon_deallocate()
...
Definition spherical_harmonics.F:220
xas_tdp_atom
This module deals with all the integrals done on local atomic grids in xas_tdp. This is mostly used t...
Definition xas_tdp_atom.F:15
xas_tdp_atom::compute_sphi_so
subroutine, public compute_sphi_so(ikind, basis_type, sphi_so, qs_env)
Computes the contraction coefficients to go from spherical orbitals to sgf for a given atomic kind an...
Definition xas_tdp_atom.F:453
xas_tdp_atom::integrate_soc_atoms
subroutine, public integrate_soc_atoms(matrix_soc, xas_atom_env, qs_env, soc_atom_env)
Computes the SOC matrix elements with respect to the ORB basis set for each atomic kind and put them ...
Definition xas_tdp_atom.F:3372
xas_tdp_atom::truncate_radial_grid
subroutine, public truncate_radial_grid(grid_atom, max_radius)
Reduces the radial extension of an atomic grid such that it only covers a given radius.
Definition xas_tdp_atom.F:406
xas_tdp_utils
Utilities for X-ray absorption spectroscopy using TDDFPT.
Definition xas_tdp_utils.F:13
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:2884
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
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_control_types::dft_control_type
Definition cp_control_types.F:570
cp_control_types::qs_control_type
Definition cp_control_types.F:358
cp_control_types::tddfpt2_control_type
Definition cp_control_types.F:469
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:180
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_grid_atom::grid_atom_type
Definition qs_grid_atom.F:44
qs_harmonics_atom::harmonics_atom_type
Definition qs_harmonics_atom.F:32
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_tddfpt2_soc_types::soc_atom_env_type
Definition qs_tddfpt2_soc_types.F:72
qs_tddfpt2_soc_types::soc_env_type
Definition qs_tddfpt2_soc_types.F:49
qs_tddfpt2_types::tddfpt_ground_state_mos
Ground state molecular orbitals.
Definition qs_tddfpt2_types.F:86