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