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