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