Line data Source code
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 : !--------------------------------------------------------------------------------------------------!
7 : MODULE 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 gaunt_coefficients, ONLY: real_gaunt_coefficients
44 : USE input_constants, ONLY: tddfpt_dipole_length
45 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
46 : section_vals_type,&
47 : section_vals_val_get
48 : USE kinds, ONLY: default_string_length,&
49 : dp
50 : USE lebedev, ONLY: deallocate_lebedev_grids,&
51 : get_number_of_lebedev_grid,&
52 : init_lebedev_grids,&
53 : lebedev_grid
54 : USE mathlib, ONLY: get_diag
55 : USE message_passing, ONLY: mp_para_env_type
56 : USE orbital_pointers, ONLY: nsoset
57 : USE parallel_gemm_api, ONLY: parallel_gemm
58 : USE particle_types, ONLY: particle_type
59 : USE physcon, ONLY: evolt,&
60 : wavenumbers
61 : USE qs_environment_types, ONLY: get_qs_env,&
62 : qs_environment_type
63 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
64 : create_grid_atom,&
65 : grid_atom_type
66 : USE qs_harmonics_atom, ONLY: allocate_harmonics_atom,&
67 : create_harmonics_atom,&
68 : get_maxl_CG,&
69 : harmonics_atom_type
70 : USE qs_kind_types, ONLY: get_qs_kind,&
71 : get_qs_kind_set,&
72 : qs_kind_type
73 : USE qs_mo_types, ONLY: mo_set_type
74 : USE qs_tddfpt2_soc_types, ONLY: soc_atom_create,&
75 : soc_atom_env_type,&
76 : soc_atom_release,&
77 : soc_env_create,&
78 : soc_env_release,&
79 : soc_env_type
80 : USE qs_tddfpt2_soc_utils, ONLY: dip_vel_op,&
81 : resort_evects,&
82 : soc_contract_evect,&
83 : soc_dipole_operator
84 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos
85 : USE soc_pseudopotential_methods, ONLY: V_SOC_xyz_from_pseudopotential
86 : USE xas_tdp_atom, ONLY: compute_sphi_so,&
87 : integrate_soc_atoms,&
88 : truncate_radial_grid
89 : USE xas_tdp_utils, ONLY: rcs_amew_soc_elements
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 :
119 : CONTAINS
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 10 : 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 10 : CALL timeset(routineN, handle)
151 :
152 10 : logger => cp_get_default_logger()
153 10 : IF (logger%para_env%is_source()) THEN
154 5 : log_unit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
155 : ELSE
156 5 : log_unit = -1
157 : END IF
158 :
159 10 : IF (log_unit > 0) THEN
160 5 : WRITE (log_unit, "(1X,A)") "", &
161 5 : "-------------------------------------------------------------------------------", &
162 5 : "- START SOC CALCULATIONS -", &
163 10 : "-------------------------------------------------------------------------------"
164 : END IF
165 :
166 10 : NULLIFY (dft_control)
167 10 : CALL get_qs_env(qs_env, dft_control=dft_control)
168 10 : do_os = dft_control%uks .OR. dft_control%roks
169 0 : IF (do_os) CPABORT("SOC only implemented for closed shell.")
170 :
171 10 : IF (log_unit > 0) THEN
172 5 : WRITE (log_unit, '(A)') " Starting from TDDFPT Excited States:"
173 5 : WRITE (log_unit, '(A)') " STATE SINGLET/eV TRIPLET/eV"
174 19 : DO istate = 1, SIZE(evals_a)
175 19 : 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 10 : IF (log_unit > 0) WRITE (log_unit, '(A)') " Starting restricted closed shell:"
180 10 : CALL tddfpt_soc_rcs(qs_env, evals_a, evals_b, evects_a, evects_b, log_unit, gs_mos)
181 :
182 10 : IF (log_unit > 0) THEN
183 5 : WRITE (log_unit, '(A,/,A)') "SOC Calculation terminated", &
184 10 : "Returning to TDDFPT for Force calculation and deallocations"
185 : END IF
186 :
187 10 : CALL timestop(handle)
188 :
189 10 : 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 10 : 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 10 : INTEGER, ALLOCATABLE, DIMENSION(:) :: evects_sort
220 10 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, col_dist, row_blk_size, &
221 10 : row_dist, row_dist_new
222 10 : 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 10 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag, tmp_evals
227 10 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: img_tmp, real_tmp
228 10 : 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 10 : 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 10 : TYPE(soc_env_type), TARGET :: soc_env
250 :
251 10 : CALL timeset(routineN, handle)
252 :
253 10 : NULLIFY (logger, pgrid, row_dist, row_dist_new)
254 10 : NULLIFY (col_dist, col_blk_size, row_blk_size, matrix_s)
255 10 : NULLIFY (gstp_struct, tddfpt_print_section, soc_print_section, soc_section)
256 :
257 10 : logger => cp_get_default_logger()
258 10 : tddfpt_print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT")
259 10 : soc_print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT%SOC_PRINT")
260 10 : soc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%SOC")
261 :
262 10 : CALL section_vals_val_get(soc_section, "EPS_FILTER", r_val=eps_filter)
263 :
264 10 : nsg = SIZE(evals_sing) ! Number of excited singlet states
265 10 : ntp = SIZE(evals_trip) ! Number of excited triplet states
266 10 : nex = nsg ! Number of excited states of each multiplicity
267 10 : 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 10 : evects_sing, evects_trip, dipole_form, gs_mos)
272 :
273 10 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
274 10 : nactive = gs_mos(1)%nmo_active
275 10 : 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 10 : CALL integrate_soc_atoms(soc_env%orb_soc, qs_env=qs_env, soc_atom_env=soc_atom_env)
279 10 : CALL soc_atom_release(soc_atom_env)
280 :
281 : !! Point at H^SOC and MOs for better readablity
282 10 : NULLIFY (tp_coeffs, orb_soc_x, orb_soc_y, orb_soc_z)
283 10 : tp_coeffs => soc_env%b_coeff
284 10 : soc_env%evals_a => evals_sing
285 10 : soc_env%evals_b => evals_trip
286 10 : orb_soc_x => soc_env%orb_soc(1)%matrix
287 10 : orb_soc_y => soc_env%orb_soc(2)%matrix
288 10 : orb_soc_z => soc_env%orb_soc(3)%matrix
289 :
290 : !! Create a matrix-structure, which links all states in this calculation
291 10 : NULLIFY (full_struct)
292 10 : CALL cp_fm_struct_create(full_struct, context=blacs_env, para_env=para_env, nrow_global=ntot, ncol_global=ntot)
293 10 : CALL cp_fm_create(real_fm, full_struct)
294 10 : CALL cp_fm_create(img_fm, full_struct)
295 10 : CALL cp_fm_set_all(real_fm, 0.0_dp)
296 10 : CALL cp_fm_set_all(img_fm, 0.0_dp)
297 :
298 : ! Put the excitation energies on the diagonal of the real matrix
299 38 : DO isg = 1, nsg
300 38 : CALL cp_fm_set_element(real_fm, 1 + isg, 1 + isg, evals_sing(isg))
301 : END DO
302 38 : DO itp = 1, ntp
303 : ! first T^-1, then T^0, then T^+1
304 28 : CALL cp_fm_set_element(real_fm, 1 + itp + nsg, 1 + itp + nsg, evals_trip(itp))
305 28 : CALL cp_fm_set_element(real_fm, 1 + itp + ntp + nsg, 1 + itp + ntp + nsg, evals_trip(itp))
306 38 : 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 10 : 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 10 : npcols=npcols, nprows=nprows)
313 40 : ALLOCATE (col_dist(nex), row_dist_new(nex)) ! Split for each excitation
314 38 : DO iex = 1, nex
315 28 : col_dist(iex) = MODULO(npcols - iex, npcols)
316 38 : row_dist_new(iex) = MODULO(nprows - iex, nprows)
317 : END DO
318 10 : ALLOCATE (coeffs_dist, prod_dist)
319 : CALL dbcsr_distribution_new(coeffs_dist, group=group, pgrid=pgrid, row_dist=row_dist, &
320 10 : col_dist=col_dist)
321 : CALL dbcsr_distribution_new(prod_dist, group=group, pgrid=pgrid, row_dist=row_dist_new, &
322 10 : col_dist=col_dist)
323 :
324 : !! Create the matrices
325 20 : ALLOCATE (col_blk_size(nex))
326 38 : col_blk_size = nactive
327 :
328 10 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
329 10 : 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 10 : 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 10 : 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 10 : 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 10 : 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 10 : 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 10 : dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
343 :
344 38 : col_blk_size = 1
345 : CALL dbcsr_create(matrix=dbcsr_tmp, name="TMP", matrix_type=dbcsr_type_no_symmetry, &
346 10 : dist=prod_dist, row_blk_size=col_blk_size, col_blk_size=col_blk_size)
347 10 : 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 10 : dbcsr_soc_package%dbcsr_sg => dbcsr_sg
358 10 : dbcsr_soc_package%dbcsr_tp => dbcsr_tp
359 10 : dbcsr_soc_package%dbcsr_work => dbcsr_work
360 10 : dbcsr_soc_package%dbcsr_ovlp => dbcsr_ovlp
361 10 : dbcsr_soc_package%dbcsr_prod => dbcsr_prod
362 10 : dbcsr_soc_package%dbcsr_tmp => dbcsr_tmp
363 :
364 : !Filling the coeffs matrices by copying from the stored fms
365 10 : CALL copy_fm_to_dbcsr(soc_env%a_coeff, dbcsr_sg)
366 10 : CALL copy_fm_to_dbcsr(soc_env%b_coeff, dbcsr_tp)
367 :
368 : !Create the work and helper fms
369 10 : gs_coeffs => gs_mos(1)%mos_active
370 10 : CALL cp_fm_get_info(gs_coeffs, matrix_struct=vec_struct)
371 10 : CALL cp_fm_create(vec_soc_x, vec_struct)
372 10 : CALL cp_fm_create(vec_soc_y, vec_struct)
373 10 : 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 10 : nrow_global=nactive, ncol_global=nactive)
377 10 : 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 10 : nrow_global=nex, ncol_global=nex)
381 10 : CALL cp_fm_create(work_fm, work_struct)
382 10 : 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 10 : sqrt2 = SQRT(2.0_dp)
394 :
395 : !! Precompute the <phi_i^0|H^SOC|phi_j^0> matrix elements
396 30 : ALLOCATE (diag(nactive))
397 80 : 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 10 : 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 10 : 0.0_dp, prod_fm)
405 10 : CALL cp_fm_get_submatrix(prod_fm, mo_soc_x)
406 :
407 10 : 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 10 : 0.0_dp, prod_fm)
412 10 : CALL cp_fm_get_submatrix(prod_fm, mo_soc_y)
413 :
414 10 : 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 10 : 0.0_dp, prod_fm)
419 10 : 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 10 : nrow_global=ntp*nactive, ncol_global=nactive)
432 10 : CALL cp_fm_create(gstp_fm, gstp_struct)
433 30 : 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 10 : 0.0_dp, gstp_fm)
442 :
443 : !! Seperate them into the different states again (nactive x nactive)
444 38 : 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 28 : n_rows=nactive, n_cols=nactive)
448 28 : diag(:) = get_diag(gstp_block)
449 200 : soc_gst = SUM(diag)
450 : !! <0|H_x|T^-1>
451 28 : CALL cp_fm_set_element(img_fm, 1, 1 + nsg + itp, -1_dp*soc_gst)
452 : !! <0|H_x|T^+1>
453 38 : 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 10 : 0.0_dp, gstp_fm)
467 :
468 38 : DO itp = 1, ntp
469 : CALL cp_fm_get_submatrix(fm=gstp_fm, target_m=gstp_block, start_row=(itp - 1)*nactive + 1, &
470 28 : start_col=1, n_rows=nactive, n_cols=nactive)
471 28 : diag(:) = get_diag(gstp_block)
472 200 : soc_gst = SUM(diag)
473 : ! <0|H_y|T^-1>
474 28 : CALL cp_fm_set_element(real_fm, 1, 1 + nsg + itp, -1.0_dp*soc_gst)
475 : ! <0|H_y|T^+1>
476 38 : 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 10 : 0.0_dp, gstp_fm)
490 :
491 38 : 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 28 : n_cols=nactive)
498 28 : diag(:) = get_diag(gstp_block)
499 200 : soc_gst = sqrt2*SUM(diag)
500 38 : 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 10 : CALL cp_fm_release(prod_fm)
514 10 : CALL cp_fm_release(vec_soc_x)
515 10 : CALL cp_fm_release(vec_soc_y)
516 10 : CALL cp_fm_release(vec_soc_z)
517 10 : CALL cp_fm_release(gstp_fm)
518 10 : CALL cp_fm_struct_release(gstp_struct)
519 10 : CALL cp_fm_struct_release(prod_struct)
520 10 : 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 10 : 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 10 : 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 10 : 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 10 : 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 10 : pref_overall=-0.5_dp*sqrt2)
561 :
562 : !<S|H_x|T^-1>
563 : !! Convert to fm for transfer to img_fm
564 10 : 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 10 : 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 10 : 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 10 : 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 10 : 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 10 : 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 10 : pref_overall=-0.5_dp*sqrt2)
617 :
618 : !<S|H_y|T^-1>
619 10 : 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 10 : 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 10 : 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 10 : 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 10 : 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 10 : pref_overall=1.0_dp)
669 :
670 : !<S|H_z|T^0>
671 10 : 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 10 : 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 10 : 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 10 : 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 10 : 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 10 : 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 10 : pref_overall=-0.5_dp*sqrt2)
720 :
721 : !<T^0|H_x|T^+1>
722 10 : 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 10 : 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 10 : CALL cp_fm_transpose(tmp_fm, work_fm)
741 10 : 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 10 : 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 10 : 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 10 : 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 10 : pref_overall=0.5_dp*sqrt2)
774 :
775 : !<T^0|H_y|T^+1>
776 10 : 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 10 : 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 10 : CALL cp_fm_transpose(tmp_fm, work_fm)
794 10 : 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 10 : 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 10 : 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 10 : 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 10 : pref_overall=1.0_dp)
826 :
827 : !<T^+1|H_z|T^+1>
828 10 : 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 10 : 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 10 : 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 10 : 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 10 : CALL cp_fm_struct_release(work_struct)
865 10 : CALL cp_fm_release(work_fm)
866 10 : CALL cp_fm_release(tmp_fm)
867 10 : DEALLOCATE (diag, mo_soc_x, mo_soc_y, mo_soc_z)
868 :
869 10 : IF (log_unit > 0) THEN
870 5 : WRITE (log_unit, '(A)') "Diagonlzing SOC-Matrix"
871 : END IF
872 :
873 : !Set-up the complex hermitian perturbation matrix
874 10 : CALL cp_cfm_create(hami_cfm, full_struct)
875 10 : CALL cp_fm_to_cfm(real_fm, img_fm, hami_cfm)
876 :
877 : !!Optinal Output: SOC-Matrix
878 10 : 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 5 : file_status="UNKNOWN")
886 : ELSE
887 5 : log_unit = -1
888 : END IF
889 :
890 : !! Get the requested energy unit and optional outputs
891 10 : CALL section_vals_val_get(soc_print_section, "UNIT_eV", l_val=print_ev)
892 10 : CALL section_vals_val_get(soc_print_section, "UNIT_wn", l_val=print_wn)
893 10 : CALL section_vals_val_get(soc_print_section, "SPLITTING", l_val=print_splitting)
894 10 : 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 10 : IF (print_wn) THEN
898 2 : print_ev = .FALSE.
899 2 : unit_scale = wavenumbers
900 8 : ELSE IF (print_ev) THEN
901 : unit_scale = evolt
902 : ELSE
903 0 : CPWARN("No output unit was selected, printout will be in Hartree!")
904 0 : 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 10 : IF (print_some) THEN
910 2 : IF (log_unit > 0) THEN
911 1 : WRITE (log_unit, '(A)') "____________________________________________________"
912 1 : WRITE (log_unit, '(A)') " FULL SOC-Matrix "
913 1 : IF (print_ev) THEN
914 0 : WRITE (log_unit, '(A)') "STATE STATE REAL-PART[eV] IMG-PART[eV] "
915 1 : ELSE IF (print_wn) THEN
916 1 : WRITE (log_unit, '(A)') "STATE STATE REAL-PART[cm-1] IMG-PART[cm-1]"
917 : ELSE
918 0 : WRITE (log_unit, '(A)') "STATE STATE REAL-PART[Hartree] IMG-PART[Hartree]"
919 : END IF
920 1 : WRITE (log_unit, '(A)') "____________________________________________________"
921 : END IF
922 12 : ALLOCATE (real_tmp(ntot, ntot), img_tmp(ntot, ntot))
923 2 : real_tmp = 0_dp
924 2 : img_tmp = 0_dp
925 :
926 2 : CALL cp_fm_get_submatrix(real_fm, real_tmp, 1, 1, ntot, ntot)
927 2 : CALL cp_fm_get_submatrix(img_fm, img_tmp, 1, 1, ntot, ntot)
928 :
929 2 : IF (log_unit > 0) THEN
930 10 : DO jj = 1, ntot
931 91 : DO ii = 1, ntot
932 81 : WRITE (unit=log_unit, fmt="(I3,5X,I3,5X,f16.8,5X,f16.8)") ii, jj, &
933 81 : real_tmp(ii, jj)*unit_scale, &
934 171 : img_tmp(ii, jj)*unit_scale
935 : END DO !! jj
936 : END DO !! ii
937 1 : WRITE (log_unit, '(A)') " END SOC-MATRIX "
938 1 : WRITE (log_unit, '(A)') "____________________________________________________"
939 1 : WRITE (log_unit, '(A)') "____________________________________________________"
940 : END IF !(log_unit)
941 2 : DEALLOCATE (real_tmp, img_tmp)
942 : END IF !print_some
943 :
944 10 : CALL cp_fm_release(real_fm)
945 10 : CALL cp_fm_release(img_fm)
946 :
947 : ! Diagonalize the Hamiltonian
948 30 : ALLOCATE (tmp_evals(ntot))
949 10 : CALL cp_cfm_create(evecs_cfm, full_struct)
950 10 : 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 30 : ALLOCATE (soc_env%soc_evals(ntot - 1))
954 122 : soc_env%soc_evals(:) = tmp_evals(2:ntot) - tmp_evals(1)
955 :
956 : !! We may be interested in the ground state stabilisation energy
957 10 : IF (logger%para_env%is_source()) THEN
958 5 : log_unit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
959 : ELSE
960 5 : log_unit = -1
961 : END IF
962 :
963 10 : IF (log_unit > 0) THEN
964 5 : 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 10 : soc_env%evals_a => evals_sing
972 10 : 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 10 : 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 10 : IF (logger%para_env%is_source()) THEN
980 5 : log_unit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
981 : ELSE
982 5 : log_unit = -1
983 : END IF
984 :
985 10 : IF (log_unit > 0) THEN
986 5 : WRITE (log_unit, '(A)') "------------------------------------------------------------------------------"
987 5 : WRITE (log_unit, '(A)') "- SOC CORRECTED SPECTRUM -"
988 5 : WRITE (log_unit, '(A)') "------------------------------------------------------------------------------"
989 5 : IF (print_ev) THEN
990 4 : WRITE (log_unit, '(A)') "- STATE SOC-corrected exc. energies [eV] Oscillator strengths [a.u.] -"
991 1 : ELSE IF (print_wn) THEN
992 1 : WRITE (log_unit, '(A)') "- STATE SOC-corrected exc. energies [cm-1] Oscillator strengths [a.u.]-"
993 : ELSE
994 0 : WRITE (log_unit, '(A)') "- STATE SOC-corrected exc.energies [Hartree] Oscillator strengths [a.u.]-"
995 : END IF
996 5 : WRITE (log_unit, '(A)') "------------------------------------------------------------------------------"
997 61 : DO istate = 1, SIZE(soc_env%soc_evals)
998 56 : WRITE (log_unit, '(6X,I5,11X,2F16.5)') istate, soc_env%soc_evals(istate)*unit_scale, &
999 117 : soc_env%soc_osc(istate)
1000 : END DO
1001 : !! This is used for regtesting
1002 61 : 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 10 : IF (print_splitting .OR. debug_this_module) THEN
1009 8 : CALL resort_evects(evecs_cfm, evects_sort)
1010 8 : IF (log_unit > 0) THEN
1011 4 : WRITE (log_unit, '(A)') "-----------------------------------------------------------------------------"
1012 4 : WRITE (log_unit, '(A)') "- Splittings -"
1013 4 : WRITE (log_unit, '(A)') "-----------------------------------------------------------------------------"
1014 4 : IF (print_ev) THEN
1015 3 : WRITE (log_unit, '(A)') "- STATE exc.energies[eV] splittings[eV] -"
1016 1 : ELSE IF (print_wn) THEN
1017 1 : WRITE (log_unit, '(A)') "- STATE exc.energies[cm-1] splittings[cm-1] -"
1018 : ELSE
1019 0 : WRITE (log_unit, '(A)') "- STATE exc.energies[Hartree] splittings[Hartree] -"
1020 : END IF
1021 4 : WRITE (log_unit, '(A)') "-----------------------------------------------------------------------------"
1022 4 : mult = " "
1023 48 : DO iex = 1, SIZE(soc_env%soc_evals)
1024 44 : IF (evects_sort(iex + 1) > nsg + ntp*2 + 1) THEN
1025 11 : mult = "T+1"
1026 33 : ELSE IF (evects_sort(iex + 1) > nsg + ntp + 1) THEN
1027 11 : mult = "T0"
1028 22 : ELSE IF (evects_sort(iex + 1) > nsg + 1) THEN
1029 11 : mult = "T-1"
1030 : ELSE
1031 11 : mult = "S "
1032 : END IF
1033 48 : IF (mult == "T-1") THEN
1034 11 : WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
1035 11 : evects_sort(iex + 1) - (1 + nsg), soc_env%soc_evals(iex)*unit_scale, &
1036 22 : (soc_env%soc_evals(iex) - evals_trip(evects_sort(iex + 1) - (1 + nsg)))*unit_scale
1037 33 : ELSE IF (mult == "T0") THEN
1038 11 : WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
1039 11 : evects_sort(iex + 1) - (1 + nsg + ntp), soc_env%soc_evals(iex)*unit_scale, &
1040 22 : (soc_env%soc_evals(iex) - evals_trip(evects_sort(iex + 1) - (1 + nsg + ntp)))*unit_scale
1041 22 : ELSE IF (mult == "T+1") THEN
1042 11 : WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
1043 11 : evects_sort(iex + 1) - (1 + nsg + ntp*2), soc_env%soc_evals(iex)*unit_scale, &
1044 22 : (soc_env%soc_evals(iex) - evals_trip(evects_sort(iex + 1) - (1 + nsg + ntp*2)))*unit_scale
1045 11 : ELSE IF (mult == "S ") THEN
1046 11 : WRITE (log_unit, "(5X,A,I5,6X,F18.10,6X,F18.10)") mult, &
1047 11 : evects_sort(iex + 1), soc_env%soc_evals(iex)*unit_scale, &
1048 22 : (soc_env%soc_evals(iex) - evals_sing(evects_sort(iex + 1) - 1))*unit_scale
1049 : END IF
1050 : END DO
1051 4 : WRITE (log_unit, '(A)') "----------------------------------------------------------------------------"
1052 4 : DEALLOCATE (evects_sort)
1053 : END IF
1054 : END IF
1055 :
1056 : !! clean up
1057 10 : CALL soc_env_release(soc_env)
1058 10 : CALL cp_fm_struct_release(full_struct)
1059 10 : CALL cp_cfm_release(hami_cfm)
1060 10 : CALL cp_cfm_release(evecs_cfm)
1061 10 : CALL dbcsr_distribution_release(coeffs_dist)
1062 10 : CALL dbcsr_distribution_release(prod_dist)
1063 10 : CALL dbcsr_release(dbcsr_sg)
1064 10 : CALL dbcsr_release(dbcsr_tp)
1065 10 : CALL dbcsr_release(dbcsr_prod)
1066 10 : CALL dbcsr_release(dbcsr_ovlp)
1067 10 : CALL dbcsr_release(dbcsr_tmp)
1068 10 : 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 10 : DEALLOCATE (dbcsr_work, dbcsr_prod, dbcsr_ovlp, dbcsr_tmp)
1075 10 : DEALLOCATE (coeffs_dist, prod_dist, col_dist, col_blk_size, row_dist_new)
1076 10 : DEALLOCATE (dbcsr_sg, dbcsr_tp, tmp_evals)
1077 :
1078 10 : CALL timestop(handle)
1079 :
1080 140 : 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 10 : SUBROUTINE inititialize_soc(qs_env, soc_atom_env, soc_env, &
1093 10 : evects_a, evects_b, &
1094 10 : 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 10 : 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 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1118 : TYPE(dft_control_type), POINTER :: dft_control
1119 10 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1120 : TYPE(mp_para_env_type), POINTER :: para_env
1121 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1122 10 : 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 10 : CALL timeset(routineN, handle)
1127 :
1128 10 : NULLIFY (qs_kind_set, particle_set, blacs_env, para_env, &
1129 10 : a_coeff, b_coeff, tddfpt_control, soc_section)
1130 :
1131 : !! Open_shell is not implemented yet
1132 10 : do_os = .FALSE.
1133 :
1134 : !! soc_env will pass most of the larger matrices to the different modules
1135 10 : 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 10 : dft_control=dft_control)
1141 :
1142 : !! The Dipole form shall be consistent with the tddfpt-module!
1143 10 : tddfpt_control => dft_control%tddfpt2_control
1144 10 : dipole_form = tddfpt_control%dipole_form
1145 :
1146 10 : 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 10 : CALL section_vals_val_get(soc_section, "GRID", n_rep_val=nrep)
1150 20 : ALLOCATE (grid_info(nrep, 3))
1151 10 : DO irep = 1, nrep
1152 : CALL section_vals_val_get(soc_section, &
1153 : "GRID", &
1154 : i_rep_val=irep, &
1155 0 : c_vals=k_list)
1156 10 : grid_info(irep, :) = k_list
1157 : END DO
1158 :
1159 10 : CALL soc_dipole_operator(soc_env, tddfpt_control, qs_env, gs_mos)
1160 :
1161 10 : nspins = SIZE(evects_a, 1)
1162 20 : DO ispin = 1, nspins
1163 20 : 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 10 : nstates = SIZE(evects_a, 2)
1169 :
1170 10 : a_coeff => soc_env%a_coeff
1171 10 : b_coeff => soc_env%b_coeff
1172 :
1173 10 : NULLIFY (fm_struct)
1174 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
1175 10 : nrow_global=nao, ncol_global=nspins*nactive*nstates)
1176 10 : CALL cp_fm_create(a_coeff, fm_struct)
1177 10 : CALL cp_fm_create(b_coeff, fm_struct)
1178 10 : CALL cp_fm_struct_release(fm_struct)
1179 :
1180 10 : CALL soc_contract_evect(evects_a, a_coeff)
1181 10 : CALL soc_contract_evect(evects_b, b_coeff)
1182 :
1183 40 : ALLOCATE (soc_env%orb_soc(3))
1184 40 : DO i = 1, 3
1185 40 : 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 10 : CALL soc_atom_create(soc_atom_env)
1191 : IF (do_os) THEN
1192 : soc_atom_env%nspins = 2
1193 : ELSE
1194 10 : soc_atom_env%nspins = 1
1195 : END IF
1196 :
1197 10 : nkind = SIZE(qs_kind_set)
1198 40 : ALLOCATE (soc_atom_env%grid_atom_set(nkind))
1199 30 : ALLOCATE (soc_atom_env%harmonics_atom_set(nkind))
1200 40 : ALLOCATE (soc_atom_env%orb_sphi_so(nkind))
1201 :
1202 10 : CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present)
1203 10 : IF (.NOT. all_potential_present) THEN
1204 2 : 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 10 : CALL init_atom_grid(soc_atom_env, grid_info, qs_env)
1210 :
1211 20 : DO ikind = 1, nkind
1212 10 : NULLIFY (soc_atom_env%orb_sphi_so(ikind)%array)
1213 20 : CALL compute_sphi_so(ikind, "ORB", soc_atom_env%orb_sphi_so(ikind)%array, qs_env)
1214 : END DO !ikind
1215 :
1216 10 : DEALLOCATE (grid_info)
1217 :
1218 10 : CALL timestop(handle)
1219 :
1220 50 : 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 20 : 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 10 : 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 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1250 :
1251 10 : CALL timeset(routineN, handle)
1252 :
1253 10 : NULLIFY (my_CG, qs_kind_set, dft_control, qs_control)
1254 10 : 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 10 : dft_control=dft_control)
1260 :
1261 10 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="ORB")
1262 10 : qs_control => dft_control%qs_control
1263 :
1264 : !!To be sure:: Is the basis grid present?
1265 10 : IF (maxlgto < 0) THEN
1266 0 : CPABORT("tmp_basis is not Present!")
1267 : END IF
1268 :
1269 : !maximum expansion
1270 10 : llmax = 2*maxlgto
1271 10 : max_s_harm = nsoset(llmax)
1272 10 : max_s_set = nsoset(maxlgto)
1273 : ! Compute the Gaunt coefficients
1274 10 : CALL real_gaunt_coefficients(my_CG, maxlgto, maxlgto, llmax)
1275 :
1276 : ! Create the Lebedev grids and compute the spherical harmonics
1277 10 : CALL init_lebedev_grids()
1278 10 : quadrature = qs_control%gapw_control%quadrature
1279 :
1280 20 : DO ikind = 1, SIZE(soc_atom_env%grid_atom_set)
1281 :
1282 : ! Allocate the grid and the harmonics for this kind
1283 10 : NULLIFY (soc_atom_env%grid_atom_set(ikind)%grid_atom)
1284 10 : NULLIFY (soc_atom_env%harmonics_atom_set(ikind)%harmonics_atom)
1285 10 : CALL allocate_grid_atom(soc_atom_env%grid_atom_set(ikind)%grid_atom)
1286 : CALL allocate_harmonics_atom( &
1287 10 : soc_atom_env%harmonics_atom_set(ikind)%harmonics_atom)
1288 :
1289 : NULLIFY (grid_atom, harmonics)
1290 10 : grid_atom => soc_atom_env%grid_atom_set(ikind)%grid_atom
1291 10 : 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 10 : name=kind_name)
1298 :
1299 : ! take the grid dimension given as input, if none, take the GAPW ones above
1300 10 : DO igrid = 1, SIZE(grid_info, 1)
1301 10 : IF (grid_info(igrid, 1) == kind_name) THEN
1302 0 : READ (grid_info(igrid, 2), *, iostat=stat) na
1303 0 : IF (stat /= 0) CPABORT("The 'na' value for the GRID keyword must be an integer")
1304 0 : READ (grid_info(igrid, 3), *, iostat=stat) nr
1305 0 : 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 :
1310 10 : ll = get_number_of_lebedev_grid(n=na)
1311 10 : na = lebedev_grid(ll)%n
1312 10 : la = lebedev_grid(ll)%l
1313 10 : grid_atom%ng_sphere = na
1314 10 : 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 10 : basis_type="ORB")
1320 : CALL get_gto_basis_set(gto_basis_set=tmp_basis, &
1321 : maxl=maxl, &
1322 10 : kind_radius=kind_radius)
1323 :
1324 10 : CALL create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
1325 10 : CALL truncate_radial_grid(grid_atom, kind_radius)
1326 :
1327 10 : 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 10 : grid_atom%pol)
1335 40 : CALL get_maxl_CG(harmonics, tmp_basis, llmax, max_s_harm)
1336 : END DO !ikind
1337 :
1338 10 : CALL deallocate_lebedev_grids()
1339 10 : DEALLOCATE (my_CG)
1340 :
1341 10 : CALL timestop(handle)
1342 :
1343 10 : 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 10 : SUBROUTINE soc_dipol(qs_env, dbcsr_soc_package, soc_env, &
1358 : soc_evecs_cfm, eps_filter, dipole_form, gs_mos, &
1359 10 : 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 10 : 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 10 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: amew_dip
1379 : TYPE(mp_para_env_type), POINTER :: para_env
1380 :
1381 10 : CALL timeset(routineN, handle)
1382 :
1383 10 : NULLIFY (para_env, blacs_env, dip_struct, full_struct, osc_str)
1384 10 : NULLIFY (soc_evals)
1385 :
1386 10 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1387 :
1388 10 : soc_evals => soc_env%soc_evals
1389 10 : nosc = SIZE(soc_evals)
1390 10 : ntot = nosc + 1
1391 30 : ALLOCATE (soc_env%soc_osc(nosc))
1392 10 : osc_str => soc_env%soc_osc
1393 122 : 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 10 : nrow_global=ntot, ncol_global=1)
1399 :
1400 10 : CALL cp_cfm_get_info(soc_evecs_cfm, matrix_struct=full_struct)
1401 10 : CALL cp_cfm_create(dip_cfm, dip_struct, set_zero=.TRUE.)
1402 10 : CALL cp_cfm_create(work1_cfm, full_struct, set_zero=.TRUE.)
1403 10 : CALL cp_cfm_create(work2_cfm, full_struct, set_zero=.TRUE.)
1404 :
1405 30 : 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 10 : evects_sing, evects_trip)
1410 :
1411 40 : DO i = 1, 3 !cartesian coord x, y, z
1412 :
1413 : !Convert the real dipole into the cfm format for calculations
1414 30 : 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 30 : 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 30 : dip_cfm)
1431 :
1432 30 : CALL cp_cfm_get_submatrix(dip_cfm, transdip)
1433 :
1434 : !transition dipoles are real numbers
1435 376 : 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 10 : IF (dipole_form == tddfpt_dipole_length) THEN
1441 208 : osc_str(:) = 2.0_dp/3.0_dp*soc_evals(:)*osc_str(:)
1442 : ELSE
1443 36 : osc_str(:) = 2.0_dp/3.0_dp*soc_evals(:)*osc_str(:)
1444 : END IF
1445 :
1446 : !clean-up
1447 10 : CALL cp_fm_struct_release(dip_struct)
1448 10 : CALL cp_cfm_release(work1_cfm)
1449 10 : CALL cp_cfm_release(work2_cfm)
1450 10 : CALL cp_cfm_release(dip_cfm)
1451 40 : DO i = 1, 3
1452 40 : CALL cp_fm_release(amew_dip(i))
1453 : END DO
1454 10 : DEALLOCATE (amew_dip, transdip)
1455 :
1456 10 : CALL timestop(handle)
1457 :
1458 10 : 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 10 : SUBROUTINE get_rcs_amew_op(amew_op, soc_env, dbcsr_soc_package, eps_filter, qs_env, gs_mos, &
1472 10 : 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 10 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: diag, gs_diag, gsgs_op
1491 10 : 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 10 : 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 10 : CALL timeset(routineN, handle)
1504 :
1505 10 : NULLIFY (gs_coeffs)
1506 10 : NULLIFY (matrix_s)
1507 10 : NULLIFY (para_env, blacs_env)
1508 10 : NULLIFY (full_struct, gsgs_struct, std_struct, tmp_struct, sggs_struct)
1509 10 : NULLIFY (ao_op_i, ao_op)
1510 10 : 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 10 : IF (ASSOCIATED(soc_env%dipmat)) THEN
1515 8 : ao_op => soc_env%dipmat
1516 8 : vel_rep = .FALSE.
1517 : ELSE
1518 2 : ao_op => soc_env%dipmat_ao
1519 2 : vel_rep = .TRUE.
1520 : END IF
1521 10 : nsg = SIZE(soc_env%evals_a)
1522 10 : ntp = nsg; nex = nsg !all the same by construction, keep them separate for clarity
1523 10 : ntot = 1 + nsg + 3*ntp
1524 :
1525 10 : CALL get_qs_env(qs_env, matrix_s=matrix_s, para_env=para_env, blacs_env=blacs_env)
1526 10 : sqrt2 = SQRT(2.0_dp)
1527 10 : dim_op = SIZE(ao_op)
1528 :
1529 10 : dbcsr_sg => dbcsr_soc_package%dbcsr_sg
1530 10 : dbcsr_tp => dbcsr_soc_package%dbcsr_tp
1531 10 : dbcsr_work => dbcsr_soc_package%dbcsr_work
1532 10 : dbcsr_prod => dbcsr_soc_package%dbcsr_prod
1533 10 : dbcsr_ovlp => dbcsr_soc_package%dbcsr_ovlp
1534 10 : 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 10 : nrow_global=ntot, ncol_global=ntot)
1540 60 : ALLOCATE (amew_op(dim_op))
1541 40 : DO i = 1, dim_op
1542 40 : 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 10 : nactive = gs_mos(1)%nmo_active
1547 10 : gs_coeffs => gs_mos(1)%mos_active
1548 10 : 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 10 : nrow_global=nactive, ncol_global=nactive)
1553 :
1554 10 : CALL cp_fm_create(gs_fm, gsgs_struct) !nactive nactive
1555 10 : CALL cp_fm_create(work_fm, std_struct) !nao nactive
1556 :
1557 30 : ALLOCATE (gsgs_op(dim_op))
1558 30 : ALLOCATE (gs_diag(nactive))
1559 :
1560 40 : DO i = 1, dim_op
1561 :
1562 30 : ao_op_i => ao_op(i)%matrix
1563 30 : 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 30 : 1.0_dp, gs_coeffs, work_fm, 0.0_dp, gs_fm)
1566 30 : CALL cp_fm_get_diag(gs_fm, gs_diag)
1567 220 : gsgs_op(i) = 2.0_dp*SUM(gs_diag)
1568 :
1569 : END DO !i
1570 10 : DEALLOCATE (gs_diag)
1571 :
1572 10 : CALL cp_fm_release(work_fm)
1573 :
1574 : ! Create the work and helper fms
1575 10 : CALL cp_fm_create(vec_op, std_struct) !nao nactive
1576 10 : 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 10 : nrow_global=nex, ncol_global=nex)
1581 : CALL cp_fm_struct_create(sggs_struct, &
1582 : context=blacs_env, para_env=para_env, &
1583 10 : nrow_global=nactive*nsg, ncol_global=nactive)
1584 :
1585 10 : CALL cp_fm_create(tmp_fm, tmp_struct)
1586 10 : CALL cp_fm_create(work_fm, full_struct)
1587 10 : CALL cp_fm_create(sggs_fm, sggs_struct)
1588 :
1589 20 : ALLOCATE (diag(nactive))
1590 40 : ALLOCATE (mo_op(nactive, nactive))
1591 30 : 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 40 : DO i = 1, dim_op
1596 :
1597 30 : ao_op_i => ao_op(i)%matrix
1598 :
1599 : ! The GS-GS contribution
1600 30 : CALL cp_fm_set_element(amew_op(i), 1, 1, gsgs_op(i))
1601 :
1602 : ! Compute the operator in the MO basis
1603 30 : CALL cp_dbcsr_sm_fm_multiply(ao_op_i, gs_coeffs, vec_op, ncol=nactive)
1604 30 : 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 30 : IF (vel_rep) THEN
1608 : CALL dip_vel_op(soc_env, qs_env, evects_sing, dbcsr_prod, i, .FALSE., &
1609 6 : gs_coeffs, sggs_fm)
1610 : ELSE
1611 : CALL parallel_gemm('T', 'N', nactive*nsg, nactive, nao, &
1612 24 : 1.0_dp, soc_env%a_coeff, vec_op, 0.0_dp, sggs_fm)
1613 : END IF
1614 :
1615 114 : 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 84 : n_cols=nactive)
1622 84 : diag(:) = get_diag(sggs_block)
1623 600 : op = sqrt2*SUM(diag)
1624 114 : 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 30 : 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 30 : filter_eps=eps_filter)
1637 :
1638 30 : IF (vel_rep) THEN
1639 6 : 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 24 : 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 24 : 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 30 : symmetric=.TRUE.)
1661 :
1662 30 : 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 30 : 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 30 : 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 30 : filter_eps=eps_filter)
1683 :
1684 : !the operator in the LR orbital basis
1685 30 : IF (vel_rep) THEN
1686 6 : 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 24 : 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 24 : 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 30 : symmetric=.TRUE.)
1706 :
1707 30 : 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 30 : 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 30 : 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 30 : t_firstcol=1 + nsg + 2*ntp + 1)
1736 :
1737 : ! Symmetrize the matrix (only upper triangle built)
1738 40 : CALL cp_fm_uplo_to_full(amew_op(i), work_fm)
1739 :
1740 : END DO !i
1741 :
1742 : ! Clean-up
1743 10 : CALL cp_fm_release(gs_fm)
1744 10 : CALL cp_fm_release(work_fm)
1745 10 : CALL cp_fm_release(tmp_fm)
1746 10 : CALL cp_fm_release(vec_op)
1747 10 : CALL cp_fm_release(sggs_fm)
1748 10 : CALL cp_fm_struct_release(full_struct)
1749 10 : CALL cp_fm_struct_release(tmp_struct)
1750 10 : CALL cp_fm_struct_release(sggs_struct)
1751 :
1752 10 : DEALLOCATE (sggs_block)
1753 10 : DEALLOCATE (diag, gsgs_op, mo_op)
1754 10 : CALL timestop(handle)
1755 :
1756 50 : END SUBROUTINE get_rcs_amew_op
1757 :
1758 0 : END MODULE qs_tddfpt2_soc
|
