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