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