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