Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
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_moments, ONLY: build_local_moment_matrix
53 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
54 : USE qs_operators_ao, ONLY: build_lin_mom_matrix
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 : CALL get_reference_point(reference_point, qs_env=qs_env, &
130 : reference=tddfpt_control%dipole_reference, &
131 8 : ref_point=tddfpt_control%dipole_ref_point)
132 :
133 : CALL build_local_moment_matrix(qs_env, soc_env%dipmat_ao, 1, &
134 8 : ref_point=reference_point, all_images=.TRUE.)
135 : !! This will lead to S C^virt C^virt,T Q_q (vgl Strand et al., J. Chem Phys. 150, 044702, 2019)
136 8 : CALL length_rep(qs_env, gs_mos, soc_env)
137 : CASE (tddfpt_dipole_velocity)
138 : !!This Routine calcluates the dipole Operator within the velocity-form within the ao basis
139 : !!This operation is only used in xas_tdp and qs_tddfpt_soc.
140 2 : CALL build_lin_mom_matrix(qs_env, soc_env%dipmat_ao, minimum_image=.FALSE.)
141 : !! This will precomute SC^virt, (omega^a-omega^i)^-1 and C^virt dS/dq
142 2 : CALL velocity_rep(qs_env, gs_mos, soc_env)
143 : CASE DEFAULT
144 10 : CPABORT("Unimplemented form of the dipole operator")
145 : END SELECT
146 :
147 10 : CALL timestop(handle)
148 :
149 20 : END SUBROUTINE soc_dipole_operator
150 :
151 : ! **************************************************************************************************
152 : !> \brief ...
153 : !> \param qs_env ...
154 : !> \param gs_mos ...
155 : !> \param soc_env ...
156 : ! **************************************************************************************************
157 8 : SUBROUTINE length_rep(qs_env, gs_mos, soc_env)
158 : TYPE(qs_environment_type), POINTER :: qs_env
159 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
160 : INTENT(in) :: gs_mos
161 : TYPE(soc_env_type), TARGET :: soc_env
162 :
163 : INTEGER :: ideriv, ispin, nao, nderivs, nspins
164 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nmo_virt
165 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
166 : TYPE(cp_fm_struct_type), POINTER :: dip_struct, fm_struct
167 8 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: S_mos_virt
168 8 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: dipole_op_mos_occ
169 : TYPE(cp_fm_type), POINTER :: dipmat_tmp, wfm_ao_ao
170 8 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
171 : TYPE(dbcsr_type), POINTER :: symm_tmp
172 : TYPE(mp_para_env_type), POINTER :: para_env
173 :
174 8 : CALL get_qs_env(qs_env, matrix_s=matrix_s, blacs_env=blacs_env, para_env=para_env)
175 :
176 8 : nderivs = 3
177 8 : nspins = 1 !!We only account for rcs, will be changed in the future
178 8 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
179 : ALLOCATE (S_mos_virt(nspins), dipole_op_mos_occ(3, nspins), &
180 48 : wfm_ao_ao, nmo_virt(nspins), symm_tmp, dipmat_tmp)
181 :
182 8 : CALL cp_fm_struct_create(dip_struct, context=blacs_env, ncol_global=nao, nrow_global=nao, para_env=para_env)
183 :
184 8 : CALL dbcsr_allocate_matrix_set(soc_env%dipmat, nderivs)
185 8 : CALL dbcsr_desymmetrize(matrix_s(1)%matrix, symm_tmp)
186 32 : DO ideriv = 1, nderivs
187 24 : ALLOCATE (soc_env%dipmat(ideriv)%matrix)
188 : CALL dbcsr_create(soc_env%dipmat(ideriv)%matrix, template=symm_tmp, &
189 24 : name="contracted operator", matrix_type="N")
190 56 : DO ispin = 1, nspins
191 48 : CALL cp_fm_create(dipole_op_mos_occ(ideriv, ispin), matrix_struct=dip_struct)
192 : END DO
193 : END DO
194 :
195 8 : CALL dbcsr_release(symm_tmp)
196 8 : DEALLOCATE (symm_tmp)
197 :
198 16 : DO ispin = 1, nspins
199 8 : nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
200 8 : CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct)
201 8 : CALL cp_fm_create(wfm_ao_ao, dip_struct)
202 8 : CALL cp_fm_create(S_mos_virt(ispin), fm_struct)
203 :
204 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
205 : gs_mos(ispin)%mos_virt, &
206 : S_mos_virt(ispin), &
207 8 : ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
208 : CALL parallel_gemm('N', 'T', nao, nao, nmo_virt(ispin), &
209 : 1.0_dp, S_mos_virt(ispin), gs_mos(ispin)%mos_virt, &
210 8 : 0.0_dp, wfm_ao_ao)
211 :
212 32 : DO ideriv = 1, nderivs
213 24 : CALL cp_fm_create(dipmat_tmp, dip_struct)
214 24 : CALL copy_dbcsr_to_fm(soc_env%dipmat_ao(ideriv)%matrix, dipmat_tmp)
215 : CALL parallel_gemm('N', 'T', nao, nao, nao, &
216 : 1.0_dp, wfm_ao_ao, dipmat_tmp, &
217 24 : 0.0_dp, dipole_op_mos_occ(ideriv, ispin))
218 24 : CALL copy_fm_to_dbcsr(dipole_op_mos_occ(ideriv, ispin), soc_env%dipmat(ideriv)%matrix)
219 56 : CALL cp_fm_release(dipmat_tmp)
220 : END DO
221 8 : CALL cp_fm_release(wfm_ao_ao)
222 24 : DEALLOCATE (wfm_ao_ao)
223 : END DO
224 :
225 8 : CALL cp_fm_struct_release(dip_struct)
226 16 : DO ispin = 1, nspins
227 8 : CALL cp_fm_release(S_mos_virt(ispin))
228 40 : DO ideriv = 1, nderivs
229 32 : CALL cp_fm_release(dipole_op_mos_occ(ideriv, ispin))
230 : END DO
231 : END DO
232 8 : DEALLOCATE (S_mos_virt, dipole_op_mos_occ, nmo_virt, dipmat_tmp)
233 :
234 8 : END SUBROUTINE length_rep
235 :
236 : ! **************************************************************************************************
237 : !> \brief ...
238 : !> \param qs_env ...
239 : !> \param gs_mos ...
240 : !> \param soc_env ...
241 : ! **************************************************************************************************
242 2 : SUBROUTINE velocity_rep(qs_env, gs_mos, soc_env)
243 : TYPE(qs_environment_type), POINTER :: qs_env
244 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
245 : INTENT(in) :: gs_mos
246 : TYPE(soc_env_type), TARGET :: soc_env
247 :
248 : INTEGER :: ici, icol, ideriv, irow, ispin, n_act, &
249 : n_virt, nao, ncols_local, nderivs, &
250 : nrows_local, nspins
251 2 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
252 : REAL(kind=dp) :: eval_occ
253 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
254 2 : POINTER :: local_data_ediff
255 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
256 : TYPE(cp_fm_struct_type), POINTER :: ao_cvirt_struct, cvirt_ao_struct, &
257 : fm_struct, scrm_struct
258 : TYPE(cp_fm_type) :: scrm_fm
259 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, scrm
260 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
261 2 : POINTER :: sab_orb
262 : TYPE(qs_ks_env_type), POINTER :: ks_env
263 :
264 2 : NULLIFY (scrm, scrm_struct, blacs_env, matrix_s, ao_cvirt_struct, cvirt_ao_struct)
265 2 : nspins = 1
266 2 : nderivs = 3
267 18 : ALLOCATE (soc_env%SC(nspins), soc_env%CdS(nspins, nderivs), soc_env%ediff(nspins))
268 :
269 2 : CALL get_qs_env(qs_env, ks_env=ks_env, sab_orb=sab_orb, blacs_env=blacs_env, matrix_s=matrix_s)
270 2 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
271 : CALL cp_fm_struct_create(scrm_struct, nrow_global=nao, ncol_global=nao, &
272 2 : context=blacs_env)
273 2 : CALL cp_fm_get_info(gs_mos(1)%mos_virt, matrix_struct=ao_cvirt_struct)
274 :
275 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, nderivative=1, &
276 : basis_type_a="ORB", basis_type_b="ORB", &
277 2 : sab_nl=sab_orb)
278 :
279 4 : DO ispin = 1, nspins
280 2 : NULLIFY (fm_struct)
281 : !deb n_occ = SIZE(gs_mos(ispin)%evals_occ)
282 2 : n_act = gs_mos(ispin)%nmo_active
283 2 : n_virt = SIZE(gs_mos(ispin)%evals_virt)
284 : CALL cp_fm_struct_create(fm_struct, nrow_global=n_virt, &
285 2 : ncol_global=n_act, context=blacs_env)
286 : CALL cp_fm_struct_create(cvirt_ao_struct, nrow_global=n_virt, &
287 2 : ncol_global=nao, context=blacs_env)
288 2 : CALL cp_fm_create(soc_env%ediff(ispin), fm_struct)
289 2 : CALL cp_fm_create(soc_env%SC(ispin), ao_cvirt_struct)
290 :
291 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
292 : gs_mos(ispin)%mos_virt, &
293 : soc_env%SC(ispin), &
294 2 : ncol=n_virt, alpha=1.0_dp, beta=0.0_dp)
295 :
296 : CALL cp_fm_get_info(soc_env%ediff(ispin), nrow_local=nrows_local, ncol_local=ncols_local, &
297 2 : row_indices=row_indices, col_indices=col_indices, local_data=local_data_ediff)
298 :
299 : !$OMP PARALLEL DO DEFAULT(NONE), &
300 : !$OMP PRIVATE(eval_occ, ici, icol, irow), &
301 2 : !$OMP SHARED(col_indices, gs_mos, ispin, local_data_ediff, ncols_local, nrows_local, row_indices)
302 : DO icol = 1, ncols_local
303 : ! E_occ_i ; imo_occ = col_indices(icol)
304 : ici = gs_mos(ispin)%index_active(col_indices(icol))
305 : eval_occ = gs_mos(ispin)%evals_occ(ici)
306 :
307 : DO irow = 1, nrows_local
308 : ! ediff_inv_weights(a, i) = 1.0 / (E_virt_a - E_occ_i)
309 : ! imo_virt = row_indices(irow)
310 : local_data_ediff(irow, icol) = 1.0_dp/(gs_mos(ispin)%evals_virt(row_indices(irow)) - eval_occ)
311 : END DO
312 : END DO
313 : !$OMP END PARALLEL DO
314 :
315 8 : DO ideriv = 1, nderivs
316 6 : CALL cp_fm_create(soc_env%CdS(ispin, ideriv), cvirt_ao_struct)
317 6 : CALL cp_fm_create(scrm_fm, scrm_struct)
318 6 : CALL copy_dbcsr_to_fm(scrm(ideriv + 1)%matrix, scrm_fm)
319 : CALL parallel_gemm('T', 'N', n_virt, nao, nao, 1.0_dp, gs_mos(ispin)%mos_virt, &
320 6 : scrm_fm, 0.0_dp, soc_env%CdS(ispin, ideriv))
321 14 : CALL cp_fm_release(scrm_fm)
322 :
323 : END DO
324 :
325 6 : CALL cp_fm_struct_release(fm_struct)
326 : END DO
327 2 : CALL dbcsr_deallocate_matrix_set(scrm)
328 2 : CALL cp_fm_struct_release(scrm_struct)
329 2 : CALL cp_fm_struct_release(cvirt_ao_struct)
330 :
331 4 : END SUBROUTINE velocity_rep
332 :
333 : ! **************************************************************************************************
334 : !> \brief This routine will construct the dipol operator within velocity representation
335 : !> \param soc_env ..
336 : !> \param qs_env ...
337 : !> \param evec_fm ...
338 : !> \param op ...
339 : !> \param ideriv ...
340 : !> \param tp ...
341 : !> \param gs_coeffs ...
342 : !> \param sggs_fm ...
343 : ! **************************************************************************************************
344 18 : SUBROUTINE dip_vel_op(soc_env, qs_env, evec_fm, op, ideriv, tp, gs_coeffs, sggs_fm)
345 : TYPE(soc_env_type), TARGET :: soc_env
346 : TYPE(qs_environment_type), POINTER :: qs_env
347 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evec_fm
348 : TYPE(dbcsr_type), INTENT(INOUT) :: op
349 : INTEGER, INTENT(IN) :: ideriv
350 : LOGICAL, INTENT(IN) :: tp
351 : TYPE(cp_fm_type), OPTIONAL, POINTER :: gs_coeffs
352 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: sggs_fm
353 :
354 : INTEGER :: iex, ispin, n_act, n_virt, nao, nex
355 : LOGICAL :: sggs
356 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
357 : TYPE(cp_fm_struct_type), POINTER :: op_struct, virt_occ_struct
358 : TYPE(cp_fm_type) :: CdSC, op_fm, SCWCdSC, WCdSC
359 18 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: WCdSC_tmp
360 : TYPE(cp_fm_type), POINTER :: coeff
361 : TYPE(mp_para_env_type), POINTER :: para_env
362 :
363 18 : NULLIFY (virt_occ_struct, virt_occ_struct, op_struct, blacs_env, para_env, coeff)
364 :
365 18 : IF (tp) THEN
366 6 : coeff => soc_env%b_coeff
367 : ELSE
368 12 : coeff => soc_env%a_coeff
369 : END IF
370 :
371 18 : sggs = .FALSE.
372 18 : IF (PRESENT(gs_coeffs)) sggs = .TRUE.
373 :
374 18 : ispin = 1 !! only rcs availble
375 18 : nex = SIZE(evec_fm, 2)
376 90 : IF (.NOT. sggs) ALLOCATE (WCdSC_tmp(ispin, nex))
377 18 : CALL get_qs_env(qs_env, blacs_env=blacs_env, para_env=para_env)
378 18 : CALL cp_fm_get_info(soc_env%CdS(ispin, ideriv), ncol_global=nao, nrow_global=n_virt)
379 18 : CALL cp_fm_get_info(evec_fm(1, 1), ncol_global=n_act)
380 :
381 18 : IF (sggs) THEN
382 : CALL cp_fm_struct_create(virt_occ_struct, context=blacs_env, para_env=para_env, nrow_global=n_virt, &
383 6 : ncol_global=n_act)
384 : CALL cp_fm_struct_create(op_struct, context=blacs_env, para_env=para_env, nrow_global=n_act*nex, &
385 6 : ncol_global=n_act)
386 : ELSE
387 : CALL cp_fm_struct_create(virt_occ_struct, context=blacs_env, para_env=para_env, nrow_global=n_virt, &
388 12 : ncol_global=n_act*nex)
389 : CALL cp_fm_struct_create(op_struct, context=blacs_env, para_env=para_env, nrow_global=n_act*nex, &
390 12 : ncol_global=n_act*nex)
391 : END IF
392 :
393 18 : CALL cp_fm_create(CdSC, soc_env%ediff(ispin)%matrix_struct)
394 18 : CALL cp_fm_create(op_fm, op_struct)
395 :
396 18 : IF (sggs) THEN
397 6 : CALL cp_fm_create(SCWCdSC, gs_coeffs%matrix_struct)
398 6 : CALL cp_fm_create(WCdSC, soc_env%ediff(ispin)%matrix_struct)
399 : CALL parallel_gemm('N', 'N', n_virt, n_act, nao, 1.0_dp, soc_env%CdS(ispin, ideriv), &
400 6 : gs_coeffs, 0.0_dp, CdSC)
401 6 : CALL cp_fm_schur_product(CdSC, soc_env%ediff(ispin), WCdSC)
402 : ELSE
403 12 : CALL cp_fm_create(SCWCdSC, coeff%matrix_struct)
404 36 : DO iex = 1, nex
405 24 : CALL cp_fm_create(WCdSC_tmp(ispin, iex), soc_env%ediff(ispin)%matrix_struct)
406 : CALL parallel_gemm('N', 'N', n_virt, n_act, nao, 1.0_dp, soc_env%CdS(ispin, ideriv), &
407 24 : evec_fm(ispin, iex), 0.0_dp, CdSC)
408 36 : CALL cp_fm_schur_product(CdSC, soc_env%ediff(ispin), WCdSC_tmp(ispin, iex))
409 : END DO
410 12 : CALL cp_fm_create(WCdSC, virt_occ_struct)
411 12 : CALL soc_contract_evect(WCdSC_tmp, WCdSC)
412 36 : DO iex = 1, nex
413 36 : CALL cp_fm_release(WCdSC_tmp(ispin, iex))
414 : END DO
415 12 : DEALLOCATE (WCdSC_tmp)
416 : END IF
417 :
418 18 : IF (sggs) THEN
419 6 : CALL parallel_gemm('N', 'N', nao, n_act, n_virt, 1.0_dp, soc_env%SC(ispin), WCdSC, 0.0_dp, SCWCdSC)
420 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)
421 : ELSE
422 12 : CALL parallel_gemm('N', 'N', nao, n_act*nex, n_virt, 1.0_dp, soc_env%SC(ispin), WCdSC, 0.0_dp, SCWCdSC)
423 12 : CALL parallel_gemm('T', 'N', n_act*nex, n_act*nex, nao, 1.0_dp, coeff, SCWCdSC, 0.0_dp, op_fm)
424 : END IF
425 :
426 18 : IF (sggs) THEN
427 6 : CALL cp_fm_to_fm(op_fm, sggs_fm)
428 : ELSE
429 12 : CALL copy_fm_to_dbcsr(op_fm, op)
430 : END IF
431 :
432 18 : CALL cp_fm_release(op_fm)
433 18 : CALL cp_fm_release(WCdSC)
434 18 : CALL cp_fm_release(SCWCdSC)
435 18 : CALL cp_fm_release(CdSC)
436 18 : CALL cp_fm_struct_release(virt_occ_struct)
437 18 : CALL cp_fm_struct_release(op_struct)
438 :
439 36 : END SUBROUTINE dip_vel_op
440 :
441 : ! **************************************************************************************************
442 : !> \brief ...
443 : !> \param fm_start ...
444 : !> \param fm_res ...
445 : ! **************************************************************************************************
446 32 : SUBROUTINE soc_contract_evect(fm_start, fm_res)
447 :
448 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: fm_start
449 : TYPE(cp_fm_type), INTENT(inout) :: fm_res
450 :
451 : CHARACTER(len=*), PARAMETER :: routineN = 'soc_contract_evect'
452 :
453 : INTEGER :: handle, ii, jj, nactive, nao, nspins, &
454 : nstates, ntmp1, ntmp2
455 :
456 32 : CALL timeset(routineN, handle)
457 :
458 32 : nstates = SIZE(fm_start, 2)
459 32 : nspins = SIZE(fm_start, 1)
460 :
461 32 : CALL cp_fm_set_all(fm_res, 0.0_dp)
462 : !! Evects are written into one matrix.
463 112 : DO ii = 1, nstates
464 192 : DO jj = 1, nspins
465 80 : CALL cp_fm_get_info(fm_start(jj, ii), nrow_global=nao, ncol_global=nactive)
466 80 : CALL cp_fm_get_info(fm_res, nrow_global=ntmp1, ncol_global=ntmp2)
467 : CALL cp_fm_to_fm_submat(fm_start(jj, ii), &
468 : fm_res, &
469 : nao, nactive, &
470 : 1, 1, 1, &
471 160 : 1 + nactive*(ii - 1) + (jj - 1)*nao*nstates)
472 : END DO !nspins
473 : END DO !nsstates
474 :
475 32 : CALL timestop(handle)
476 :
477 32 : END SUBROUTINE soc_contract_evect
478 :
479 : ! **************************************************************************************************
480 : !> \brief ...
481 : !> \param vec ...
482 : !> \param new_entry ...
483 : !> \param res ...
484 : !> \param res_int ...
485 : ! **************************************************************************************************
486 466 : SUBROUTINE test_repetition(vec, new_entry, res, res_int)
487 : INTEGER, DIMENSION(:), INTENT(IN) :: vec
488 : INTEGER, INTENT(IN) :: new_entry
489 : LOGICAL, INTENT(OUT) :: res
490 : INTEGER, INTENT(OUT), OPTIONAL :: res_int
491 :
492 : INTEGER :: i
493 :
494 466 : res = .TRUE.
495 466 : IF (PRESENT(res_int)) res_int = -1
496 :
497 4266 : DO i = 1, SIZE(vec)
498 4266 : IF (vec(i) == new_entry) THEN
499 204 : res = .FALSE.
500 204 : IF (PRESENT(res_int)) res_int = i
501 : EXIT
502 : END IF
503 : END DO
504 :
505 466 : END SUBROUTINE test_repetition
506 :
507 : ! **************************************************************************************************
508 : !> \brief Used to find out, which state has which spin-multiplicity
509 : !> \param evects_cfm ...
510 : !> \param sort ...
511 : ! **************************************************************************************************
512 8 : SUBROUTINE resort_evects(evects_cfm, sort)
513 : TYPE(cp_cfm_type), INTENT(INOUT) :: evects_cfm
514 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: sort
515 :
516 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: cpl_tmp
517 : INTEGER :: i_rep, ii, jj, ntot, tmp
518 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep_int
519 : LOGICAL :: rep
520 : REAL(dp) :: max_dev, max_wfn, wfn_sq
521 :
522 8 : CALL cp_cfm_get_info(evects_cfm, nrow_global=ntot)
523 32 : ALLOCATE (cpl_tmp(ntot, ntot))
524 32 : ALLOCATE (sort(ntot), rep_int(ntot))
525 8 : cpl_tmp = 0_dp
526 8 : sort = 0
527 8 : max_dev = 0.5
528 8 : CALL cp_cfm_get_submatrix(evects_cfm, cpl_tmp)
529 :
530 104 : DO jj = 1, ntot
531 96 : rep_int = 0
532 96 : tmp = 0
533 96 : max_wfn = 0_dp
534 1272 : DO ii = 1, ntot
535 1176 : wfn_sq = ABS(REAL(cpl_tmp(ii, jj)**2 - AIMAG(cpl_tmp(ii, jj)**2)))
536 1272 : IF (max_wfn <= wfn_sq) THEN
537 466 : CALL test_repetition(sort, ii, rep, rep_int(ii))
538 466 : IF (rep) THEN
539 262 : max_wfn = wfn_sq
540 262 : tmp = ii
541 : END IF
542 : END IF
543 : END DO
544 104 : IF (tmp > 0) THEN
545 96 : sort(jj) = tmp
546 : ELSE
547 0 : DO i_rep = 1, ntot
548 0 : IF (rep_int(i_rep) > 0) THEN
549 0 : max_wfn = ABS(REAL(cpl_tmp(sort(i_rep), jj)**2 - AIMAG(cpl_tmp(sort(i_rep), jj)**2))) - max_dev
550 0 : DO ii = 1, ntot
551 0 : wfn_sq = ABS(REAL(cpl_tmp(ii, jj)**2 - AIMAG(cpl_tmp(ii, jj)**2)))
552 0 : IF ((max_wfn - wfn_sq)/max_wfn <= max_dev) THEN
553 0 : CALL test_repetition(sort, ii, rep)
554 0 : IF (rep .AND. ii /= i_rep) THEN
555 0 : sort(jj) = sort(i_rep)
556 0 : sort(i_rep) = ii
557 : END IF
558 : END IF
559 : END DO
560 : END IF
561 : END DO
562 : END IF
563 : END DO
564 :
565 8 : DEALLOCATE (cpl_tmp, rep_int)
566 :
567 8 : END SUBROUTINE resort_evects
568 :
569 0 : END MODULE qs_tddfpt2_soc_utils
|
