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 :
8 : ! **************************************************************************************************
9 : !> \brief Utilities absorption spectroscopy using TDDFPT with SOC
10 : !> \author JRVogt (12.2023)
11 : ! **************************************************************************************************
12 :
13 : MODULE qs_tddfpt2_soc_utils
14 : USE cp_blacs_env, ONLY: cp_blacs_env_type
15 : USE cp_cfm_types, ONLY: cp_cfm_get_info,&
16 : cp_cfm_get_submatrix,&
17 : cp_cfm_type
18 : USE cp_control_types, ONLY: tddfpt2_control_type
19 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
20 : dbcsr_create,&
21 : dbcsr_desymmetrize,&
22 : dbcsr_get_info,&
23 : dbcsr_p_type,&
24 : dbcsr_release,&
25 : dbcsr_type
26 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
27 : copy_fm_to_dbcsr,&
28 : cp_dbcsr_sm_fm_multiply,&
29 : dbcsr_allocate_matrix_set,&
30 : dbcsr_deallocate_matrix_set
31 : USE cp_fm_basic_linalg, ONLY: cp_fm_schur_product
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: cp_fm_create,&
36 : cp_fm_get_info,&
37 : cp_fm_release,&
38 : cp_fm_set_all,&
39 : cp_fm_to_fm,&
40 : cp_fm_to_fm_submat,&
41 : cp_fm_type
42 : USE input_constants, ONLY: tddfpt_dipole_berry,&
43 : tddfpt_dipole_length,&
44 : tddfpt_dipole_velocity
45 : USE kinds, ONLY: dp
46 : USE message_passing, ONLY: mp_para_env_type
47 : USE moments_utils, ONLY: get_reference_point
48 : USE parallel_gemm_api, ONLY: parallel_gemm
49 : USE qs_environment_types, ONLY: get_qs_env,&
50 : qs_environment_type
51 : USE qs_ks_types, ONLY: qs_ks_env_type
52 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
53 : USE qs_operators_ao, ONLY: p_xyz_ao,&
54 : rRc_xyz_ao
55 : USE qs_overlap, ONLY: build_overlap_matrix
56 : USE qs_tddfpt2_soc_types, ONLY: soc_env_type
57 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos
58 :
59 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
60 : #include "./base/base_uses.f90"
61 :
62 : IMPLICIT NONE
63 : PRIVATE
64 :
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_soc_utils'
66 :
67 : PUBLIC :: soc_dipole_operator, soc_contract_evect, resort_evects, dip_vel_op
68 :
69 : !A helper type for SOC
70 : TYPE dbcsr_soc_package_type
71 : TYPE(dbcsr_type), POINTER :: dbcsr_sg => Null()
72 : TYPE(dbcsr_type), POINTER :: dbcsr_tp => Null()
73 : TYPE(dbcsr_type), POINTER :: dbcsr_sc => Null()
74 : TYPE(dbcsr_type), POINTER :: dbcsr_sf => Null()
75 : TYPE(dbcsr_type), POINTER :: dbcsr_prod => Null()
76 : TYPE(dbcsr_type), POINTER :: dbcsr_ovlp => Null()
77 : TYPE(dbcsr_type), POINTER :: dbcsr_tmp => Null()
78 : TYPE(dbcsr_type), POINTER :: dbcsr_work => Null()
79 : END TYPE dbcsr_soc_package_type
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief Build the atomic dipole operator
85 : !> \param soc_env ...
86 : !> \param tddfpt_control informations on how to build the operaot
87 : !> \param qs_env Qucikstep environment
88 : !> \param gs_mos ...
89 : ! **************************************************************************************************
90 10 : SUBROUTINE soc_dipole_operator(soc_env, tddfpt_control, qs_env, gs_mos)
91 : TYPE(soc_env_type), TARGET :: soc_env
92 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
93 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
94 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
95 : INTENT(in) :: gs_mos
96 :
97 : CHARACTER(len=*), PARAMETER :: routineN = 'soc_dipole_operator'
98 :
99 : INTEGER :: dim_op, handle, i_dim, nao, nspin
100 : REAL(kind=dp), DIMENSION(3) :: reference_point
101 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
102 :
103 10 : CALL timeset(routineN, handle)
104 :
105 10 : NULLIFY (matrix_s)
106 :
107 10 : IF (tddfpt_control%dipole_form == tddfpt_dipole_berry) THEN
108 0 : CPABORT("BERRY DIPOLE FORM NOT IMPLEMENTED FOR SOC")
109 : END IF
110 : !! ONLY RCS have been implemented, Therefore, nspin sould always be 1!
111 10 : nspin = 1
112 : !! Number of dimensions should be 3, unless multipole is implemented in the future
113 10 : dim_op = 3
114 :
115 : !! Initzilize the dipmat structure
116 10 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
117 10 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
118 :
119 40 : ALLOCATE (soc_env%dipmat_ao(dim_op))
120 40 : DO i_dim = 1, dim_op
121 30 : ALLOCATE (soc_env%dipmat_ao(i_dim)%matrix)
122 : CALL dbcsr_copy(soc_env%dipmat_ao(i_dim)%matrix, &
123 : matrix_s(1)%matrix, &
124 40 : name="dipole operator matrix")
125 : END DO
126 :
127 18 : SELECT CASE (tddfpt_control%dipole_form)
128 : CASE (tddfpt_dipole_length)
129 : !!This routine is analog to qs_tddfpt_prperties but only until the rRc_xyz_ao routine
130 : !! This will lead to an operator within the nao x nao basis
131 : !! qs_tddpft_properies uses nvirt x nocc
132 : CALL get_reference_point(reference_point, qs_env=qs_env, &
133 : reference=tddfpt_control%dipole_reference, &
134 8 : ref_point=tddfpt_control%dipole_ref_point)
135 :
136 : CALL rRc_xyz_ao(op=soc_env%dipmat_ao, qs_env=qs_env, rc=reference_point, order=1, &
137 8 : minimum_image=.FALSE., soft=.FALSE.)
138 : !! This will lead to S C^virt C^virt,T Q_q (vgl Strand et al., J. Chem Phys. 150, 044702, 2019)
139 8 : CALL length_rep(qs_env, gs_mos, soc_env)
140 : CASE (tddfpt_dipole_velocity)
141 : !!This Routine calcluates the dipole Operator within the velocity-form within the ao basis
142 : !!This Operation is only used in xas_tdp and qs_tddfpt_soc, lines uses rmc_x_p_xyz_ao
143 2 : CALL p_xyz_ao(soc_env%dipmat_ao, qs_env, minimum_image=.FALSE.)
144 : !! This will precomute SC^virt, (omega^a-omega^i)^-1 and C^virt dS/dq
145 2 : CALL velocity_rep(qs_env, gs_mos, soc_env)
146 : CASE DEFAULT
147 10 : CPABORT("Unimplemented form of the dipole operator")
148 : END SELECT
149 :
150 10 : CALL timestop(handle)
151 :
152 20 : END SUBROUTINE soc_dipole_operator
153 :
154 : ! **************************************************************************************************
155 : !> \brief ...
156 : !> \param qs_env ...
157 : !> \param gs_mos ...
158 : !> \param soc_env ...
159 : ! **************************************************************************************************
160 8 : SUBROUTINE length_rep(qs_env, gs_mos, soc_env)
161 : TYPE(qs_environment_type), POINTER :: qs_env
162 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
163 : INTENT(in) :: gs_mos
164 : TYPE(soc_env_type), TARGET :: soc_env
165 :
166 : INTEGER :: ideriv, ispin, nao, nderivs, nspins
167 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nmo_virt
168 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
169 : TYPE(cp_fm_struct_type), POINTER :: dip_struct, fm_struct
170 8 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: S_mos_virt
171 8 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: dipole_op_mos_occ
172 : TYPE(cp_fm_type), POINTER :: dipmat_tmp, wfm_ao_ao
173 8 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
174 : TYPE(dbcsr_type), POINTER :: symm_tmp
175 : TYPE(mp_para_env_type), POINTER :: para_env
176 :
177 8 : CALL get_qs_env(qs_env, matrix_s=matrix_s, blacs_env=blacs_env, para_env=para_env)
178 :
179 8 : nderivs = 3
180 8 : nspins = 1 !!We only account for rcs, will be changed in the future
181 8 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
182 : ALLOCATE (S_mos_virt(nspins), dipole_op_mos_occ(3, nspins), &
183 48 : wfm_ao_ao, nmo_virt(nspins), symm_tmp, dipmat_tmp)
184 :
185 8 : CALL cp_fm_struct_create(dip_struct, context=blacs_env, ncol_global=nao, nrow_global=nao, para_env=para_env)
186 :
187 8 : CALL dbcsr_allocate_matrix_set(soc_env%dipmat, nderivs)
188 8 : CALL dbcsr_desymmetrize(matrix_s(1)%matrix, symm_tmp)
189 32 : DO ideriv = 1, nderivs
190 24 : ALLOCATE (soc_env%dipmat(ideriv)%matrix)
191 : CALL dbcsr_create(soc_env%dipmat(ideriv)%matrix, template=symm_tmp, &
192 24 : name="contracted operator", matrix_type="N")
193 56 : DO ispin = 1, nspins
194 48 : CALL cp_fm_create(dipole_op_mos_occ(ideriv, ispin), matrix_struct=dip_struct)
195 : END DO
196 : END DO
197 :
198 8 : CALL dbcsr_release(symm_tmp)
199 8 : DEALLOCATE (symm_tmp)
200 :
201 16 : DO ispin = 1, nspins
202 8 : nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
203 8 : CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct)
204 8 : CALL cp_fm_create(wfm_ao_ao, dip_struct)
205 8 : CALL cp_fm_create(S_mos_virt(ispin), fm_struct)
206 :
207 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
208 : gs_mos(ispin)%mos_virt, &
209 : S_mos_virt(ispin), &
210 8 : ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
211 : CALL parallel_gemm('N', 'T', nao, nao, nmo_virt(ispin), &
212 : 1.0_dp, S_mos_virt(ispin), gs_mos(ispin)%mos_virt, &
213 8 : 0.0_dp, wfm_ao_ao)
214 :
215 32 : DO ideriv = 1, nderivs
216 24 : CALL cp_fm_create(dipmat_tmp, dip_struct)
217 24 : CALL copy_dbcsr_to_fm(soc_env%dipmat_ao(ideriv)%matrix, dipmat_tmp)
218 : CALL parallel_gemm('N', 'T', nao, nao, nao, &
219 : 1.0_dp, wfm_ao_ao, dipmat_tmp, &
220 24 : 0.0_dp, dipole_op_mos_occ(ideriv, ispin))
221 24 : CALL copy_fm_to_dbcsr(dipole_op_mos_occ(ideriv, ispin), soc_env%dipmat(ideriv)%matrix)
222 56 : CALL cp_fm_release(dipmat_tmp)
223 : END DO
224 8 : CALL cp_fm_release(wfm_ao_ao)
225 24 : DEALLOCATE (wfm_ao_ao)
226 : END DO
227 :
228 8 : CALL cp_fm_struct_release(dip_struct)
229 16 : DO ispin = 1, nspins
230 8 : CALL cp_fm_release(S_mos_virt(ispin))
231 40 : DO ideriv = 1, nderivs
232 32 : CALL cp_fm_release(dipole_op_mos_occ(ideriv, ispin))
233 : END DO
234 : END DO
235 8 : DEALLOCATE (S_mos_virt, dipole_op_mos_occ, nmo_virt, dipmat_tmp)
236 :
237 8 : END SUBROUTINE length_rep
238 :
239 : ! **************************************************************************************************
240 : !> \brief ...
241 : !> \param qs_env ...
242 : !> \param gs_mos ...
243 : !> \param soc_env ...
244 : ! **************************************************************************************************
245 2 : SUBROUTINE velocity_rep(qs_env, gs_mos, soc_env)
246 : TYPE(qs_environment_type), POINTER :: qs_env
247 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
248 : INTENT(in) :: gs_mos
249 : TYPE(soc_env_type), TARGET :: soc_env
250 :
251 : INTEGER :: ici, icol, ideriv, irow, ispin, n_act, &
252 : n_virt, nao, ncols_local, nderivs, &
253 : nrows_local, nspins
254 2 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
255 : REAL(kind=dp) :: eval_occ
256 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
257 2 : POINTER :: local_data_ediff
258 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
259 : TYPE(cp_fm_struct_type), POINTER :: ao_cvirt_struct, cvirt_ao_struct, &
260 : fm_struct, scrm_struct
261 : TYPE(cp_fm_type) :: scrm_fm
262 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, scrm
263 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
264 2 : POINTER :: sab_orb
265 : TYPE(qs_ks_env_type), POINTER :: ks_env
266 :
267 2 : NULLIFY (scrm, scrm_struct, blacs_env, matrix_s, ao_cvirt_struct, cvirt_ao_struct)
268 2 : nspins = 1
269 2 : nderivs = 3
270 18 : ALLOCATE (soc_env%SC(nspins), soc_env%CdS(nspins, nderivs), soc_env%ediff(nspins))
271 :
272 2 : CALL get_qs_env(qs_env, ks_env=ks_env, sab_orb=sab_orb, blacs_env=blacs_env, matrix_s=matrix_s)
273 2 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
274 : CALL cp_fm_struct_create(scrm_struct, nrow_global=nao, ncol_global=nao, &
275 2 : context=blacs_env)
276 2 : CALL cp_fm_get_info(gs_mos(1)%mos_virt, matrix_struct=ao_cvirt_struct)
277 :
278 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, nderivative=1, &
279 : basis_type_a="ORB", basis_type_b="ORB", &
280 2 : sab_nl=sab_orb)
281 :
282 4 : DO ispin = 1, nspins
283 2 : NULLIFY (fm_struct)
284 : !deb n_occ = SIZE(gs_mos(ispin)%evals_occ)
285 2 : n_act = gs_mos(ispin)%nmo_active
286 2 : n_virt = SIZE(gs_mos(ispin)%evals_virt)
287 : CALL cp_fm_struct_create(fm_struct, nrow_global=n_virt, &
288 2 : ncol_global=n_act, context=blacs_env)
289 : CALL cp_fm_struct_create(cvirt_ao_struct, nrow_global=n_virt, &
290 2 : ncol_global=nao, context=blacs_env)
291 2 : CALL cp_fm_create(soc_env%ediff(ispin), fm_struct)
292 2 : CALL cp_fm_create(soc_env%SC(ispin), ao_cvirt_struct)
293 :
294 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
295 : gs_mos(ispin)%mos_virt, &
296 : soc_env%SC(ispin), &
297 2 : ncol=n_virt, alpha=1.0_dp, beta=0.0_dp)
298 :
299 : CALL cp_fm_get_info(soc_env%ediff(ispin), nrow_local=nrows_local, ncol_local=ncols_local, &
300 2 : row_indices=row_indices, col_indices=col_indices, local_data=local_data_ediff)
301 :
302 : !$OMP PARALLEL DO DEFAULT(NONE), &
303 : !$OMP PRIVATE(eval_occ, ici, icol, irow), &
304 2 : !$OMP SHARED(col_indices, gs_mos, ispin, local_data_ediff, ncols_local, nrows_local, row_indices)
305 : DO icol = 1, ncols_local
306 : ! E_occ_i ; imo_occ = col_indices(icol)
307 : ici = gs_mos(ispin)%index_active(col_indices(icol))
308 : eval_occ = gs_mos(ispin)%evals_occ(ici)
309 :
310 : DO irow = 1, nrows_local
311 : ! ediff_inv_weights(a, i) = 1.0 / (E_virt_a - E_occ_i)
312 : ! imo_virt = row_indices(irow)
313 : local_data_ediff(irow, icol) = 1.0_dp/(gs_mos(ispin)%evals_virt(row_indices(irow)) - eval_occ)
314 : END DO
315 : END DO
316 : !$OMP END PARALLEL DO
317 :
318 8 : DO ideriv = 1, nderivs
319 6 : CALL cp_fm_create(soc_env%CdS(ispin, ideriv), cvirt_ao_struct)
320 6 : CALL cp_fm_create(scrm_fm, scrm_struct)
321 6 : CALL copy_dbcsr_to_fm(scrm(ideriv + 1)%matrix, scrm_fm)
322 : CALL parallel_gemm('T', 'N', n_virt, nao, nao, 1.0_dp, gs_mos(ispin)%mos_virt, &
323 6 : scrm_fm, 0.0_dp, soc_env%CdS(ispin, ideriv))
324 14 : CALL cp_fm_release(scrm_fm)
325 :
326 : END DO
327 :
328 6 : CALL cp_fm_struct_release(fm_struct)
329 : END DO
330 2 : CALL dbcsr_deallocate_matrix_set(scrm)
331 2 : CALL cp_fm_struct_release(scrm_struct)
332 2 : CALL cp_fm_struct_release(cvirt_ao_struct)
333 :
334 4 : END SUBROUTINE velocity_rep
335 :
336 : ! **************************************************************************************************
337 : !> \brief This routine will construct the dipol operator within velocity representation
338 : !> \param soc_env ..
339 : !> \param qs_env ...
340 : !> \param evec_fm ...
341 : !> \param op ...
342 : !> \param ideriv ...
343 : !> \param tp ...
344 : !> \param gs_coeffs ...
345 : !> \param sggs_fm ...
346 : ! **************************************************************************************************
347 18 : SUBROUTINE dip_vel_op(soc_env, qs_env, evec_fm, op, ideriv, tp, gs_coeffs, sggs_fm)
348 : TYPE(soc_env_type), TARGET :: soc_env
349 : TYPE(qs_environment_type), POINTER :: qs_env
350 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evec_fm
351 : TYPE(dbcsr_type), INTENT(INOUT) :: op
352 : INTEGER, INTENT(IN) :: ideriv
353 : LOGICAL, INTENT(IN) :: tp
354 : TYPE(cp_fm_type), OPTIONAL, POINTER :: gs_coeffs
355 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: sggs_fm
356 :
357 : INTEGER :: iex, ispin, n_act, n_virt, nao, nex
358 : LOGICAL :: sggs
359 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
360 : TYPE(cp_fm_struct_type), POINTER :: op_struct, virt_occ_struct
361 : TYPE(cp_fm_type) :: CdSC, op_fm, SCWCdSC, WCdSC
362 18 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: WCdSC_tmp
363 : TYPE(cp_fm_type), POINTER :: coeff
364 : TYPE(mp_para_env_type), POINTER :: para_env
365 :
366 18 : NULLIFY (virt_occ_struct, virt_occ_struct, op_struct, blacs_env, para_env, coeff)
367 :
368 18 : IF (tp) THEN
369 6 : coeff => soc_env%b_coeff
370 : ELSE
371 12 : coeff => soc_env%a_coeff
372 : END IF
373 :
374 18 : sggs = .FALSE.
375 18 : IF (PRESENT(gs_coeffs)) sggs = .TRUE.
376 :
377 18 : ispin = 1 !! only rcs availble
378 18 : nex = SIZE(evec_fm, 2)
379 90 : IF (.NOT. sggs) ALLOCATE (WCdSC_tmp(ispin, nex))
380 18 : CALL get_qs_env(qs_env, blacs_env=blacs_env, para_env=para_env)
381 18 : CALL cp_fm_get_info(soc_env%CdS(ispin, ideriv), ncol_global=nao, nrow_global=n_virt)
382 18 : CALL cp_fm_get_info(evec_fm(1, 1), ncol_global=n_act)
383 :
384 18 : IF (sggs) THEN
385 : CALL cp_fm_struct_create(virt_occ_struct, context=blacs_env, para_env=para_env, nrow_global=n_virt, &
386 6 : ncol_global=n_act)
387 : CALL cp_fm_struct_create(op_struct, context=blacs_env, para_env=para_env, nrow_global=n_act*nex, &
388 6 : ncol_global=n_act)
389 : ELSE
390 : CALL cp_fm_struct_create(virt_occ_struct, context=blacs_env, para_env=para_env, nrow_global=n_virt, &
391 12 : ncol_global=n_act*nex)
392 : CALL cp_fm_struct_create(op_struct, context=blacs_env, para_env=para_env, nrow_global=n_act*nex, &
393 12 : ncol_global=n_act*nex)
394 : END IF
395 :
396 18 : CALL cp_fm_create(CdSC, soc_env%ediff(ispin)%matrix_struct)
397 18 : CALL cp_fm_create(op_fm, op_struct)
398 :
399 18 : IF (sggs) THEN
400 6 : CALL cp_fm_create(SCWCdSC, gs_coeffs%matrix_struct)
401 6 : CALL cp_fm_create(WCdSC, soc_env%ediff(ispin)%matrix_struct)
402 : CALL parallel_gemm('N', 'N', n_virt, n_act, nao, 1.0_dp, soc_env%CdS(ispin, ideriv), &
403 6 : gs_coeffs, 0.0_dp, CdSC)
404 6 : CALL cp_fm_schur_product(CdSC, soc_env%ediff(ispin), WCdSC)
405 : ELSE
406 12 : CALL cp_fm_create(SCWCdSC, coeff%matrix_struct)
407 36 : DO iex = 1, nex
408 24 : CALL cp_fm_create(WCdSC_tmp(ispin, iex), soc_env%ediff(ispin)%matrix_struct)
409 : CALL parallel_gemm('N', 'N', n_virt, n_act, nao, 1.0_dp, soc_env%CdS(ispin, ideriv), &
410 24 : evec_fm(ispin, iex), 0.0_dp, CdSC)
411 36 : CALL cp_fm_schur_product(CdSC, soc_env%ediff(ispin), WCdSC_tmp(ispin, iex))
412 : END DO
413 12 : CALL cp_fm_create(WCdSC, virt_occ_struct)
414 12 : CALL soc_contract_evect(WCdSC_tmp, WCdSC)
415 36 : DO iex = 1, nex
416 36 : CALL cp_fm_release(WCdSC_tmp(ispin, iex))
417 : END DO
418 12 : DEALLOCATE (WCdSC_tmp)
419 : END IF
420 :
421 18 : IF (sggs) THEN
422 6 : CALL parallel_gemm('N', 'N', nao, n_act, n_virt, 1.0_dp, soc_env%SC(ispin), WCdSC, 0.0_dp, SCWCdSC)
423 6 : CALL parallel_gemm('T', 'N', n_act*nex, n_act, nao, 1.0_dp, soc_env%a_coeff, SCWCdSC, 0.0_dp, op_fm)
424 : ELSE
425 12 : CALL parallel_gemm('N', 'N', nao, n_act*nex, n_virt, 1.0_dp, soc_env%SC(ispin), WCdSC, 0.0_dp, SCWCdSC)
426 12 : CALL parallel_gemm('T', 'N', n_act*nex, n_act*nex, nao, 1.0_dp, coeff, SCWCdSC, 0.0_dp, op_fm)
427 : END IF
428 :
429 18 : IF (sggs) THEN
430 6 : CALL cp_fm_to_fm(op_fm, sggs_fm)
431 : ELSE
432 12 : CALL copy_fm_to_dbcsr(op_fm, op)
433 : END IF
434 :
435 18 : CALL cp_fm_release(op_fm)
436 18 : CALL cp_fm_release(WCdSC)
437 18 : CALL cp_fm_release(SCWCdSC)
438 18 : CALL cp_fm_release(CdSC)
439 18 : CALL cp_fm_struct_release(virt_occ_struct)
440 18 : CALL cp_fm_struct_release(op_struct)
441 :
442 36 : END SUBROUTINE dip_vel_op
443 :
444 : ! **************************************************************************************************
445 : !> \brief ...
446 : !> \param fm_start ...
447 : !> \param fm_res ...
448 : ! **************************************************************************************************
449 32 : SUBROUTINE soc_contract_evect(fm_start, fm_res)
450 :
451 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: fm_start
452 : TYPE(cp_fm_type), INTENT(inout) :: fm_res
453 :
454 : CHARACTER(len=*), PARAMETER :: routineN = 'soc_contract_evect'
455 :
456 : INTEGER :: handle, ii, jj, nactive, nao, nspins, &
457 : nstates, ntmp1, ntmp2
458 :
459 32 : CALL timeset(routineN, handle)
460 :
461 32 : nstates = SIZE(fm_start, 2)
462 32 : nspins = SIZE(fm_start, 1)
463 :
464 32 : CALL cp_fm_set_all(fm_res, 0.0_dp)
465 : !! Evects are written into one matrix.
466 112 : DO ii = 1, nstates
467 192 : DO jj = 1, nspins
468 80 : CALL cp_fm_get_info(fm_start(jj, ii), nrow_global=nao, ncol_global=nactive)
469 80 : CALL cp_fm_get_info(fm_res, nrow_global=ntmp1, ncol_global=ntmp2)
470 : CALL cp_fm_to_fm_submat(fm_start(jj, ii), &
471 : fm_res, &
472 : nao, nactive, &
473 : 1, 1, 1, &
474 160 : 1 + nactive*(ii - 1) + (jj - 1)*nao*nstates)
475 : END DO !nspins
476 : END DO !nsstates
477 :
478 32 : CALL timestop(handle)
479 :
480 32 : END SUBROUTINE soc_contract_evect
481 :
482 : ! **************************************************************************************************
483 : !> \brief ...
484 : !> \param vec ...
485 : !> \param new_entry ...
486 : !> \param res ...
487 : !> \param res_int ...
488 : ! **************************************************************************************************
489 462 : SUBROUTINE test_repetition(vec, new_entry, res, res_int)
490 : INTEGER, DIMENSION(:), INTENT(IN) :: vec
491 : INTEGER, INTENT(IN) :: new_entry
492 : LOGICAL, INTENT(OUT) :: res
493 : INTEGER, INTENT(OUT), OPTIONAL :: res_int
494 :
495 : INTEGER :: i
496 :
497 462 : res = .TRUE.
498 462 : IF (PRESENT(res_int)) res_int = -1
499 :
500 4262 : DO i = 1, SIZE(vec)
501 4262 : IF (vec(i) == new_entry) THEN
502 198 : res = .FALSE.
503 198 : IF (PRESENT(res_int)) res_int = i
504 : EXIT
505 : END IF
506 : END DO
507 :
508 462 : END SUBROUTINE test_repetition
509 :
510 : ! **************************************************************************************************
511 : !> \brief Used to find out, which state has which spin-multiplicity
512 : !> \param evects_cfm ...
513 : !> \param sort ...
514 : ! **************************************************************************************************
515 8 : SUBROUTINE resort_evects(evects_cfm, sort)
516 : TYPE(cp_cfm_type), INTENT(INOUT) :: evects_cfm
517 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: sort
518 :
519 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: cpl_tmp
520 : INTEGER :: i_rep, ii, jj, ntot, tmp
521 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep_int
522 : LOGICAL :: rep
523 : REAL(dp) :: max_dev, max_wfn, wfn_sq
524 :
525 8 : CALL cp_cfm_get_info(evects_cfm, nrow_global=ntot)
526 32 : ALLOCATE (cpl_tmp(ntot, ntot))
527 32 : ALLOCATE (sort(ntot), rep_int(ntot))
528 1280 : cpl_tmp = 0_dp
529 104 : sort = 0
530 8 : max_dev = 0.5
531 8 : CALL cp_cfm_get_submatrix(evects_cfm, cpl_tmp)
532 :
533 104 : DO jj = 1, ntot
534 1272 : rep_int = 0
535 96 : tmp = 0
536 96 : max_wfn = 0_dp
537 1272 : DO ii = 1, ntot
538 1176 : wfn_sq = ABS(REAL(cpl_tmp(ii, jj)**2 - AIMAG(cpl_tmp(ii, jj)**2)))
539 1272 : IF (max_wfn <= wfn_sq) THEN
540 462 : CALL test_repetition(sort, ii, rep, rep_int(ii))
541 462 : IF (rep) THEN
542 264 : max_wfn = wfn_sq
543 264 : tmp = ii
544 : END IF
545 : END IF
546 : END DO
547 104 : IF (tmp > 0) THEN
548 96 : sort(jj) = tmp
549 : ELSE
550 0 : DO i_rep = 1, ntot
551 0 : IF (rep_int(i_rep) > 0) THEN
552 0 : max_wfn = ABS(REAL(cpl_tmp(sort(i_rep), jj)**2 - AIMAG(cpl_tmp(sort(i_rep), jj)**2))) - max_dev
553 0 : DO ii = 1, ntot
554 0 : wfn_sq = ABS(REAL(cpl_tmp(ii, jj)**2 - AIMAG(cpl_tmp(ii, jj)**2)))
555 0 : IF ((max_wfn - wfn_sq)/max_wfn <= max_dev) THEN
556 0 : CALL test_repetition(sort, ii, rep)
557 0 : IF (rep .AND. ii /= i_rep) THEN
558 0 : sort(jj) = sort(i_rep)
559 0 : sort(i_rep) = ii
560 : END IF
561 : END IF
562 : END DO
563 : END IF
564 : END DO
565 : END IF
566 : END DO
567 :
568 8 : DEALLOCATE (cpl_tmp, rep_int)
569 :
570 8 : END SUBROUTINE resort_evects
571 :
572 0 : END MODULE qs_tddfpt2_soc_utils
|
