Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
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 8 : 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 8 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
102 :
103 8 : CALL timeset(routineN, handle)
104 :
105 8 : NULLIFY (matrix_s)
106 :
107 8 : 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 8 : nspin = 1
112 : !! Number of dimensions should be 3, unless multipole is implemented in the future
113 8 : dim_op = 3
114 :
115 : !! Initzilize the dipmat structure
116 8 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
117 8 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
118 :
119 32 : ALLOCATE (soc_env%dipmat_ao(dim_op))
120 32 : DO i_dim = 1, dim_op
121 24 : 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 32 : name="dipole operator matrix")
125 : END DO
126 :
127 14 : 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 6 : 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 6 : 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 6 : 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 8 : CPABORT("Unimplemented form of the dipole operator")
148 : END SELECT
149 :
150 8 : CALL timestop(handle)
151 :
152 16 : END SUBROUTINE soc_dipole_operator
153 :
154 : ! **************************************************************************************************
155 : !> \brief ...
156 : !> \param qs_env ...
157 : !> \param gs_mos ...
158 : !> \param soc_env ...
159 : ! **************************************************************************************************
160 6 : 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 6 : 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 6 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: S_mos_virt
171 6 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: dipole_op_mos_occ
172 : TYPE(cp_fm_type), POINTER :: dipmat_tmp, wfm_ao_ao
173 6 : 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 6 : CALL get_qs_env(qs_env, matrix_s=matrix_s, blacs_env=blacs_env, para_env=para_env)
178 :
179 6 : nderivs = 3
180 6 : nspins = 1 !!We only account for rcs, will be changed in the future
181 6 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
182 : ALLOCATE (S_mos_virt(nspins), dipole_op_mos_occ(3, nspins), &
183 36 : wfm_ao_ao, nmo_virt(nspins), symm_tmp, dipmat_tmp)
184 :
185 6 : CALL cp_fm_struct_create(dip_struct, context=blacs_env, ncol_global=nao, nrow_global=nao, para_env=para_env)
186 :
187 6 : CALL dbcsr_allocate_matrix_set(soc_env%dipmat, nderivs)
188 6 : CALL dbcsr_desymmetrize(matrix_s(1)%matrix, symm_tmp)
189 24 : DO ideriv = 1, nderivs
190 18 : ALLOCATE (soc_env%dipmat(ideriv)%matrix)
191 : CALL dbcsr_create(soc_env%dipmat(ideriv)%matrix, template=symm_tmp, &
192 18 : name="contracted operator", matrix_type="N")
193 42 : DO ispin = 1, nspins
194 36 : CALL cp_fm_create(dipole_op_mos_occ(ideriv, ispin), matrix_struct=dip_struct)
195 : END DO
196 : END DO
197 :
198 6 : CALL dbcsr_release(symm_tmp)
199 6 : DEALLOCATE (symm_tmp)
200 :
201 12 : DO ispin = 1, nspins
202 6 : nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
203 6 : CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct)
204 6 : CALL cp_fm_create(wfm_ao_ao, dip_struct)
205 6 : 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 6 : 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 6 : 0.0_dp, wfm_ao_ao)
214 :
215 24 : DO ideriv = 1, nderivs
216 18 : CALL cp_fm_create(dipmat_tmp, dip_struct)
217 18 : 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 18 : 0.0_dp, dipole_op_mos_occ(ideriv, ispin))
221 18 : CALL copy_fm_to_dbcsr(dipole_op_mos_occ(ideriv, ispin), soc_env%dipmat(ideriv)%matrix)
222 42 : CALL cp_fm_release(dipmat_tmp)
223 : END DO
224 6 : CALL cp_fm_release(wfm_ao_ao)
225 18 : DEALLOCATE (wfm_ao_ao)
226 : END DO
227 :
228 6 : CALL cp_fm_struct_release(dip_struct)
229 12 : DO ispin = 1, nspins
230 6 : CALL cp_fm_release(S_mos_virt(ispin))
231 30 : DO ideriv = 1, nderivs
232 24 : CALL cp_fm_release(dipole_op_mos_occ(ideriv, ispin))
233 : END DO
234 : END DO
235 6 : DEALLOCATE (S_mos_virt, dipole_op_mos_occ, nmo_virt, dipmat_tmp)
236 :
237 6 : 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 :: icol, ideriv, irow, ispin, n_occ, &
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, local_data_wfm
258 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
259 : TYPE(cp_fm_struct_type), POINTER :: ao_cvirt_struct, ao_nocc_struct, &
260 : cvirt_ao_struct, fm_struct, scrm_struct
261 : TYPE(cp_fm_type) :: scrm_fm, wfm_mo_virt_mo_occ
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, ao_nocc_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 2 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, matrix_struct=ao_nocc_struct)
278 :
279 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, nderivative=1, &
280 : basis_type_a="ORB", basis_type_b="ORB", &
281 2 : sab_nl=sab_orb)
282 :
283 4 : DO ispin = 1, nspins
284 2 : NULLIFY (fm_struct)
285 2 : n_occ = SIZE(gs_mos(ispin)%evals_occ)
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_occ, 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(wfm_mo_virt_mo_occ, fm_struct)
293 2 : CALL cp_fm_create(soc_env%SC(ispin), ao_cvirt_struct)
294 :
295 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
296 : gs_mos(ispin)%mos_virt, &
297 : soc_env%SC(ispin), &
298 2 : ncol=n_virt, alpha=1.0_dp, beta=0.0_dp)
299 :
300 : CALL cp_fm_get_info(soc_env%ediff(ispin), nrow_local=nrows_local, ncol_local=ncols_local, &
301 2 : row_indices=row_indices, col_indices=col_indices, local_data=local_data_ediff)
302 2 : CALL cp_fm_get_info(wfm_mo_virt_mo_occ, local_data=local_data_wfm)
303 :
304 : !$OMP PARALLEL DO DEFAULT(NONE), &
305 : !$OMP PRIVATE(eval_occ, icol, irow), &
306 2 : !$OMP SHARED(col_indices, gs_mos, ispin, local_data_ediff, ncols_local, nrows_local, row_indices)
307 : DO icol = 1, ncols_local
308 : ! E_occ_i ; imo_occ = col_indices(icol)
309 : eval_occ = gs_mos(ispin)%evals_occ(col_indices(icol))
310 :
311 : DO irow = 1, nrows_local
312 : ! ediff_inv_weights(a, i) = 1.0 / (E_virt_a - E_occ_i)
313 : ! imo_virt = row_indices(irow)
314 : local_data_ediff(irow, icol) = 1.0_dp/(gs_mos(ispin)%evals_virt(row_indices(irow)) - eval_occ)
315 : END DO
316 : END DO
317 : !$OMP END PARALLEL DO
318 :
319 8 : DO ideriv = 1, nderivs
320 6 : CALL cp_fm_create(soc_env%CdS(ispin, ideriv), cvirt_ao_struct)
321 6 : CALL cp_fm_create(scrm_fm, scrm_struct)
322 6 : CALL copy_dbcsr_to_fm(scrm(ideriv + 1)%matrix, scrm_fm)
323 : CALL parallel_gemm('T', 'N', n_virt, nao, nao, 1.0_dp, gs_mos(ispin)%mos_virt, &
324 6 : scrm_fm, 0.0_dp, soc_env%CdS(ispin, ideriv))
325 14 : CALL cp_fm_release(scrm_fm)
326 :
327 : END DO
328 :
329 2 : CALL cp_fm_release(wfm_mo_virt_mo_occ)
330 8 : CALL cp_fm_struct_release(fm_struct)
331 : END DO
332 2 : CALL dbcsr_deallocate_matrix_set(scrm)
333 2 : CALL cp_fm_struct_release(scrm_struct)
334 2 : CALL cp_fm_struct_release(cvirt_ao_struct)
335 :
336 4 : END SUBROUTINE velocity_rep
337 :
338 : ! **************************************************************************************************
339 : !> \brief This routine will construct the dipol operator within velocity representation
340 : !> \param soc_env ..
341 : !> \param qs_env ...
342 : !> \param evec_fm ...
343 : !> \param op ...
344 : !> \param ideriv ...
345 : !> \param tp ...
346 : !> \param gs_coeffs ...
347 : !> \param sggs_fm ...
348 : ! **************************************************************************************************
349 18 : SUBROUTINE dip_vel_op(soc_env, qs_env, evec_fm, op, ideriv, tp, gs_coeffs, sggs_fm)
350 : TYPE(soc_env_type), TARGET :: soc_env
351 : TYPE(qs_environment_type), POINTER :: qs_env
352 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evec_fm
353 : TYPE(dbcsr_type), INTENT(INOUT) :: op
354 : INTEGER, INTENT(IN) :: ideriv
355 : LOGICAL, INTENT(IN) :: tp
356 : TYPE(cp_fm_type), OPTIONAL, POINTER :: gs_coeffs
357 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: sggs_fm
358 :
359 : INTEGER :: iex, ispin, n_occ, n_virt, nao, nex
360 : LOGICAL :: sggs
361 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
362 : TYPE(cp_fm_struct_type), POINTER :: op_struct, virt_occ_struct
363 : TYPE(cp_fm_type) :: CdSC, op_fm, SCWCdSC, WCdSC
364 18 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: WCdSC_tmp
365 : TYPE(cp_fm_type), POINTER :: coeff
366 : TYPE(mp_para_env_type), POINTER :: para_env
367 :
368 18 : NULLIFY (virt_occ_struct, virt_occ_struct, op_struct, blacs_env, para_env, coeff)
369 :
370 18 : IF (tp) THEN
371 6 : coeff => soc_env%b_coeff
372 : ELSE
373 12 : coeff => soc_env%a_coeff
374 : END IF
375 :
376 18 : sggs = .FALSE.
377 18 : IF (PRESENT(gs_coeffs)) sggs = .TRUE.
378 :
379 18 : ispin = 1 !! only rcs availble
380 18 : nex = SIZE(evec_fm, 2)
381 90 : IF (.NOT. sggs) ALLOCATE (WCdSC_tmp(ispin, nex))
382 18 : CALL get_qs_env(qs_env, blacs_env=blacs_env, para_env=para_env)
383 18 : CALL cp_fm_get_info(soc_env%CdS(ispin, ideriv), ncol_global=nao, nrow_global=n_virt)
384 18 : CALL cp_fm_get_info(evec_fm(1, 1), ncol_global=n_occ)
385 :
386 18 : IF (sggs) THEN
387 : CALL cp_fm_struct_create(virt_occ_struct, context=blacs_env, para_env=para_env, nrow_global=n_virt, &
388 6 : ncol_global=n_occ)
389 : CALL cp_fm_struct_create(op_struct, context=blacs_env, para_env=para_env, nrow_global=n_occ*nex, &
390 6 : ncol_global=n_occ)
391 : ELSE
392 : CALL cp_fm_struct_create(virt_occ_struct, context=blacs_env, para_env=para_env, nrow_global=n_virt, &
393 12 : ncol_global=n_occ*nex)
394 : CALL cp_fm_struct_create(op_struct, context=blacs_env, para_env=para_env, nrow_global=n_occ*nex, &
395 12 : ncol_global=n_occ*nex)
396 : END IF
397 :
398 18 : CALL cp_fm_create(CdSC, soc_env%ediff(ispin)%matrix_struct)
399 18 : CALL cp_fm_create(op_fm, op_struct)
400 :
401 18 : IF (sggs) THEN
402 6 : CALL cp_fm_create(SCWCdSC, gs_coeffs%matrix_struct)
403 6 : CALL cp_fm_create(WCdSC, soc_env%ediff(ispin)%matrix_struct)
404 : CALL parallel_gemm('N', 'N', n_virt, n_occ, nao, 1.0_dp, soc_env%CdS(ispin, ideriv), &
405 6 : gs_coeffs, 0.0_dp, CdSC)
406 6 : CALL cp_fm_schur_product(CdSC, soc_env%ediff(ispin), WCdSC)
407 : ELSE
408 12 : CALL cp_fm_create(SCWCdSC, coeff%matrix_struct)
409 36 : DO iex = 1, nex
410 24 : CALL cp_fm_create(WCdSC_tmp(ispin, iex), soc_env%ediff(ispin)%matrix_struct)
411 : CALL parallel_gemm('N', 'N', n_virt, n_occ, nao, 1.0_dp, soc_env%CdS(ispin, ideriv), &
412 24 : evec_fm(ispin, iex), 0.0_dp, CdSC)
413 36 : CALL cp_fm_schur_product(CdSC, soc_env%ediff(ispin), WCdSC_tmp(ispin, iex))
414 : END DO
415 12 : CALL cp_fm_create(WCdSC, virt_occ_struct)
416 12 : CALL soc_contract_evect(WCdSC_tmp, WCdSC)
417 36 : DO iex = 1, nex
418 36 : CALL cp_fm_release(WCdSC_tmp(ispin, iex))
419 : END DO
420 12 : DEALLOCATE (WCdSC_tmp)
421 : END IF
422 :
423 18 : IF (sggs) THEN
424 6 : CALL parallel_gemm('N', 'N', nao, n_occ, n_virt, 1.0_dp, soc_env%SC(ispin), WCdSC, 0.0_dp, SCWCdSC)
425 6 : CALL parallel_gemm('T', 'N', n_occ*nex, n_occ, nao, 1.0_dp, soc_env%a_coeff, SCWCdSC, 0.0_dp, op_fm)
426 : ELSE
427 12 : CALL parallel_gemm('N', 'N', nao, n_occ*nex, n_virt, 1.0_dp, soc_env%SC(ispin), WCdSC, 0.0_dp, SCWCdSC)
428 12 : CALL parallel_gemm('T', 'N', n_occ*nex, n_occ*nex, nao, 1.0_dp, coeff, SCWCdSC, 0.0_dp, op_fm)
429 : END IF
430 :
431 18 : IF (sggs) THEN
432 6 : CALL cp_fm_to_fm(op_fm, sggs_fm)
433 : ELSE
434 12 : CALL copy_fm_to_dbcsr(op_fm, op)
435 : END IF
436 :
437 18 : CALL cp_fm_release(op_fm)
438 18 : CALL cp_fm_release(WCdSC)
439 18 : CALL cp_fm_release(SCWCdSC)
440 18 : CALL cp_fm_release(CdSC)
441 18 : CALL cp_fm_struct_release(virt_occ_struct)
442 18 : CALL cp_fm_struct_release(op_struct)
443 :
444 36 : END SUBROUTINE dip_vel_op
445 :
446 : ! **************************************************************************************************
447 : !> \brief ...
448 : !> \param fm_start ...
449 : !> \param fm_res ...
450 : ! **************************************************************************************************
451 28 : SUBROUTINE soc_contract_evect(fm_start, fm_res)
452 :
453 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: fm_start
454 : TYPE(cp_fm_type), INTENT(inout) :: fm_res
455 :
456 : CHARACTER(len=*), PARAMETER :: routineN = 'soc_contract_evect'
457 :
458 : INTEGER :: handle, ii, jj, nactive, nao, nspins, &
459 : nstates, ntmp1, ntmp2
460 :
461 28 : CALL timeset(routineN, handle)
462 :
463 28 : nstates = SIZE(fm_start, 2)
464 28 : nspins = SIZE(fm_start, 1)
465 :
466 28 : CALL cp_fm_set_all(fm_res, 0.0_dp)
467 : !! Evects are written into one matrix.
468 96 : DO ii = 1, nstates
469 164 : DO jj = 1, nspins
470 68 : CALL cp_fm_get_info(fm_start(jj, ii), nrow_global=nao, ncol_global=nactive)
471 68 : CALL cp_fm_get_info(fm_res, nrow_global=ntmp1, ncol_global=ntmp2)
472 : CALL cp_fm_to_fm_submat(fm_start(jj, ii), &
473 : fm_res, &
474 : nao, nactive, &
475 : 1, 1, 1, &
476 136 : 1 + nactive*(ii - 1) + (jj - 1)*nao*nstates)
477 : END DO !nspins
478 : END DO !nsstates
479 :
480 28 : CALL timestop(handle)
481 :
482 28 : END SUBROUTINE soc_contract_evect
483 :
484 : ! **************************************************************************************************
485 : !> \brief ...
486 : !> \param vec ...
487 : !> \param new_entry ...
488 : !> \param res ...
489 : !> \param res_int ...
490 : ! **************************************************************************************************
491 478 : SUBROUTINE test_repetition(vec, new_entry, res, res_int)
492 : INTEGER, DIMENSION(:), INTENT(IN) :: vec
493 : INTEGER, INTENT(IN) :: new_entry
494 : LOGICAL, INTENT(OUT) :: res
495 : INTEGER, INTENT(OUT), OPTIONAL :: res_int
496 :
497 : INTEGER :: i
498 :
499 478 : res = .TRUE.
500 478 : IF (PRESENT(res_int)) res_int = -1
501 :
502 4226 : DO i = 1, SIZE(vec)
503 4226 : IF (vec(i) == new_entry) THEN
504 222 : res = .FALSE.
505 222 : IF (PRESENT(res_int)) res_int = i
506 : EXIT
507 : END IF
508 : END DO
509 :
510 478 : END SUBROUTINE test_repetition
511 :
512 : ! **************************************************************************************************
513 : !> \brief Used to find out, which state has which spin-multiplicity
514 : !> \param evects_cfm ...
515 : !> \param sort ...
516 : ! **************************************************************************************************
517 8 : SUBROUTINE resort_evects(evects_cfm, sort)
518 : TYPE(cp_cfm_type), INTENT(INOUT) :: evects_cfm
519 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: sort
520 :
521 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: cpl_tmp
522 : INTEGER :: i_rep, ii, jj, ntot, tmp
523 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep_int
524 : LOGICAL :: rep
525 : REAL(dp) :: max_dev, max_wfn, wfn_sq
526 :
527 8 : CALL cp_cfm_get_info(evects_cfm, nrow_global=ntot)
528 32 : ALLOCATE (cpl_tmp(ntot, ntot))
529 32 : ALLOCATE (sort(ntot), rep_int(ntot))
530 1280 : cpl_tmp = 0_dp
531 104 : sort = 0
532 8 : max_dev = 0.5
533 8 : CALL cp_cfm_get_submatrix(evects_cfm, cpl_tmp)
534 :
535 104 : DO jj = 1, ntot
536 1272 : rep_int = 0
537 96 : tmp = 0
538 96 : max_wfn = 0_dp
539 1272 : DO ii = 1, ntot
540 1176 : wfn_sq = ABS(REAL(cpl_tmp(ii, jj)**2 - AIMAG(cpl_tmp(ii, jj)**2)))
541 1272 : IF (max_wfn .LE. wfn_sq) THEN
542 478 : CALL test_repetition(sort, ii, rep, rep_int(ii))
543 478 : IF (rep) THEN
544 256 : max_wfn = wfn_sq
545 256 : tmp = ii
546 : END IF
547 : END IF
548 : END DO
549 104 : IF (tmp > 0) THEN
550 96 : sort(jj) = tmp
551 : ELSE
552 0 : DO i_rep = 1, ntot
553 0 : IF (rep_int(i_rep) > 0) THEN
554 0 : max_wfn = ABS(REAL(cpl_tmp(sort(i_rep), jj)**2 - AIMAG(cpl_tmp(sort(i_rep), jj)**2))) - max_dev
555 0 : DO ii = 1, ntot
556 0 : wfn_sq = ABS(REAL(cpl_tmp(ii, jj)**2 - AIMAG(cpl_tmp(ii, jj)**2)))
557 0 : IF ((max_wfn - wfn_sq)/max_wfn .LE. max_dev) THEN
558 0 : CALL test_repetition(sort, ii, rep)
559 0 : IF (rep .AND. ii /= i_rep) THEN
560 0 : sort(jj) = sort(i_rep)
561 0 : sort(i_rep) = ii
562 : END IF
563 : END IF
564 : END DO
565 : END IF
566 : END DO
567 : END IF
568 : END DO
569 :
570 8 : DEALLOCATE (cpl_tmp, rep_int)
571 :
572 8 : END SUBROUTINE resort_evects
573 0 : END MODULE qs_tddfpt2_soc_utils
|
