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 Floquet stuff
10 : !> \par History
11 : !> \author Shridhar Shanbhag (27.01.2026)
12 : ! **************************************************************************************************
13 : MODULE floquet_utils
14 : USE cell_types, ONLY: cell_type,&
15 : get_cell
16 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_lu_invert,&
17 : cp_cfm_scale_and_add
18 : USE cp_cfm_types, ONLY: cp_cfm_create,&
19 : cp_cfm_get_submatrix,&
20 : cp_cfm_release,&
21 : cp_cfm_set_all,&
22 : cp_cfm_set_submatrix,&
23 : cp_cfm_type
24 : USE cp_control_types, ONLY: dft_control_type
25 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
26 : USE cp_files, ONLY: close_file,&
27 : open_file
28 : USE kinds, ONLY: default_string_length,&
29 : dp
30 : USE kpoint_k_r_trafo_simple, ONLY: replicate_rs_matrices,&
31 : rs_to_kp
32 : USE kpoint_methods, ONLY: kpoint_init_cell_index
33 : USE kpoint_types, ONLY: get_kpoint_info,&
34 : kpoint_create,&
35 : kpoint_release,&
36 : kpoint_type
37 : USE mathconstants, ONLY: gaussi,&
38 : pi,&
39 : z_one,&
40 : z_zero
41 : USE mathlib, ONLY: geeig_right,&
42 : gemm_square
43 : USE message_passing, ONLY: mp_para_env_type
44 : USE physcon, ONLY: evolt
45 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
46 : USE qs_environment_types, ONLY: get_qs_env,&
47 : qs_environment_type
48 : USE qs_mo_types, ONLY: get_mo_set,&
49 : mo_set_type
50 : USE qs_moments, ONLY: qs_moment_kpoints_deep
51 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
52 : USE util, ONLY: sort
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 :
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'floquet_utils'
60 :
61 : ! Public subroutines
62 : PUBLIC :: build_floquet_matrix, &
63 : calculate_epsilon_derivative, &
64 : build_momentum_matrix, &
65 : check_floquet_conversion, &
66 : calculate_floquet_spectral_function, &
67 : write_quasi_energy, &
68 : write_floquet_spectral_function
69 :
70 : CONTAINS
71 :
72 : ! **************************************************************************************************
73 : !> \brief ...
74 : !> \param qs_env ...
75 : !> \param xkp ...
76 : !> \param e_k ...
77 : !> \param de_dk ...
78 : ! **************************************************************************************************
79 204 : SUBROUTINE calculate_epsilon_derivative(qs_env, xkp, e_k, de_dk)
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp
82 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
83 : OPTIONAL :: e_k
84 : REAL(KIND=dp), ALLOCATABLE, &
85 : DIMENSION(:, :, :, :), OPTIONAL :: de_dk
86 :
87 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_epsilon_derivative'
88 :
89 204 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: C_dH_C, C_dS_C, C_k, dH_dk_i, dS_dk_i, &
90 204 : H_k, S_k
91 : INTEGER :: handle, i_dir, ikp, ispin, n, n_img_all, &
92 : n_spin, nao, nkp
93 204 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_all
94 204 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_all
95 : LOGICAL :: present_dedk, present_ek
96 204 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvals
97 204 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: H_rs, S_rs
98 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
99 : TYPE(cell_type), POINTER :: cell
100 204 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
101 : TYPE(kpoint_type), POINTER :: kpoints_all, kpoints_scf
102 204 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
103 : TYPE(mp_para_env_type), POINTER :: para_env
104 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
105 204 : POINTER :: sab_all
106 :
107 204 : CALL timeset(routineN, handle)
108 :
109 204 : present_ek = PRESENT(e_k) ! calculate band energies, ε_k for all kpoints xkp
110 204 : present_dedk = PRESENT(de_dk) ! calculate derivative, ∇_k ε_k of band energies for all kpoints xkp
111 :
112 204 : IF (.NOT. (present_ek .OR. present_dedk)) CPABORT("Subroutine needs either e_k or de_dk")
113 :
114 : CALL get_qs_env(qs_env, &
115 : matrix_ks_kp=matrix_ks_kp, &
116 : matrix_s_kp=matrix_s_kp, &
117 : sab_all=sab_all, &
118 : cell=cell, &
119 : kpoints=kpoints_scf, &
120 : para_env=para_env, &
121 204 : mos=mos)
122 :
123 204 : CALL get_mo_set(mo_set=mos(1), nao=nao)
124 204 : CALL get_cell(cell=cell, h=hmat)
125 :
126 204 : n_spin = SIZE(matrix_ks_kp, 1)
127 204 : nkp = SIZE(xkp, 2)
128 :
129 : ! create kpoint environment kpoints_all which contains all neighbor cells R
130 : ! without considering any lattice symmetry
131 204 : NULLIFY (kpoints_all)
132 204 : CALL kpoint_create(kpoints_all)
133 204 : CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, n_img_all)
134 204 : CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index_all, index_to_cell=index_to_cell_all)
135 :
136 2040 : ALLOCATE (S_rs(1, nao, nao, n_img_all), H_rs(n_spin, nao, nao, n_img_all), source=0.0_dp)
137 :
138 : ! Convert real-space dbcsr matrices into arrays
139 204 : CALL replicate_rs_matrices(matrix_s_kp, kpoints_scf, S_rs, cell_to_index_all)
140 204 : CALL replicate_rs_matrices(matrix_ks_kp, kpoints_scf, H_rs, cell_to_index_all)
141 :
142 214 : IF (present_dedk) ALLOCATE (de_dk(n_spin, nkp, 3, nao), source=0.0_dp)
143 1020 : IF (present_ek) ALLOCATE (e_k(n_spin, nkp, nao), source=0.0_dp)
144 :
145 : !$OMP PARALLEL DEFAULT(NONE) &
146 : !$OMP PRIVATE(ikp, ispin, S_k, H_k, eigenvals, C_k, dS_dk_i, dH_dk_i, C_dS_C, C_dH_C) &
147 : !$OMP SHARED(nao, n_spin, de_dk, e_k, present_ek, present_dedk, &
148 204 : !$OMP nkp, xkp, S_rs, H_rs, index_to_cell_all, hmat)
149 : IF (present_dedk) ALLOCATE (dS_dk_i(nao, nao), C_dS_C(nao, nao), &
150 : dH_dk_i(nao, nao), C_dH_C(nao, nao), source=z_zero)
151 : ALLOCATE (C_k(nao, nao), S_k(nao, nao), H_k(nao, nao), source=z_zero)
152 : ALLOCATE (eigenvals(nao), source=0.0_dp)
153 : !$OMP DO COLLAPSE(2)
154 : DO ispin = 1, n_spin
155 : DO ikp = 1, nkp
156 :
157 : ! S^R -> S(k), H^R -> H(k)
158 : S_k = 0
159 : H_k = 0
160 : CALL rs_to_kp(S_rs(1, :, :, :), S_k, index_to_cell_all, xkp(:, ikp))
161 : CALL rs_to_kp(H_rs(ispin, :, :, :), H_k, index_to_cell_all, xkp(:, ikp))
162 :
163 : ! Diagonalize H(k)C(k) = S(k)C(k)ε(k)
164 : CALL geeig_right(H_k, S_k, eigenvals, C_k)
165 : IF (present_ek) e_k(ispin, ikp, :) = eigenvals(:)
166 :
167 : IF (present_dedk) THEN
168 : ! Evaluate the derivatives using
169 : ! ∇ ε_k = C^H(k) ∇ H_k C(k) - ε_k C^H(k) ∇ S_k C(k)
170 : DO i_dir = 1, 3
171 : CALL rs_to_kp(S_rs(1, :, :, :), dS_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
172 : CALL rs_to_kp(H_rs(ispin, :, :, :), dH_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
173 :
174 : CALL gemm_square(C_k, 'C', dS_dk_i, 'N', C_k, 'N', C_dS_C)
175 : CALL gemm_square(C_k, 'C', dH_dk_i, 'N', C_k, 'N', C_dH_C)
176 :
177 : DO n = 1, nao
178 : de_dk(ispin, ikp, i_dir, n) = DBLE(C_dH_C(n, n)) - DBLE(eigenvals(n)*C_dS_C(n, n))
179 : END DO
180 : END DO
181 : END IF
182 : END DO
183 : END DO
184 : !$OMP END DO
185 : IF (present_dedk) DEALLOCATE (dS_dk_i, C_dS_C, dH_dk_i, C_dH_C)
186 : DEALLOCATE (S_k, H_k, C_k, eigenvals)
187 : !$OMP END PARALLEL
188 204 : DEALLOCATE (S_rs, H_rs)
189 204 : CALL kpoint_release(kpoints_all)
190 :
191 204 : CALL timestop(handle)
192 :
193 816 : END SUBROUTINE calculate_epsilon_derivative
194 :
195 : ! **************************************************************************************************
196 : !> \brief ...
197 : !> \param qs_env ...
198 : !> \param xkp ...
199 : !> \param momentum ...
200 : ! **************************************************************************************************
201 2 : SUBROUTINE build_momentum_matrix(qs_env, xkp, momentum)
202 : TYPE(qs_environment_type), POINTER :: qs_env
203 : REAL(KIND=dp), DIMENSION(3) :: xkp
204 : COMPLEX(KIND=dp), ALLOCATABLE, &
205 : DIMENSION(:, :, :, :) :: momentum
206 :
207 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_momentum_matrix'
208 :
209 : COMPLEX(KIND=dp), ALLOCATABLE, &
210 2 : DIMENSION(:, :, :, :, :) :: dipole
211 : INTEGER :: handle, i, i_dir, ispin, j, n_spin, nao
212 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp1
213 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: e_k
214 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: de_dk
215 : TYPE(dft_control_type), POINTER :: dft_control
216 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
217 :
218 2 : CALL timeset(routineN, handle)
219 :
220 : ! We calculate momentum matrix elements p_nm = <ψ_n|-iħ∇_r|ψ_m>
221 : ! p_nm = <u_n|(ħk - iħ∇_r)|u_m>
222 : ! p_nm = i d_nm (ε_n - ε_m) + ∇_k ε_n δ_nm
223 :
224 2 : NULLIFY (dft_control)
225 2 : CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
226 2 : CALL get_mo_set(mo_set=mos(1), nao=nao)
227 :
228 2 : ALLOCATE (xkp1(3, 1), source=0.0_dp)
229 8 : xkp1(:, 1) = xkp(:)
230 2 : CALL calculate_epsilon_derivative(qs_env, xkp1, e_k, de_dk)
231 2 : n_spin = dft_control%nspins
232 2 : CALL qs_moment_kpoints_deep(qs_env, xkp1, dipole)
233 :
234 : !$OMP PARALLEL DEFAULT(NONE) &
235 : !$OMP PRIVATE(ispin, i_dir, i, j) &
236 2 : !$OMP SHARED(nao, n_spin, momentum, dipole, e_k, de_dk)
237 : !$OMP DO COLLAPSE(4)
238 : DO ispin = 1, n_spin
239 : DO i_dir = 1, 3
240 : DO i = 1, nao
241 : DO j = 1, nao
242 : IF (j == i) THEN
243 : momentum(ispin, i_dir, i, j) = de_dk(ispin, 1, i_dir, i)
244 : ELSE
245 : momentum(ispin, i_dir, i, j) = gaussi*(e_k(ispin, 1, i) - e_k(ispin, 1, j))* &
246 : dipole(ispin, 1, i_dir, i, j)
247 : END IF
248 : END DO
249 : END DO
250 : END DO
251 : END DO
252 : !$OMP END DO
253 : !$OMP END PARALLEL
254 2 : DEALLOCATE (xkp1, dipole, e_k, de_dk)
255 :
256 2 : CALL timestop(handle)
257 6 : END SUBROUTINE build_momentum_matrix
258 :
259 : ! **************************************************************************************************
260 : !> \brief ...
261 : !> \param qs_env ...
262 : !> \param bs_env ...
263 : !> \param xkp ...
264 : !> \param ispin ...
265 : !> \param off_diag_m ...
266 : ! **************************************************************************************************
267 2 : SUBROUTINE build_off_diagonal_matrix(qs_env, bs_env, xkp, ispin, off_diag_m)
268 : TYPE(qs_environment_type), POINTER :: qs_env
269 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
270 : REAL(KIND=dp), DIMENSION(3) :: xkp
271 : INTEGER :: ispin
272 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: off_diag_m
273 :
274 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_off_diagonal_matrix'
275 :
276 : COMPLEX(KIND=dp), ALLOCATABLE, &
277 : DIMENSION(:, :, :, :) :: momentum
278 : COMPLEX(KIND=dp), DIMENSION(3) :: efactor
279 : INTEGER :: handle, i, i_dir, j, n_spin, nao
280 : REAL(KIND=dp) :: amplitude, omega
281 : REAL(KIND=dp), DIMENSION(3) :: e_vec, phi, polarisation
282 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
283 :
284 2 : CALL timeset(routineN, handle)
285 :
286 : ! Builds the matrix that occupies the off diagonal blocks in the Floquet Matrix H_F
287 :
288 8 : polarisation(:) = bs_env%floquet_polarisation(:)
289 8 : IF (SQRT(SUM(polarisation**2)) < EPSILON(0.0_dp)) THEN
290 0 : CPABORT("Invalid (too small) polarisation vector specified for POLARISATION_PUMP")
291 : END IF
292 :
293 2 : amplitude = bs_env%floquet_amplitude
294 8 : e_vec(:) = amplitude*polarisation
295 8 : phi(:) = pi*bs_env%floquet_phi(:)
296 2 : omega = bs_env%floquet_omega
297 2 : n_spin = bs_env%n_spin
298 :
299 2 : CALL get_qs_env(qs_env, mos=mos)
300 2 : CALL get_mo_set(mo_set=mos(1), nao=nao)
301 10 : ALLOCATE (momentum(n_spin, 3, nao, nao), source=z_zero)
302 2 : CALL build_momentum_matrix(qs_env, xkp, momentum)
303 :
304 : ! E_factor(α) = (i E(α) exp(i·φ(α)))/(2ω) where α = x,y,z
305 8 : DO i_dir = 1, 3
306 8 : efactor(i_dir) = gaussi*e_vec(i_dir)*CMPLX(DCOS(phi(i_dir)), DSIN(phi(i_dir)), KIND=dp)/(2*omega)
307 : END DO
308 :
309 8 : DO i_dir = 1, 3
310 56 : DO i = 1, nao
311 438 : DO j = 1, nao
312 : off_diag_m(i, j) = off_diag_m(i, j) + &
313 432 : momentum(ispin, i_dir, i, j)*efactor(i_dir)
314 : END DO
315 : END DO
316 : END DO
317 2 : DEALLOCATE (momentum)
318 :
319 2 : CALL timestop(handle)
320 :
321 4 : END SUBROUTINE build_off_diagonal_matrix
322 :
323 : ! **************************************************************************************************
324 : !> \brief ...
325 : !> \param qs_env ...
326 : !> \param bs_env ...
327 : !> \param xkp ...
328 : !> \param ispin ...
329 : !> \param f_index ...
330 : !> \param diag_e ...
331 : ! **************************************************************************************************
332 202 : SUBROUTINE build_diagonal_matrix(qs_env, bs_env, xkp, ispin, f_index, diag_e)
333 : TYPE(qs_environment_type), POINTER :: qs_env
334 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
335 : REAL(KIND=dp), DIMENSION(3) :: xkp
336 : INTEGER :: ispin, f_index
337 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: diag_e
338 :
339 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_diagonal_matrix'
340 :
341 : INTEGER :: handle, i, n_spin, nao
342 : REAL(KIND=dp) :: omega
343 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp1
344 202 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: e_k
345 202 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
346 :
347 202 : CALL timeset(routineN, handle)
348 :
349 : ! Builds the matrix that occupies the off diagonal blocks in the Floquet Matrix H_F
350 :
351 202 : omega = bs_env%floquet_omega
352 202 : n_spin = bs_env%n_spin
353 :
354 202 : CALL get_qs_env(qs_env, mos=mos)
355 202 : CALL get_mo_set(mo_set=mos(1), nao=nao)
356 202 : ALLOCATE (xkp1(3, 1), source=0.0_dp)
357 808 : xkp1(:, 1) = xkp(:)
358 202 : CALL calculate_epsilon_derivative(qs_env, xkp1, e_k)
359 :
360 1818 : DO i = 1, nao
361 1818 : diag_e(i, i) = e_k(ispin, 1, i) + f_index*omega
362 : END DO
363 202 : DEALLOCATE (xkp1, e_k)
364 202 : CALL timestop(handle)
365 :
366 404 : END SUBROUTINE build_diagonal_matrix
367 :
368 : ! **************************************************************************************************
369 : !> \brief ...
370 : !> \param qs_env ...
371 : !> \param bs_env ...
372 : !> \param xkp ...
373 : !> \param ispin ...
374 : !> \param floquet_matrix ...
375 : ! **************************************************************************************************
376 2 : SUBROUTINE build_floquet_matrix(qs_env, bs_env, xkp, ispin, floquet_matrix)
377 : TYPE(qs_environment_type), POINTER :: qs_env
378 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
379 : REAL(KIND=dp), DIMENSION(3) :: xkp
380 : INTEGER :: ispin
381 : TYPE(cp_cfm_type) :: floquet_matrix
382 :
383 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_floquet_matrix'
384 :
385 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: conj_off_diag_m, diag_e, off_diag_m
386 : INTEGER :: f_index, handle, i, i_f, max_f_index, &
387 : n_fbands, n_spin, nao
388 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp1
389 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
390 :
391 2 : CALL timeset(routineN, handle)
392 :
393 : CALL get_qs_env(qs_env, &
394 2 : mos=mos)
395 2 : CALL get_mo_set(mo_set=mos(1), nao=nao)
396 2 : max_f_index = bs_env%max_floquet_index
397 2 : n_fbands = 1 + 2*max_f_index
398 2 : n_spin = bs_env%n_spin
399 :
400 : ! Creates the floquet matrix H_F by placing the diagonal and off diagonal blocks
401 :
402 2 : ALLOCATE (xkp1(3, 1), source=0.0_dp)
403 8 : ALLOCATE (diag_e(nao, nao), source=z_zero)
404 6 : ALLOCATE (off_diag_m(nao, nao), source=z_zero)
405 6 : ALLOCATE (conj_off_diag_m(nao, nao), source=z_zero)
406 8 : xkp1(:, 1) = xkp(:)
407 2 : CALL build_off_diagonal_matrix(qs_env, bs_env, xkp, ispin, off_diag_m)
408 146 : conj_off_diag_m(:, :) = CONJG(TRANSPOSE(off_diag_m(:, :)))
409 :
410 654482 : floquet_matrix%local_data(:, :) = z_zero
411 204 : DO i = 1, n_fbands
412 202 : i_f = 1 + (i - 1)*nao
413 202 : f_index = i - max_f_index - 1
414 202 : CALL build_diagonal_matrix(qs_env, bs_env, xkp, ispin, f_index, diag_e)
415 202 : CALL cp_cfm_set_submatrix(floquet_matrix, diag_e, i_f, i_f)
416 204 : IF (i > 1) THEN
417 200 : CALL cp_cfm_set_submatrix(floquet_matrix, off_diag_m, i_f, i_f - nao)
418 200 : CALL cp_cfm_set_submatrix(floquet_matrix, conj_off_diag_m, i_f - nao, i_f)
419 : END IF
420 : END DO
421 2 : DEALLOCATE (diag_e, off_diag_m, conj_off_diag_m, xkp1)
422 :
423 2 : CALL timestop(handle)
424 4 : END SUBROUTINE build_floquet_matrix
425 :
426 : ! **************************************************************************************************
427 : !> \brief ...
428 : !> \param qs_env ...
429 : !> \param bs_env ...
430 : !> \param cfm_eigenvectors ...
431 : ! **************************************************************************************************
432 2 : SUBROUTINE check_floquet_conversion(qs_env, bs_env, cfm_eigenvectors)
433 : TYPE(qs_environment_type), POINTER :: qs_env
434 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
435 : TYPE(cp_cfm_type) :: cfm_eigenvectors
436 :
437 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_floquet_conversion'
438 :
439 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vector_block_0, vector_block_1
440 : INTEGER :: handle, max_f_index, nao, start_row_0, &
441 : start_row_1
442 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
443 :
444 2 : CALL timeset(routineN, handle)
445 :
446 2 : CALL get_qs_env(qs_env, mos=mos)
447 2 : CALL get_mo_set(mo_set=mos(1), nao=nao)
448 :
449 2 : max_f_index = bs_env%max_floquet_index
450 :
451 2 : start_row_0 = 1 + nao*max_f_index
452 2 : start_row_1 = 1 + nao*(max_f_index + 1)
453 :
454 12 : ALLOCATE (vector_block_0(3*nao, nao), vector_block_1(3*nao, nao), source=z_zero)
455 2 : CALL cp_cfm_get_submatrix(cfm_eigenvectors, vector_block_0, start_row_0, start_row_0)
456 2 : CALL cp_cfm_get_submatrix(cfm_eigenvectors, vector_block_1, start_row_1, start_row_1)
457 2 : IF (bs_env%eps_floquet > 0) THEN
458 802 : IF (ABS(MAXVAL(ABS(vector_block_0)) - MAXVAL(ABS(vector_block_1))) > bs_env%eps_floquet) THEN
459 0 : CPABORT("FLOQUET DIAGONALIZATION DID NOT CONVERGE. MAX_FLOQUET_INDEX TOO SMALL")
460 : END IF
461 : END IF
462 2 : DEALLOCATE (vector_block_0, vector_block_1)
463 2 : CALL timestop(handle)
464 :
465 4 : END SUBROUTINE check_floquet_conversion
466 :
467 : ! **************************************************************************************************
468 : !> \brief ...
469 : !> \param qs_env ...
470 : !> \param bs_env ...
471 : !> \param floquet_matrix ...
472 : !> \param a_k ...
473 : ! **************************************************************************************************
474 2 : SUBROUTINE calculate_floquet_spectral_function(qs_env, bs_env, floquet_matrix, a_k)
475 : TYPE(qs_environment_type), POINTER :: qs_env
476 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
477 : TYPE(cp_cfm_type) :: floquet_matrix
478 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: a_k
479 :
480 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_floquet_spectral_function'
481 :
482 : COMPLEX(KIND=dp) :: diag
483 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: g00
484 : INTEGER :: handle, i, i_E, max_f_index, n_E, nao, &
485 : start_row_g00
486 : REAL(KIND=dp) :: broad, E_min, energy, energy_step, mu
487 : TYPE(cp_cfm_type) :: g_r
488 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
489 :
490 2 : CALL timeset(routineN, handle)
491 :
492 2 : CALL get_qs_env(qs_env, mos=mos)
493 :
494 2 : CALL get_mo_set(mo_set=mos(1), mu=mu, nao=nao)
495 :
496 2 : broad = bs_env%broadening_floquet
497 2 : energy_step = bs_env%energy_step_floquet
498 2 : E_min = mu - bs_env%energy_window_floquet
499 2 : max_f_index = bs_env%max_floquet_index
500 2 : n_E = SIZE(a_k)
501 2 : start_row_g00 = 1 + nao*max_f_index
502 :
503 8 : ALLOCATE (g00(nao, nao))
504 2 : CALL cp_cfm_create(g_r, floquet_matrix%matrix_struct, set_zero=.TRUE.)
505 :
506 10 : DO i_E = 1, n_E
507 8 : energy = E_min + i_E*energy_step
508 8 : diag = energy + gaussi*broad/2.0_dp
509 8 : CALL cp_cfm_set_all(g_r, z_zero, diag)
510 8 : CALL cp_cfm_scale_and_add(z_one, g_r, -z_one, floquet_matrix)
511 8 : CALL cp_cfm_lu_invert(g_r)
512 8 : g00(:, :) = 0.0_dp
513 8 : CALL cp_cfm_get_submatrix(g_r, g00, start_row_g00, start_row_g00)
514 146 : a_k(i_E) = -AIMAG(SUM([(g00(i, i), i=1, nao)]))/pi
515 : END DO
516 :
517 2 : DEALLOCATE (g00)
518 2 : CALL cp_cfm_release(g_r)
519 :
520 2 : CALL timestop(handle)
521 :
522 4 : END SUBROUTINE calculate_floquet_spectral_function
523 :
524 : ! **************************************************************************************************
525 : !> \brief ...
526 : !> \param qs_env ...
527 : !> \param bs_env ...
528 : !> \param ispin ...
529 : !> \param ikp_for_file ...
530 : !> \param xkp ...
531 : !> \param eigenvalues ...
532 : ! **************************************************************************************************
533 2 : SUBROUTINE write_quasi_energy(qs_env, bs_env, ispin, ikp_for_file, xkp, eigenvalues)
534 : TYPE(qs_environment_type), POINTER :: qs_env
535 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
536 : INTEGER :: ispin, ikp_for_file
537 : REAL(KIND=dp), DIMENSION(3) :: xkp
538 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
539 :
540 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_quasi_energy'
541 :
542 : CHARACTER(LEN=default_string_length) :: fname
543 : INTEGER :: handle, i, j, n, nao, nkp_start, qunit
544 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: indices, work
545 : REAL(KIND=dp) :: fbz_max, fbz_min, omega, qe
546 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: quasi_energies
547 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
548 :
549 2 : CALL timeset(routineN, handle)
550 :
551 2 : CALL get_qs_env(qs_env, mos=mos)
552 2 : CALL get_mo_set(mo_set=mos(1), nao=nao)
553 :
554 2 : nkp_start = bs_env%nkp_only_DOS
555 2 : omega = bs_env%floquet_omega
556 2 : fbz_max = omega/2.0_dp
557 2 : fbz_min = -omega/2.0_dp
558 :
559 2 : n = SIZE(eigenvalues)
560 12 : ALLOCATE (indices(nao), work(nao), quasi_energies(nao))
561 34 : indices(:) = [(i, i=1, nao)]
562 18 : DO i = 1, nao
563 12946 : DO j = 1, n
564 115776 : IF (ANY(indices == j)) CYCLE
565 12816 : IF (ABS(eigenvalues(j)) < ABS(eigenvalues(indices(i)))) indices(i) = j
566 : END DO
567 : END DO
568 2 : CALL sort(indices, nao, work)
569 :
570 18 : quasi_energies(:) = eigenvalues(indices(:))
571 :
572 2 : IF (bs_env%para_env%is_source()) THEN
573 1 : WRITE (fname, "(2A)") TRIM(bs_env%floquet_qe_file), ".bs"
574 1 : IF (ikp_for_file == 1) THEN
575 : CALL open_file(TRIM(fname), unit_number=qunit, file_status="REPLACE", &
576 1 : file_action="WRITE")
577 : ELSE
578 : CALL open_file(TRIM(fname), unit_number=qunit, file_status="OLD", &
579 0 : file_action="WRITE", file_position="APPEND")
580 : END IF
581 :
582 1 : WRITE (qunit, "(A)") "# Quasi-energies obtained by diagonalising the Floquet Hamiltonian"
583 1 : WRITE (qunit, "(A)") "# (in units of eV, folded to Floquet Brillouin zone)"
584 : WRITE (qunit, "(A,I0,T10,A,I0,A,T24,3(1X,F14.8))") &
585 1 : "# Spin ", ispin, " Point ", ikp_for_file, ": ", xkp(1:3)
586 1 : WRITE (qunit, "(A)") "# Floquet band Quasi-energy [eV]"
587 9 : DO i = 1, nao
588 8 : qe = quasi_energies(i)
589 9 : WRITE (qunit, "(I8,F21.8)") i, qe*evolt
590 : END DO
591 1 : CALL close_file(qunit)
592 : END IF
593 2 : DEALLOCATE (indices, work, quasi_energies)
594 :
595 2 : CALL timestop(handle)
596 :
597 4 : END SUBROUTINE write_quasi_energy
598 : ! **************************************************************************************************
599 : !> \brief ...
600 : !> \param qs_env ...
601 : !> \param bs_env ...
602 : !> \param ispin ...
603 : !> \param ikp_for_file ...
604 : !> \param xkp ...
605 : !> \param a_k ...
606 : ! **************************************************************************************************
607 4 : SUBROUTINE write_floquet_spectral_function(qs_env, bs_env, ispin, ikp_for_file, xkp, a_k)
608 : TYPE(qs_environment_type), POINTER :: qs_env
609 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
610 : INTEGER :: ispin, ikp_for_file
611 : REAL(KIND=dp), DIMENSION(3) :: xkp
612 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: a_k
613 :
614 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_floquet_spectral_function'
615 :
616 : CHARACTER(LEN=default_string_length) :: fname
617 : INTEGER :: handle, i_E, n_E, nkp_start, wunit
618 : REAL(KIND=dp) :: E_min, energy, energy_step, mu
619 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
620 :
621 2 : CALL timeset(routineN, handle)
622 :
623 2 : CALL get_qs_env(qs_env, mos=mos)
624 :
625 2 : CALL get_mo_set(mo_set=mos(1), mu=mu)
626 :
627 2 : energy_step = bs_env%energy_step_floquet
628 2 : E_min = mu - bs_env%energy_window_floquet
629 2 : n_E = SIZE(a_k)
630 2 : nkp_start = bs_env%nkp_only_DOS
631 :
632 2 : IF (bs_env%para_env%is_source()) THEN
633 1 : WRITE (fname, "(2A)") TRIM(bs_env%floquet_dos_file), ".out"
634 1 : IF (ikp_for_file == 1) THEN
635 : CALL open_file(TRIM(fname), unit_number=wunit, file_status="REPLACE", &
636 1 : file_action="WRITE")
637 : ELSE
638 : CALL open_file(TRIM(fname), unit_number=wunit, file_status="OLD", &
639 0 : file_action="WRITE", file_position="APPEND")
640 : END IF
641 :
642 1 : WRITE (wunit, "(A)") "# Floquet Density of States: D(ω,k) = -1/π*Im[Tr_KS(G^R(ω,k))]"
643 : WRITE (wunit, "(A,I0,T10,A,I0,A,T24,3(1X,F14.8))") &
644 1 : "# Spin ", ispin, " Point ", ikp_for_file, ": ", xkp(1:3)
645 1 : WRITE (wunit, "(A)") "#Energy-E_F (eV) A(ω,k) = DOS (1/eV)"
646 5 : DO i_E = 1, n_E
647 4 : energy = E_min + i_E*energy_step
648 5 : WRITE (wunit, "(2X,2G13.4)") energy*evolt, a_k(i_E)/evolt
649 : END DO
650 1 : CALL close_file(wunit)
651 : END IF
652 :
653 2 : CALL timestop(handle)
654 :
655 2 : END SUBROUTINE write_floquet_spectral_function
656 :
657 : END MODULE floquet_utils
|
